1 /*
2 SPDX-FileCopyrightText: 2001-2005 Jason Harris <jharris@30doradus.org>
3 SPDX-FileCopyrightText: 2004-2005 Pablo de Vicente <p.devicente@wanadoo.es>
4
5 SPDX-License-Identifier: GPL-2.0-or-later
6 */
7
8 #include "skypoint.h"
9
10 #include "dms.h"
11 #include "ksnumbers.h"
12 #include "kstarsdatetime.h"
13 #include "kssun.h"
14 #include "kstarsdata.h"
15 #include "Options.h"
16 #include "skyobject.h"
17 #include "skycomponents/skymapcomposite.h"
18
19 #include "config-kstars.h"
20
21 #include <KLocalizedString>
22
23 #include <QDebug>
24
25 #include <cmath>
26
27
28
29 // N.B. To use libnova depending on the availability of the package,
30 // uncomment the following:
31 /*
32 #ifdef HAVE_LIBNOVA
33 #define SKYPOINT_USE_LIBNOVA 1
34 #endif
35 */
36
37 #ifdef SKYPOINT_USE_LIBNOVA
38 #include <libnova/libnova.h>
39 bool SkyPoint::implementationIsLibnova = true;
40 #else
41 bool SkyPoint::implementationIsLibnova = false;
42 #endif
43
44 #ifdef PROFILE_COORDINATE_CONVERSION
45 #include <ctime> // For profiling, remove if not profiling.
46 long unsigned SkyPoint::eqToHzCalls = 0;
47 double SkyPoint::cpuTime_EqToHz = 0.;
48 #endif
49
50 KSSun *SkyPoint::m_Sun = nullptr;
51 const double SkyPoint::altCrit = -1.0;
52
SkyPoint()53 SkyPoint::SkyPoint()
54 {
55 // Default constructor. Set nonsense values
56 RA0.setD(-1); // RA >= 0 always :-)
57 Dec0.setD(180); // Dec is between -90 and 90 Degrees :-)
58 RA = RA0;
59 Dec = Dec0;
60 lastPrecessJD = J2000; // By convention, we use J2000 coordinates
61 }
62
set(const dms & r,const dms & d)63 void SkyPoint::set(const dms &r, const dms &d)
64 {
65 RA0 = RA = r;
66 Dec0 = Dec = d;
67 lastPrecessJD = J2000; // By convention, we use J2000 coordinates
68 }
69
EquatorialToHorizontal(const dms * LST,const dms * lat)70 void SkyPoint::EquatorialToHorizontal(const dms *LST, const dms *lat)
71 {
72 // qDebug() << "NOTE: This EquatorialToHorizontal overload (using dms pointers instead of CachingDms pointers) is deprecated and should be replaced with CachingDms prototype wherever speed is desirable!";
73 CachingDms _LST(*LST), _lat(*lat);
74 EquatorialToHorizontal(&_LST, &_lat);
75 }
76
EquatorialToHorizontal(const CachingDms * LST,const CachingDms * lat)77 void SkyPoint::EquatorialToHorizontal(const CachingDms *LST, const CachingDms *lat)
78 {
79 #ifdef PROFILE_COORDINATE_CONVERSION
80 std::clock_t start = std::clock();
81 #endif
82 //Uncomment for spherical trig version
83 double AltRad, AzRad;
84 double sindec, cosdec, sinlat, coslat, sinHA, cosHA;
85 double sinAlt, cosAlt;
86
87 CachingDms HourAngle =
88 (*LST) - ra(); // Using CachingDms subtraction operator to find cos/sin of HourAngle without calling sincos()
89
90 lat->SinCos(sinlat, coslat);
91 dec().SinCos(sindec, cosdec);
92 HourAngle.SinCos(sinHA, cosHA);
93
94 sinAlt = sindec * sinlat + cosdec * coslat * cosHA;
95 AltRad = asin(sinAlt);
96
97 cosAlt = sqrt(
98 1 -
99 sinAlt *
100 sinAlt); // Avoid trigonometric function. Return value of asin is always in [-pi/2, pi/2] and in this domain cosine is always non-negative, so we can use this.
101 if (cosAlt == 0.)
102 cosAlt = cos(AltRad);
103
104 double arg = (sindec - sinlat * sinAlt) / (coslat * cosAlt);
105 if (arg <= -1.0)
106 AzRad = dms::PI;
107 else if (arg >= 1.0)
108 AzRad = 0.0;
109 else
110 AzRad = acos(arg);
111
112 if (sinHA > 0.0 && AzRad != 0.0)
113 AzRad = 2.0 * dms::PI - AzRad; // resolve acos() ambiguity
114
115 Alt.setRadians(AltRad);
116 Az.setRadians(AzRad);
117 #ifdef PROFILE_COORDINATE_CONVERSION
118 std::clock_t stop = std::clock();
119 cpuTime_EqToHz += double(stop - start) / double(CLOCKS_PER_SEC); // Accumulate time in seconds
120 ++eqToHzCalls;
121 #endif
122
123 // //Uncomment for XYZ version
124 // double xr, yr, zr, xr1, zr1, sa, ca;
125 // //Z-axis rotation by -LST
126 // dms a = dms( -1.0*LST->Degrees() );
127 // a.SinCos( sa, ca );
128 // xr1 = m_X*ca + m_Y*sa;
129 // yr = -1.0*m_X*sa + m_Y*ca;
130 // zr1 = m_Z;
131 //
132 // //Y-axis rotation by lat - 90.
133 // a = dms( lat->Degrees() - 90.0 );
134 // a.SinCos( sa, ca );
135 // xr = xr1*ca - zr1*sa;
136 // zr = xr1*sa + zr1*ca;
137 //
138 // //FIXME: eventually, we will work with XYZ directly
139 // Alt.setRadians( asin( zr ) );
140 // Az.setRadians( atan2( yr, xr ) );
141 }
142
HorizontalToEquatorial(const dms * LST,const dms * lat)143 void SkyPoint::HorizontalToEquatorial(const dms *LST, const dms *lat)
144 {
145 double HARad, DecRad;
146 double sinlat, coslat, sinAlt, cosAlt, sinAz, cosAz;
147 double sinDec, cosDec;
148
149 lat->SinCos(sinlat, coslat);
150 alt().SinCos(sinAlt, cosAlt);
151 Az.SinCos(sinAz, cosAz);
152
153 sinDec = sinAlt * sinlat + cosAlt * coslat * cosAz;
154 DecRad = asin(sinDec);
155 cosDec = cos(DecRad);
156 Dec.setRadians(DecRad);
157
158 double x = (sinAlt - sinlat * sinDec) / (coslat * cosDec);
159
160 //Under certain circumstances, x can be very slightly less than -1.0000, or slightly
161 //greater than 1.0000, leading to a crash on acos(x). However, the value isn't
162 //*really* out of range; it's a kind of roundoff error.
163 if (x < -1.0 && x > -1.000001)
164 HARad = dms::PI;
165 else if (x > 1.0 && x < 1.000001)
166 HARad = 0.0;
167 else if (x < -1.0)
168 {
169 //qWarning() << "Coordinate out of range.";
170 HARad = dms::PI;
171 }
172 else if (x > 1.0)
173 {
174 //qWarning() << "Coordinate out of range.";
175 HARad = 0.0;
176 }
177 else
178 HARad = acos(x);
179
180 if (sinAz > 0.0)
181 HARad = 2.0 * dms::PI - HARad; // resolve acos() ambiguity
182
183 RA.setRadians(LST->radians() - HARad);
184 RA.reduceToRange(dms::ZERO_TO_2PI);
185 }
186
findEcliptic(const CachingDms * Obliquity,dms & EcLong,dms & EcLat)187 void SkyPoint::findEcliptic(const CachingDms *Obliquity, dms &EcLong, dms &EcLat)
188 {
189 double sinRA, cosRA, sinOb, cosOb, sinDec, cosDec;
190 ra().SinCos(sinRA, cosRA);
191 dec().SinCos(sinDec, cosDec);
192 Obliquity->SinCos(sinOb, cosOb);
193
194 double ycosDec = sinRA * cosOb * cosDec + sinDec * sinOb;
195 double ELongRad = atan2(ycosDec, cosDec * cosRA);
196 EcLong.setRadians(ELongRad);
197 EcLong.reduceToRange(dms::ZERO_TO_2PI);
198 EcLat.setRadians(asin(sinDec * cosOb - cosDec * sinOb * sinRA)); // FIXME: Haversine
199 }
200
setFromEcliptic(const CachingDms * Obliquity,const dms & EcLong,const dms & EcLat)201 void SkyPoint::setFromEcliptic(const CachingDms *Obliquity, const dms &EcLong, const dms &EcLat)
202 {
203 double sinLong, cosLong, sinLat, cosLat, sinObliq, cosObliq;
204 EcLong.SinCos(sinLong, cosLong);
205 EcLat.SinCos(sinLat, cosLat);
206 Obliquity->SinCos(sinObliq, cosObliq);
207
208 // double sinDec = sinLat * cosObliq + cosLat * sinObliq * sinLong;
209
210 double ycosLat = sinLong * cosObliq * cosLat - sinLat * sinObliq;
211 // double RARad = atan2( y, cosLong );
212 RA.setUsing_atan2(ycosLat, cosLong * cosLat);
213 RA.reduceToRange(dms::ZERO_TO_2PI);
214 // Dec.setUsing_asin(sinDec);
215
216 // Use Haversine to set declination
217 Dec.setRadians(dms::PI/2.0 - 2.0 * asin(sqrt(0.5 * (
218 1.0 - sinLat * cosObliq
219 - cosLat * sinObliq * sinLong
220 ))));
221 }
222
precess(const KSNumbers * num)223 void SkyPoint::precess(const KSNumbers *num)
224 {
225 double cosRA0, sinRA0, cosDec0, sinDec0;
226 const Eigen::Matrix3d &precessionMatrix = num->p2();
227 Eigen::Vector3d v, s;
228
229 RA0.SinCos(sinRA0, cosRA0);
230 Dec0.SinCos(sinDec0, cosDec0);
231
232 s[0] = cosRA0 * cosDec0;
233 s[1] = sinRA0 * cosDec0;
234 s[2] = sinDec0;
235
236 // NOTE: Rotation matrices are the fastest way to do rotations on
237 // a vector. Quaternions need more multiplications. The rotation
238 // matrix compensates in some sense by having more 'precomputed'
239 // multiplications. The matrix elements seem to cache nicely, so
240 // there isn't much overhead in accessing them.
241
242 //Multiply P2 and s to get v, the vector representing the new coords.
243 // for ( unsigned int i=0; i<3; ++i ) {
244 // v[i] = 0.0;
245 // for (uint j=0; j< 3; ++j) {
246 // v[i] += num->p2( j, i )*s[j];
247 // }
248 // }
249 v.noalias() = precessionMatrix * s;
250
251 //Extract RA, Dec from the vector:
252 RA.setUsing_atan2(v[1], v[0]);
253 RA.reduceToRange(dms::ZERO_TO_2PI);
254 Dec.setUsing_asin(v[2]);
255 }
256
deprecess(const KSNumbers * num,long double epoch)257 SkyPoint SkyPoint::deprecess(const KSNumbers *num, long double epoch)
258 {
259 SkyPoint p1(RA, Dec);
260 long double now = num->julianDay();
261 p1.precessFromAnyEpoch(now, epoch);
262 if ((std::isnan(RA0.Degrees()) || std::isnan(Dec0.Degrees())) ||
263 (!std::isnan(Dec0.Degrees()) && fabs(Dec0.Degrees()) > 90.0))
264 {
265 // We have invalid RA0 and Dec0, so set them if epoch = J2000. Otherwise, do not touch.
266 if (epoch == J2000L)
267 {
268 RA0 = p1.ra();
269 Dec0 = p1.dec();
270 }
271 }
272 return p1;
273 }
274
nutate(const KSNumbers * num,const bool reverse)275 void SkyPoint::nutate(const KSNumbers *num, const bool reverse)
276 {
277 //Step 2: Nutation
278 if (fabs(Dec.Degrees()) < 80.0) //approximate method
279 {
280 double cosRA, sinRA, cosDec, sinDec, tanDec;
281 double cosOb, sinOb;
282 double dRA, dDec;
283
284 RA.SinCos(sinRA, cosRA);
285 Dec.SinCos(sinDec, cosDec);
286
287 tanDec = sinDec / cosDec;
288
289 // Equ 23.1 in Jean Meeus' book
290 // nut_ecliptic / num->obliquity() is called epsilon in Meeus
291 // nut.longitude / num->dEcLong() is called delta psi in Meeus
292 // nut.obliquity / num->dObliq() is called delta epsilon in Meeus
293 // Meeus notes that these expressions are invalid if the star is
294 // close to one of the celestial poles
295
296 // N.B. These expressions are valid for FK5 coordinates
297 // (presumably also valid for ICRS by extension), but not for
298 // FK4. See the "Important Remark" on Page 152 of Jean Meeus'
299 // book.
300
301 #ifdef SKYPOINT_USE_LIBNOVA
302 // code lifted from libnova ln_get_equ_nut, tailored to our needs
303 // with the option to add or remove nutation
304 struct ln_nutation nut;
305 ln_get_nutation (num->julianDay(), &nut); // FIXME: Is this cached, or is it a slow call? If it is slow, we should move it to KSNumbers
306
307 double nut_ecliptic = ln_deg_to_rad(nut.ecliptic + nut.obliquity);
308 sinOb = sin(nut_ecliptic);
309 cosOb = cos(nut_ecliptic);
310
311 dRA = nut.longitude * (cosOb + sinOb * sinRA * tanDec) - nut.obliquity * cosRA * tanDec;
312 dDec = nut.longitude * (sinOb * cosRA) + nut.obliquity * sinRA;
313 #else
314 num->obliquity()->SinCos(sinOb, cosOb);
315
316 // N.B. num->dEcLong() and num->dObliq() methods return in
317 // degrees, whereby the corrections dRA and dDec are also in
318 // degrees.
319 dRA = num->dEcLong() * (cosOb + sinOb * sinRA * tanDec) - num->dObliq() * cosRA * tanDec;
320 dDec = num->dEcLong() * (sinOb * cosRA) + num->dObliq() * sinRA;
321 #endif
322 // the sign changed to remove nutation
323 if (reverse)
324 {
325 dRA = -dRA;
326 dDec = -dDec;
327 }
328 RA.setD(RA.Degrees() + dRA);
329 Dec.setD(Dec.Degrees() + dDec);
330 }
331 else //exact method
332 {
333 // NOTE: Meeus declares that you must add Δψ to the ecliptic
334 // longitude of the body to get a more accurate precession
335 // result, but fails to emphasize that the NCP of the two
336 // coordinates systems is different. The (RA, Dec) without
337 // nutation computed, i.e. the mean place of the sky point is
338 // referenced to the un-nutated geocentric frame (without the
339 // obliquity correction), whereas the (RA, Dec) after nutation
340 // is applied is referenced to the nutated geocentric
341 // frame. This is more clearly explained in the "Explanatory
342 // Supplement to the Astronomical Almanac" by
343 // K. P. Seidelmann, which can be borrowed on the internet
344 // archive as of this writing, see page 114:
345 // https://archive.org/details/explanatorysuppl00pken
346 //
347 // The rotation matrix formulation in (3.222-3) and the figure
348 // 3.222.1 make this clear
349
350 // TODO apply reverse test above 80 degrees
351 dms EcLong, EcLat;
352 CachingDms obliquityWithoutNutation(*num->obliquity() - dms(num->dObliq()));
353
354 if (reverse)
355 {
356 //Subtract dEcLong from the Ecliptic Longitude
357 findEcliptic(num->obliquity(), EcLong, EcLat);
358 dms newLong(EcLong.Degrees() - num->dEcLong());
359 setFromEcliptic(&obliquityWithoutNutation, newLong, EcLat); // FIXME: Check
360 }
361 else
362 {
363 //Add dEcLong to the Ecliptic Longitude
364 findEcliptic(&obliquityWithoutNutation, EcLong, EcLat);
365 dms newLong(EcLong.Degrees() + num->dEcLong());
366 setFromEcliptic(num->obliquity(), newLong, EcLat);
367 }
368 }
369 }
370
moveAway(const SkyPoint & from,double dist) const371 SkyPoint SkyPoint::moveAway(const SkyPoint &from, double dist) const
372 {
373 CachingDms lat1, dtheta;
374
375 if (dist == 0.0)
376 {
377 qDebug() << "moveAway called with zero distance!";
378 return *this;
379 }
380
381 double dst = fabs(dist * dms::DegToRad / 3600.0); // In radian
382
383 // Compute the bearing angle w.r.t. the RA axis ("latitude")
384 CachingDms dRA(ra() - from.ra());
385 CachingDms dDec(dec() - from.dec());
386 double bearing = atan2(dRA.sin() / dRA.cos(), dDec.sin()); // Do not use dRA = PI / 2!!
387 //double bearing = atan2( dDec.radians() , dRA.radians() );
388
389 // double dir0 = (dist >= 0 ) ? bearing : bearing + dms::PI; // in radian
390 double dir0 = bearing + std::signbit(dist) * dms::PI; // might be faster?
391 double sinDst = sin(dst), cosDst = cos(dst);
392
393 lat1.setUsing_asin(dec().sin() * cosDst + dec().cos() * sinDst * cos(dir0));
394 dtheta.setUsing_atan2(sin(dir0) * sinDst * dec().cos(), cosDst - dec().sin() * lat1.sin());
395
396 return SkyPoint(ra() + dtheta, lat1);
397 }
398
checkBendLight()399 bool SkyPoint::checkBendLight()
400 {
401 // First see if we are close enough to the sun to bother about the
402 // gravitational lensing effect. We correct for the effect at
403 // least till b = 10 solar radii, where the effect is only about
404 // 0.06". Assuming min. sun-earth distance is 200 solar radii.
405 static const dms maxAngle(1.75 * (30.0 / 200.0) / dms::DegToRad);
406
407 if (!m_Sun)
408 {
409 SkyComposite *skycomopsite = KStarsData::Instance()->skyComposite();
410
411 if (skycomopsite == nullptr)
412 return false;
413
414 m_Sun = dynamic_cast<KSSun *>(skycomopsite->findByName(i18n("Sun")));
415
416 if (m_Sun == nullptr)
417 return false;
418 }
419
420 // TODO: This can be optimized further. We only need a ballpark estimate of the distance to the sun to start with.
421 return (fabs(angularDistanceTo(static_cast<const SkyPoint *>(m_Sun)).Degrees()) <=
422 maxAngle.Degrees()); // NOTE: dynamic_cast is slow and not important here.
423 }
424
bendlight()425 bool SkyPoint::bendlight()
426 {
427 // NOTE: This should be applied before aberration
428 // NOTE: One must call checkBendLight() before unnecessarily calling this.
429 // We correct for GR effects
430
431 // NOTE: This code is buggy. The sun needs to be initialized to
432 // the current epoch -- but we are not certain that this is the
433 // case. We have, as of now, no way of telling if the sun is
434 // initialized or not. If we initialize the sun here, we will be
435 // slowing down the program rather substantially and potentially
436 // introducing bugs. Therefore, we just ignore this problem, and
437 // hope that whenever the user is interested in seeing the effects
438 // of GR, we have the sun initialized correctly. This is usually
439 // the case. When the sun is not correctly initialized, rearth()
440 // is not computed, so we just assume it is nominally equal to 1
441 // AU to get a reasonable estimate.
442 Q_ASSERT(m_Sun);
443 double corr_sec = 1.75 * m_Sun->physicalSize() /
444 ((std::isfinite(m_Sun->rearth()) ? m_Sun->rearth() : 1) * AU_KM *
445 angularDistanceTo(static_cast<const SkyPoint *>(m_Sun)).sin());
446 Q_ASSERT(corr_sec > 0);
447
448 SkyPoint sp = moveAway(*m_Sun, corr_sec);
449 setRA(sp.ra());
450 setDec(sp.dec());
451 return true;
452 }
453
aberrate(const KSNumbers * num,bool reverse)454 void SkyPoint::aberrate(const KSNumbers *num, bool reverse)
455 {
456 #ifdef SKYPOINT_USE_LIBNOVA
457 ln_equ_posn pos { RA.Degrees(), Dec.Degrees() };
458 ln_equ_posn abPos { 0, 0 };
459 ln_get_equ_aber(&pos, num->julianDay(), &abPos);
460 if (reverse)
461 {
462 RA.setD(RA.Degrees() * 2 - abPos.ra);
463 Dec.setD(Dec.Degrees() * 2 - abPos.dec);
464 }
465 else
466 {
467 RA.setD(abPos.ra);
468 Dec.setD(abPos.dec);
469 }
470
471 #else
472 // N.B. These expressions are valid for FK5 coordinates
473 // (presumably also valid for ICRS by extension), but not for
474 // FK4. See the "Important Remark" on Page 152 of Jean Meeus'
475 // book.
476
477 // N.B. Even though Meeus does not say this explicitly, these
478 // equations must not apply near the pole. Therefore, we fall-back
479 // to the expressions provided by Meeus in ecliptic coordinates
480 // (Equ 23.2) when we are near the pole.
481
482 double K = num->constAberr().Degrees(); //constant of aberration
483 double e = num->earthEccentricity(); // eccentricity of Earth's orbit
484
485 if (fabs(Dec.Degrees()) < 80.0)
486 {
487
488 double cosRA, sinRA, cosDec, sinDec;
489 double cosL, sinL, cosP, sinP;
490 double cosOb, sinOb;
491
492
493 RA.SinCos(sinRA, cosRA);
494 Dec.SinCos(sinDec, cosDec);
495
496 num->obliquity()->SinCos(sinOb, cosOb);
497 // double tanOb = sinOb/cosOb;
498
499 num->sunTrueLongitude().SinCos(sinL, cosL);
500 num->earthPerihelionLongitude().SinCos(sinP, cosP);
501
502 //Step 3: Aberration
503 // These are the expressions given in Jean Meeus, Equation (23.3)
504
505 // double dRA = -1.0 * K * ( cosRA * cosL * cosOb + sinRA * sinL )/cosDec
506 // + e * K * ( cosRA * cosP * cosOb + sinRA * sinP )/cosDec;
507
508 // double dDec = -1.0 * K * ( cosL * cosOb * ( tanOb * cosDec - sinRA * sinDec ) + cosRA * sinDec * sinL )
509 // + e * K * ( cosP * cosOb * ( tanOb * cosDec - sinRA * sinDec ) + cosRA * sinDec * sinP );
510
511 // However, we have factorized the expressions below by pulling
512 // out common factors to make it more optimized, in case the
513 // compiler fails to spot these optimizations.
514
515 // N.B. I had failed to factor out the expressions correctly,
516 // making mistakes, in c5e709bd91, which should now be
517 // fixed. --asimha
518
519 // FIXME: Check if the unit tests have sufficient coverage to
520 // check this expression
521
522 double dRA = (K / cosDec) * (
523 cosRA * cosOb * (e * cosP - cosL)
524 + sinRA * (e * sinP - sinL)
525 );
526 double dDec = K * (
527 (sinOb * cosDec - cosOb * sinDec * sinRA) * (e * cosP - cosL)
528 + cosRA * sinDec * (e * sinP - sinL)
529 );
530
531 // N.B. Meeus points out that the result has the same units as
532 // K, so the corrections are in degrees.
533
534 if (reverse)
535 {
536 dRA = -dRA;
537 dDec = -dDec;
538 }
539 RA.setD(RA.Degrees() + dRA);
540 Dec.setD(Dec.Degrees() + dDec);
541 }
542 else
543 {
544 dms EcLong, EcLat;
545 double sinEcLat, cosEcLat;
546 const auto &L = num->sunTrueLongitude();
547 const auto &P = num->earthPerihelionLongitude();
548
549 findEcliptic(num->obliquity(), EcLong, EcLat);
550 EcLat.SinCos(sinEcLat, cosEcLat);
551
552 double sin_L_minus_EcLong, cos_L_minus_EcLong;
553 double sin_P_minus_EcLong, cos_P_minus_EcLong;
554 (L - EcLong).SinCos(sin_L_minus_EcLong, cos_L_minus_EcLong);
555 (P - EcLong).SinCos(sin_P_minus_EcLong, cos_P_minus_EcLong);
556
557 // Equation (23.2) in Meeus
558 // N.B. dEcLong, dEcLat are in degrees, because K is expressed in degrees.
559 double dEcLong = (K / cosEcLat) * (e * cos_P_minus_EcLong - cos_L_minus_EcLong);
560 double dEcLat = K * sinEcLat * (e * sin_P_minus_EcLong - sin_L_minus_EcLong);
561
562 // Note: These are approximate corrections, so it is
563 // appropriate to change their sign to reverse them to first
564 // order in the corrections. As a result, the forward and
565 // reverse calculations will not be exact inverses, but only
566 // approximate inverses.
567 if (reverse)
568 {
569 dEcLong = -dEcLong;
570 dEcLat = -dEcLat;
571 }
572
573 // Update the ecliptic coordinates to their new values
574 EcLong.setD(EcLong.Degrees() + dEcLong);
575 EcLat.setD(EcLat.Degrees() + dEcLat);
576 setFromEcliptic(num->obliquity(), EcLong, EcLat);
577 }
578 #endif
579 }
580
581 // Note: This method is one of the major rate determining factors in how fast the map pans / zooms in or out
updateCoords(const KSNumbers * num,bool,const CachingDms * lat,const CachingDms * LST,bool forceRecompute)582 void SkyPoint::updateCoords(const KSNumbers *num, bool /*includePlanets*/, const CachingDms *lat, const CachingDms *LST,
583 bool forceRecompute)
584 {
585 //Correct the catalog coordinates for the time-dependent effects
586 //of precession, nutation and aberration
587 bool recompute, lens;
588
589 // NOTE: The same short-circuiting checks are also implemented in
590 // StarObject::JITUpdate(), even before calling
591 // updateCoords(). While this is code-duplication, these bits of
592 // code need to be really optimized, at least for stars. For
593 // optimization purposes, the code is left duplicated in two
594 // places. Please be wary of changing one without changing the
595 // other.
596
597 Q_ASSERT(std::isfinite(lastPrecessJD));
598
599 if (Options::useRelativistic() && checkBendLight())
600 {
601 recompute = true;
602 lens = true;
603 }
604 else
605 {
606 recompute = (Options::alwaysRecomputeCoordinates() || forceRecompute ||
607 std::abs(lastPrecessJD - num->getJD()) >= 0.00069444); // Update once per solar minute
608 lens = false;
609 }
610 if (recompute)
611 {
612 precess(num);
613 nutate(num);
614 if (lens)
615 bendlight(); // FIXME: Shouldn't we apply this on the horizontal coordinates?
616 aberrate(num);
617 lastPrecessJD = num->getJD();
618 Q_ASSERT(std::isfinite(RA.Degrees()) && std::isfinite(Dec.Degrees()));
619 }
620
621 if (lat || LST)
622 qWarning() << i18n("lat and LST parameters should only be used in KSPlanetBase objects.");
623 }
624
precessFromAnyEpoch(long double jd0,long double jdf)625 void SkyPoint::precessFromAnyEpoch(long double jd0, long double jdf)
626 {
627 double cosRA, sinRA, cosDec, sinDec;
628 double v[3], s[3];
629
630 RA = RA0;
631 Dec = Dec0; // Is this necessary?
632
633 if (jd0 == jdf)
634 return;
635
636 RA.SinCos(sinRA, cosRA);
637 Dec.SinCos(sinDec, cosDec);
638
639 if (jd0 == B1950L)
640 {
641 B1950ToJ2000();
642 jd0 = J2000L;
643 RA.SinCos(sinRA, cosRA);
644 Dec.SinCos(sinDec, cosDec);
645 }
646
647 if (jd0 != jdf)
648 {
649 // The original coordinate is referred to the FK5 system and
650 // is NOT J2000.
651 if (jd0 != J2000L)
652 {
653 //v is a column vector representing input coordinates.
654 v[0] = cosRA * cosDec;
655 v[1] = sinRA * cosDec;
656 v[2] = sinDec;
657
658 //Need to first precess to J2000.0 coords
659 //s is the product of P1 and v; s represents the
660 //coordinates precessed to J2000
661 KSNumbers num(jd0);
662 for (unsigned int i = 0; i < 3; ++i)
663 {
664 s[i] = num.p1(0, i) * v[0] + num.p1(1, i) * v[1] + num.p1(2, i) * v[2];
665 }
666
667 //Input coords already in J2000, set s accordingly.
668 }
669 else
670 {
671 s[0] = cosRA * cosDec;
672 s[1] = sinRA * cosDec;
673 s[2] = sinDec;
674 }
675
676 if (jdf == B1950L)
677 {
678 RA.setRadians(atan2(s[1], s[0]));
679 Dec.setRadians(asin(s[2]));
680 J2000ToB1950();
681
682 return;
683 }
684
685 KSNumbers num(jdf);
686 for (unsigned int i = 0; i < 3; ++i)
687 {
688 v[i] = num.p2(0, i) * s[0] + num.p2(1, i) * s[1] + num.p2(2, i) * s[2];
689 }
690
691 RA.setUsing_atan2(v[1], v[0]);
692 Dec.setUsing_asin(v[2]);
693
694 RA.reduceToRange(dms::ZERO_TO_2PI);
695
696 return;
697 }
698 }
699
apparentCoord(long double jd0,long double jdf)700 void SkyPoint::apparentCoord(long double jd0, long double jdf)
701 {
702 precessFromAnyEpoch(jd0, jdf);
703 KSNumbers num(jdf);
704 nutate(&num);
705 if (Options::useRelativistic() && checkBendLight())
706 bendlight();
707 aberrate(&num);
708 }
709
catalogueCoord(long double jdf)710 SkyPoint SkyPoint::catalogueCoord(long double jdf)
711 {
712 KSNumbers num(jdf);
713
714 // remove abberation
715 aberrate(&num, true);
716
717 // remove nutation
718 nutate(&num, true);
719
720 // remove precession
721 // the start position needs to be in RA0,Dec0
722 RA0 = RA;
723 Dec0 = Dec;
724 // from now to J2000
725 precessFromAnyEpoch(jdf, static_cast<long double>(J2000));
726 // the J2000 position is in RA,Dec, move to RA0, Dec0
727 RA0 = RA;
728 Dec0 = Dec;
729 lastPrecessJD = J2000;
730
731 SkyPoint sp(RA0, Dec0);
732 return sp;
733 }
734
Equatorial1950ToGalactic(dms & galLong,dms & galLat)735 void SkyPoint::Equatorial1950ToGalactic(dms &galLong, dms &galLat)
736 {
737 double a = 192.25;
738 double sinb, cosb, sina_RA, cosa_RA, sinDEC, cosDEC, tanDEC;
739
740 dms c(303.0);
741 dms b(27.4);
742 tanDEC = tan(Dec.radians());
743
744 b.SinCos(sinb, cosb);
745 dms(a - RA.Degrees()).SinCos(sina_RA, cosa_RA);
746 Dec.SinCos(sinDEC, cosDEC);
747
748 galLong.setRadians(c.radians() - atan2(sina_RA, cosa_RA * sinb - tanDEC * cosb));
749 galLong.reduceToRange(dms::ZERO_TO_2PI);
750
751 galLat.setRadians(asin(sinDEC * sinb + cosDEC * cosb * cosa_RA));
752 }
753
GalacticToEquatorial1950(const dms * galLong,const dms * galLat)754 void SkyPoint::GalacticToEquatorial1950(const dms *galLong, const dms *galLat)
755 {
756 double a = 123.0;
757 double sinb, cosb, singLat, cosgLat, tangLat, singLong_a, cosgLong_a;
758
759 dms c(12.25);
760 dms b(27.4);
761 tangLat = tan(galLat->radians());
762 galLat->SinCos(singLat, cosgLat);
763
764 dms(galLong->Degrees() - a).SinCos(singLong_a, cosgLong_a);
765 b.SinCos(sinb, cosb);
766
767 RA.setRadians(c.radians() + atan2(singLong_a, cosgLong_a * sinb - tangLat * cosb));
768 RA.reduceToRange(dms::ZERO_TO_2PI);
769
770 Dec.setRadians(asin(singLat * sinb + cosgLat * cosb * cosgLong_a));
771 }
772
B1950ToJ2000(void)773 void SkyPoint::B1950ToJ2000(void)
774 {
775 double cosRA, sinRA, cosDec, sinDec;
776 // double cosRA0, sinRA0, cosDec0, sinDec0;
777 double v[3], s[3];
778
779 // 1984 January 1 0h
780 KSNumbers num(2445700.5L);
781
782 // Eterms due to aberration
783 addEterms();
784 RA.SinCos(sinRA, cosRA);
785 Dec.SinCos(sinDec, cosDec);
786
787 // Precession from B1950 to J1984
788 s[0] = cosRA * cosDec;
789 s[1] = sinRA * cosDec;
790 s[2] = sinDec;
791 for (unsigned int i = 0; i < 3; ++i)
792 {
793 v[i] = num.p2b(0, i) * s[0] + num.p2b(1, i) * s[1] + num.p2b(2, i) * s[2];
794 }
795
796 // RA zero-point correction at 1984 day 1, 0h.
797 RA.setRadians(atan2(v[1], v[0]));
798 Dec.setRadians(asin(v[2]));
799
800 RA.setH(RA.Hours() + 0.06390 / 3600.);
801 RA.SinCos(sinRA, cosRA);
802 Dec.SinCos(sinDec, cosDec);
803
804 s[0] = cosRA * cosDec;
805 s[1] = sinRA * cosDec;
806 s[2] = sinDec;
807
808 // Precession from 1984 to J2000.
809
810 for (unsigned int i = 0; i < 3; ++i)
811 {
812 v[i] = num.p1(0, i) * s[0] + num.p1(1, i) * s[1] + num.p1(2, i) * s[2];
813 }
814
815 RA.setRadians(atan2(v[1], v[0]));
816 Dec.setRadians(asin(v[2]));
817 }
818
J2000ToB1950(void)819 void SkyPoint::J2000ToB1950(void)
820 {
821 double cosRA, sinRA, cosDec, sinDec;
822 // double cosRA0, sinRA0, cosDec0, sinDec0;
823 double v[3], s[3];
824
825 // 1984 January 1 0h
826 KSNumbers num(2445700.5L);
827
828 RA.SinCos(sinRA, cosRA);
829 Dec.SinCos(sinDec, cosDec);
830
831 s[0] = cosRA * cosDec;
832 s[1] = sinRA * cosDec;
833 s[2] = sinDec;
834
835 // Precession from J2000 to 1984 day, 0h.
836
837 for (unsigned int i = 0; i < 3; ++i)
838 {
839 v[i] = num.p2(0, i) * s[0] + num.p2(1, i) * s[1] + num.p2(2, i) * s[2];
840 }
841
842 RA.setRadians(atan2(v[1], v[0]));
843 Dec.setRadians(asin(v[2]));
844
845 // RA zero-point correction at 1984 day 1, 0h.
846
847 RA.setH(RA.Hours() - 0.06390 / 3600.);
848 RA.SinCos(sinRA, cosRA);
849 Dec.SinCos(sinDec, cosDec);
850
851 // Precession from B1950 to J1984
852
853 s[0] = cosRA * cosDec;
854 s[1] = sinRA * cosDec;
855 s[2] = sinDec;
856 for (unsigned int i = 0; i < 3; ++i)
857 {
858 v[i] = num.p1b(0, i) * s[0] + num.p1b(1, i) * s[1] + num.p1b(2, i) * s[2];
859 }
860
861 RA.setRadians(atan2(v[1], v[0]));
862 Dec.setRadians(asin(v[2]));
863
864 // Eterms due to aberration
865 subtractEterms();
866 }
867
Eterms(void)868 SkyPoint SkyPoint::Eterms(void)
869 {
870 double sd, cd, sinEterm, cosEterm;
871 dms raTemp, raDelta, decDelta;
872
873 Dec.SinCos(sd, cd);
874 raTemp.setH(RA.Hours() + 11.25);
875 raTemp.SinCos(sinEterm, cosEterm);
876
877 raDelta.setH(0.0227 * sinEterm / (3600. * cd));
878 decDelta.setD(0.341 * cosEterm * sd / 3600. + 0.029 * cd / 3600.);
879
880 return SkyPoint(raDelta, decDelta);
881 }
882
addEterms(void)883 void SkyPoint::addEterms(void)
884 {
885 SkyPoint spd = Eterms();
886
887 RA = RA + spd.ra();
888 Dec = Dec + spd.dec();
889 }
890
subtractEterms(void)891 void SkyPoint::subtractEterms(void)
892 {
893 SkyPoint spd = Eterms();
894
895 RA = RA - spd.ra();
896 Dec = Dec - spd.dec();
897 }
898
angularDistanceTo(const SkyPoint * sp,double * const positionAngle) const899 dms SkyPoint::angularDistanceTo(const SkyPoint *sp, double *const positionAngle) const
900 {
901 // double dalpha = sp->ra().radians() - ra().radians() ;
902 // double ddelta = sp->dec().radians() - dec().radians();
903 CachingDms dalpha = sp->ra() - ra();
904 CachingDms ddelta = sp->dec() - dec();
905
906 // double sa = sin(dalpha/2.);
907 // double sd = sin(ddelta/2.);
908
909 // double hava = sa*sa;
910 // double havd = sd*sd;
911
912 // Compute the haversin directly:
913 double hava = (1 - dalpha.cos()) / 2.;
914 double havd = (1 - ddelta.cos()) / 2.;
915
916 // Haversine law
917 double aux = havd + (sp->dec().cos()) * dec().cos() * hava;
918
919 dms angDist;
920 angDist.setRadians(2. * fabs(asin(sqrt(aux))));
921
922 if (positionAngle)
923 {
924 // Also compute the position angle of the line from this SkyPoint to sp
925 //*positionAngle = acos( tan(-ddelta)/tan( angDist.radians() ) ); // FIXME: Might fail for large ddelta / zero angDist
926 //if( -dalpha < 0 )
927 // *positionAngle = 2*M_PI - *positionAngle;
928 *positionAngle =
929 atan2f(dalpha.sin(), (dec().cos()) * tan(sp->dec().radians()) - (dec().sin()) * dalpha.cos()) * 180 / M_PI;
930 }
931 return angDist;
932 }
933
vRSun(long double jd0)934 double SkyPoint::vRSun(long double jd0)
935 {
936 double ca, sa, cd, sd, vsun;
937 double cosRA, sinRA, cosDec, sinDec;
938
939 /* Sun apex (where the sun goes) coordinates */
940
941 dms asun(270.9592); // Right ascention: 18h 3m 50.2s [J2000]
942 dms dsun(30.00467); // Declination: 30^o 0' 16.8'' [J2000]
943 vsun = 20.; // [km/s]
944
945 asun.SinCos(sa, ca);
946 dsun.SinCos(sd, cd);
947
948 /* We need an auxiliary SkyPoint since we need the
949 * source referred to the J2000 equinox and we do not want to overwrite
950 * the current values
951 */
952
953 SkyPoint aux;
954 aux.set(RA0, Dec0);
955
956 aux.precessFromAnyEpoch(jd0, J2000L);
957
958 aux.ra().SinCos(sinRA, cosRA);
959 aux.dec().SinCos(sinDec, cosDec);
960
961 /* Computation is done performing the scalar product of a unitary vector
962 in the direction of the source with the vector velocity of Sun, both being in the
963 LSR reference system:
964 Vlsr = Vhel + Vsun.u_radial =>
965 Vlsr = Vhel + vsun(cos D cos A,cos D sen A,sen D).(cos d cos a,cos d sen a,sen d)
966 Vhel = Vlsr - Vsun.u_radial
967 */
968
969 return vsun * (cd * cosDec * (cosRA * ca + sa * sinRA) + sd * sinDec);
970 }
971
vHeliocentric(double vlsr,long double jd0)972 double SkyPoint::vHeliocentric(double vlsr, long double jd0)
973 {
974 return vlsr - vRSun(jd0);
975 }
976
vHelioToVlsr(double vhelio,long double jd0)977 double SkyPoint::vHelioToVlsr(double vhelio, long double jd0)
978 {
979 return vhelio + vRSun(jd0);
980 }
981
vREarth(long double jd0)982 double SkyPoint::vREarth(long double jd0)
983 {
984 double sinRA, sinDec, cosRA, cosDec;
985
986 /* u_radial = unitary vector in the direction of the source
987 Vlsr = Vhel + Vsun.u_radial
988 = Vgeo + VEarth.u_radial + Vsun.u_radial =>
989
990 Vgeo = (Vlsr -Vsun.u_radial) - VEarth.u_radial
991 = Vhel - VEarth.u_radial
992 = Vhel - (vx, vy, vz).(cos d cos a,cos d sen a,sen d)
993 */
994
995 /* We need an auxiliary SkyPoint since we need the
996 * source referred to the J2000 equinox and we do not want to overwrite
997 * the current values
998 */
999
1000 SkyPoint aux(RA0, Dec0);
1001
1002 aux.precessFromAnyEpoch(jd0, J2000L);
1003
1004 aux.ra().SinCos(sinRA, cosRA);
1005 aux.dec().SinCos(sinDec, cosDec);
1006
1007 /* vEarth is referred to the J2000 equinox, hence we need that
1008 the source coordinates are also in the same reference system.
1009 */
1010
1011 KSNumbers num(jd0);
1012 return num.vEarth(0) * cosDec * cosRA + num.vEarth(1) * cosDec * sinRA + num.vEarth(2) * sinDec;
1013 }
1014
vGeocentric(double vhelio,long double jd0)1015 double SkyPoint::vGeocentric(double vhelio, long double jd0)
1016 {
1017 return vhelio - vREarth(jd0);
1018 }
1019
vGeoToVHelio(double vgeo,long double jd0)1020 double SkyPoint::vGeoToVHelio(double vgeo, long double jd0)
1021 {
1022 return vgeo + vREarth(jd0);
1023 }
1024
vRSite(double vsite[3])1025 double SkyPoint::vRSite(double vsite[3])
1026 {
1027 double sinRA, sinDec, cosRA, cosDec;
1028
1029 RA.SinCos(sinRA, cosRA);
1030 Dec.SinCos(sinDec, cosDec);
1031
1032 return vsite[0] * cosDec * cosRA + vsite[1] * cosDec * sinRA + vsite[2] * sinDec;
1033 }
1034
vTopoToVGeo(double vtopo,double vsite[3])1035 double SkyPoint::vTopoToVGeo(double vtopo, double vsite[3])
1036 {
1037 return vtopo + vRSite(vsite);
1038 }
1039
vTopocentric(double vgeo,double vsite[3])1040 double SkyPoint::vTopocentric(double vgeo, double vsite[3])
1041 {
1042 return vgeo - vRSite(vsite);
1043 }
1044
checkCircumpolar(const dms * gLat) const1045 bool SkyPoint::checkCircumpolar(const dms *gLat) const
1046 {
1047 return fabs(dec().Degrees()) > (90 - fabs(gLat->Degrees()));
1048 }
1049
altRefracted() const1050 dms SkyPoint::altRefracted() const
1051 {
1052 return refract(Alt, Options::useRefraction());
1053 }
1054
setAltRefracted(dms alt_apparent)1055 void SkyPoint::setAltRefracted(dms alt_apparent)
1056 {
1057 setAlt(unrefract(alt_apparent, Options::useRefraction()));
1058 }
1059
setAltRefracted(double alt_apparent)1060 void SkyPoint::setAltRefracted(double alt_apparent)
1061 {
1062 setAlt(unrefract(alt_apparent, Options::useRefraction()));
1063 }
1064
refractionCorr(double alt)1065 double SkyPoint::refractionCorr(double alt)
1066 {
1067 return 1.02 / tan(dms::DegToRad * (alt + 10.3 / (alt + 5.11))) / 60;
1068 }
1069
refract(const double alt,bool conditional)1070 double SkyPoint::refract(const double alt, bool conditional)
1071 {
1072 if (!conditional)
1073 {
1074 return alt;
1075 }
1076 static double corrCrit = SkyPoint::refractionCorr(SkyPoint::altCrit);
1077
1078 if (alt > SkyPoint::altCrit)
1079 return (alt + SkyPoint::refractionCorr(alt));
1080 else
1081 return (alt +
1082 corrCrit * (alt + 90) /
1083 (SkyPoint::altCrit + 90)); // Linear extrapolation from corrCrit at altCrit to 0 at -90 degrees
1084 }
1085
1086 // Found uncorrected value by solving equation. This is OK since
1087 // unrefract is never called in loops with the potential exception of
1088 // slewing.
1089 //
1090 // Convergence is quite fast just a few iterations.
unrefract(const double alt,bool conditional)1091 double SkyPoint::unrefract(const double alt, bool conditional)
1092 {
1093 if (!conditional)
1094 {
1095 return alt;
1096 }
1097 double h0 = alt;
1098 double h1 =
1099 alt -
1100 (refract(h0) -
1101 h0); // It's probably okay to add h0 in refract() and subtract it here, since refract() is called way more frequently.
1102
1103 while (fabs(h1 - h0) > 1e-4)
1104 {
1105 h0 = h1;
1106 h1 = alt - (refract(h0) - h0);
1107 }
1108 return h1;
1109 }
1110
findAltitude(const SkyPoint * p,const KStarsDateTime & dt,const GeoLocation * geo,const double hour)1111 dms SkyPoint::findAltitude(const SkyPoint *p, const KStarsDateTime &dt, const GeoLocation *geo, const double hour)
1112 {
1113 Q_ASSERT(p);
1114 if (!p)
1115 return dms(NaN::d);
1116
1117 // Jasem 2015-08-24 Using correct procedure to find altitude
1118 return SkyPoint::timeTransformed(p, dt, geo, hour).alt();
1119 }
1120
timeTransformed(const SkyPoint * p,const KStarsDateTime & dt,const GeoLocation * geo,const double hour)1121 SkyPoint SkyPoint::timeTransformed(const SkyPoint *p, const KStarsDateTime &dt, const GeoLocation *geo,
1122 const double hour)
1123 {
1124 Q_ASSERT(p);
1125 if (!p)
1126 return SkyPoint(NaN::d, NaN::d);
1127
1128 // Jasem 2015-08-24 Using correct procedure to find altitude
1129 SkyPoint sp = *p; // make a copy
1130 KStarsDateTime targetDateTime = dt.addSecs(hour * 3600.0);
1131 dms LST = geo->GSTtoLST(targetDateTime.gst());
1132 sp.EquatorialToHorizontal(&LST, geo->lat());
1133 return sp;
1134 }
1135
maxAlt(const dms & lat) const1136 double SkyPoint::maxAlt(const dms &lat) const
1137 {
1138 double retval = (lat.Degrees() + 90. - dec().Degrees());
1139 if (retval > 90.)
1140 retval = 180. - retval;
1141 return retval;
1142 }
1143
minAlt(const dms & lat) const1144 double SkyPoint::minAlt(const dms &lat) const
1145 {
1146 double retval = (lat.Degrees() - 90. + dec().Degrees());
1147 if (retval < -90.)
1148 retval = 180. + retval;
1149 return retval;
1150 }
1151
1152 #ifndef KSTARS_LITE
operator <<(QDBusArgument & argument,const SkyPoint & source)1153 QDBusArgument &operator<<(QDBusArgument &argument, const SkyPoint &source)
1154 {
1155 argument.beginStructure();
1156 argument << source.ra().Hours() << source.dec().Degrees();
1157 argument.endStructure();
1158 return argument;
1159 }
1160
operator >>(const QDBusArgument & argument,SkyPoint & dest)1161 const QDBusArgument &operator>>(const QDBusArgument &argument, SkyPoint &dest)
1162 {
1163 double ra, dec;
1164 argument.beginStructure();
1165 argument >> ra >> dec;
1166 argument.endStructure();
1167 dest = SkyPoint(ra, dec);
1168 return argument;
1169 }
1170 #endif
1171
1172