Kstars

skypoint.cpp
1 /*
2  SPDX-FileCopyrightText: 2001-2005 Jason Harris <[email protected]>
3  SPDX-FileCopyrightText: 2004-2005 Pablo de Vicente <[email protected]>
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 
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 
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 
70 void SkyPoint::EquatorialToHorizontal(const dms *LST, const dms *lat)
71 {
72  // qDebug() << Q_FUNC_INFO << "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 
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 
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 
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 
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 
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 
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 
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 
371 SkyPoint SkyPoint::moveAway(const SkyPoint &from, double dist) const
372 {
373  CachingDms lat1, dtheta;
374 
375  if (dist == 0.0)
376  {
377  qDebug() << Q_FUNC_INFO << "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 
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 
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 
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
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 
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 
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 
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 
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 
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 
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 
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 
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 
884 {
885  SkyPoint spd = Eterms();
886 
887  RA = RA + spd.ra();
888  Dec = Dec + spd.dec();
889 }
890 
892 {
893  SkyPoint spd = Eterms();
894 
895  RA = RA - spd.ra();
896  Dec = Dec - spd.dec();
897 }
898 
899 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 
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 
972 double SkyPoint::vHeliocentric(double vlsr, long double jd0)
973 {
974  return vlsr - vRSun(jd0);
975 }
976 
977 double SkyPoint::vHelioToVlsr(double vhelio, long double jd0)
978 {
979  return vhelio + vRSun(jd0);
980 }
981 
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 
1015 double SkyPoint::vGeocentric(double vhelio, long double jd0)
1016 {
1017  return vhelio - vREarth(jd0);
1018 }
1019 
1020 double SkyPoint::vGeoToVHelio(double vgeo, long double jd0)
1021 {
1022  return vgeo + vREarth(jd0);
1023 }
1024 
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 
1035 double SkyPoint::vTopoToVGeo(double vtopo, double vsite[3])
1036 {
1037  return vtopo + vRSite(vsite);
1038 }
1039 
1040 double SkyPoint::vTopocentric(double vgeo, double vsite[3])
1041 {
1042  return vgeo - vRSite(vsite);
1043 }
1044 
1045 bool SkyPoint::checkCircumpolar(const dms *gLat) const
1046 {
1047  return fabs(dec().Degrees()) > (90 - fabs(gLat->Degrees()));
1048 }
1049 
1051 {
1052  return refract(Alt, Options::useRefraction());
1053 }
1054 
1055 void SkyPoint::setAltRefracted(dms alt_apparent)
1056 {
1057  setAlt(unrefract(alt_apparent, Options::useRefraction()));
1058 }
1059 
1060 void SkyPoint::setAltRefracted(double alt_apparent)
1061 {
1062  setAlt(unrefract(alt_apparent, Options::useRefraction()));
1063 }
1064 
1065 double SkyPoint::refractionCorr(double alt)
1066 {
1067  return 1.02 / tan(dms::DegToRad * (alt + 10.3 / (alt + 5.11))) / 60;
1068 }
1069 
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.
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 
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 
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 
1136 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 
1144 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 
1153 {
1154  // N.B. Technically, we could use the form
1155  // cos(angle) = (sin(φ) - sin(h) sin(δ))/(cos(h) cos(δ))
1156  // where h = altitude, δ = declination, φ = latitude,
1157  // and then disambiguate the sign as
1158  // if (az().reduce() < 180°) angle = -angle;
1159  // However, acos(...) is inaccurate when cosine is nearly flat, i.e. near 0° and 180°.
1160  // It is therefore better to go through some extra pain to use atan2()
1161 
1162  // Therefore we use the form shown in Jean Meeus' book (14.1)
1163  dms HA = LST - ra();
1164  double tan_lat = lat.sin() / lat.cos();
1165  double angle = atan2( // Measured CW on sky map (See Meeus' Fig on Pg 99)
1166  HA.sin(),
1167  tan_lat * dec().cos() - HA.cos() * dec().sin()
1168  );
1169  return dms(angle / dms::DegToRad);
1170 }
1171 
1172 
1173 #ifndef KSTARS_LITE
1174 QDBusArgument &operator<<(QDBusArgument &argument, const SkyPoint &source)
1175 {
1176  argument.beginStructure();
1177  argument << source.ra().Hours() << source.dec().Degrees();
1178  argument.endStructure();
1179  return argument;
1180 }
1181 
1182 const QDBusArgument &operator>>(const QDBusArgument &argument, SkyPoint &dest)
1183 {
1184  double ra, dec;
1185  argument.beginStructure();
1186  argument >> ra >> dec;
1187  argument.endStructure();
1188  dest = SkyPoint(ra, dec);
1189  return argument;
1190 }
1191 #endif
1192 
static double unrefract(const double alt, bool conditional=true)
Remove refraction correction, depending on conditional.
Definition: skypoint.cpp:1091
const dms & alt() const
Definition: skypoint.h:281
SkyObject * findByName(const QString &name, bool exact=true) override
Search the children of this SkyComposite for a SkyObject whose name matches the argument.
double sin() const
Get the sine of this angle.
Definition: cachingdms.h:190
Extension of QDateTime for KStars KStarsDateTime can represent the date/time as a Julian Day,...
bool checkCircumpolar(const dms *gLat) const
Check if this point is circumpolar at the given geographic latitude.
Definition: skypoint.cpp:1045
void setAlt(dms alt)
Sets Alt, the Altitude.
Definition: skypoint.h:194
static constexpr double PI
PI is a const static member; it's public so that it can be used anywhere, as long as dms....
Definition: dms.h:385
bool bendlight()
Correct for the effect of "bending" of light around the sun for positions near the sun.
Definition: skypoint.cpp:425
double vREarth(long double jd0)
Computes the velocity of any object projected on the direction of the source.
Definition: skypoint.cpp:982
double p2b(int i1, int i2) const
Definition: ksnumbers.h:106
void setD(const double &x) override
Sets the angle in degrees supplied as a double.
Definition: cachingdms.h:59
virtual void setD(const double &x)
Sets floating-point value of angle, in degrees.
Definition: dms.h:179
static SkyPoint timeTransformed(const SkyPoint *p, const KStarsDateTime &dt, const GeoLocation *geo, const double hour=0)
returns a time-transformed SkyPoint.
Definition: skypoint.cpp:1121
static double refractionCorr(double alt)
Calculate refraction correction.
Definition: skypoint.cpp:1065
void setUsing_atan2(const double &y, const double &x)
Sets the angle using atan2()
Definition: cachingdms.cpp:62
void findEcliptic(const CachingDms *Obliquity, dms &EcLong, dms &EcLat)
Determine the Ecliptic coordinates of the SkyPoint, given the Julian Date.
Definition: skypoint.cpp:187
double p1b(int i1, int i2) const
Definition: ksnumbers.h:103
void GalacticToEquatorial1950(const dms *galLong, const dms *galLat)
Computes equatorial coordinates referred to 1950 from galactic ones referred to epoch B1950.
Definition: skypoint.cpp:754
Stores dms coordinates for a point in the sky. for converting between coordinate systems.
Definition: skypoint.h:44
a dms subclass that caches its sine and cosine values every time the angle is changed.
Definition: cachingdms.h:18
static constexpr double DegToRad
DegToRad is a const static member equal to the number of radians in one degree (dms::PI/180....
Definition: dms.h:390
double vTopocentric(double vgeo, double vsite[3])
Computes the radial velocity of a source referred to the observer site on the surface of the earth fr...
Definition: skypoint.cpp:1040
double minAlt(const dms &lat) const
Return the object's altitude at the lower culmination for the given latitude.
Definition: skypoint.cpp:1144
void SinCos(double &s, double &c) const
Get the sine and cosine together.
Definition: cachingdms.h:175
void SinCos(double &s, double &c) const
Compute Sine and Cosine of the angle simultaneously.
Definition: dms.h:444
CachingDms sunTrueLongitude() const
Definition: ksnumbers.h:71
SkyPoint Eterms(void)
Determine the E-terms of aberration In the past, the mean places of stars published in catalogs inclu...
Definition: skypoint.cpp:868
KCALENDARCORE_EXPORT QDataStream & operator>>(QDataStream &in, const KCalendarCore::Alarm::Ptr &)
dms parallacticAngle(const CachingDms &LST, const CachingDms &lat)
Return the Parallactic Angle.
Definition: skypoint.cpp:1152
double rearth() const
Definition: ksplanetbase.h:139
void subtractEterms(void)
Coordinates in the FK4 catalog include the effect of aberration due to the ellipticity of the orbit o...
Definition: skypoint.cpp:891
static double refract(const double alt, bool conditional=true)
Apply refraction correction to altitude, depending on conditional.
Definition: skypoint.cpp:1070
void J2000ToB1950(void)
Exact precession from epoch J2000 Besselian epoch 1950.
Definition: skypoint.cpp:819
double vGeocentric(double vhelio, long double jd)
Computes the radial velocity of a source referred to the center of the earth from the radial velocity...
Definition: skypoint.cpp:1015
void EquatorialToHorizontal(const CachingDms *LST, const CachingDms *lat)
Determine the (Altitude, Azimuth) coordinates of the SkyPoint from its (RA, Dec) coordinates,...
Definition: skypoint.cpp:77
KStarsDateTime addSecs(double s) const
dms altRefracted() const
Definition: skypoint.cpp:1050
long double getJD() const
Definition: ksnumbers.h:128
virtual void setH(const double &x)
Sets floating-point value of angle, in hours.
Definition: dms.h:210
double maxAlt(const dms &lat) const
Return the object's altitude at the upper culmination for the given latitude.
Definition: skypoint.cpp:1136
double vRSun(long double jd)
Computes the velocity of the Sun projected on the direction of the source.
Definition: skypoint.cpp:934
virtual void updateCoords(const KSNumbers *num, bool includePlanets=true, const CachingDms *lat=nullptr, const CachingDms *LST=nullptr, bool forceRecompute=false)
Determine the current coordinates (RA, Dec) from the catalog coordinates (RA0, Dec0),...
Definition: skypoint.cpp:582
Provides necessary information about the Sun.
Definition: kssun.h:23
QString i18n(const char *text, const TYPE &arg...)
Store several time-dependent astronomical quantities.
Definition: ksnumbers.h:42
const CachingDms & dec() const
Definition: skypoint.h:269
void setRadians(const double &a) override
Sets the angle in radians.
Definition: cachingdms.h:136
void precessFromAnyEpoch(long double jd0, long double jdf)
General case of precession.
Definition: skypoint.cpp:625
void setAltRefracted(dms alt_apparent)
Sets the apparent altitude, checking whether refraction corrections are enabled.
Definition: skypoint.cpp:1055
double dObliq() const
Definition: ksnumbers.h:81
double earthEccentricity() const
Definition: ksnumbers.h:77
void B1950ToJ2000(void)
Exact precession from Besselian epoch 1950 to epoch J2000.
Definition: skypoint.cpp:773
long double julianDay() const
Definition: ksnumbers.h:91
SkyPoint()
Default constructor.
Definition: skypoint.cpp:53
double dEcLong() const
Definition: ksnumbers.h:85
double vRSite(double vsite[3])
Computes the velocity of any object (observer's site) projected on the direction of the source.
Definition: skypoint.cpp:1025
double vGeoToVHelio(double vgeo, long double jd)
Computes the radial velocity of a source referred to the solar system barycenter from the velocity re...
Definition: skypoint.cpp:1020
void Equatorial1950ToGalactic(dms &galLong, dms &galLat)
Computes galactic coordinates from equatorial coordinates referred to epoch 1950.
Definition: skypoint.cpp:735
QTextStream & dec(QTextStream &stream)
void setFromEcliptic(const CachingDms *Obliquity, const dms &EcLong, const dms &EcLat)
Set the current (RA, Dec) coordinates of the SkyPoint, given pointers to its Ecliptic (Long,...
Definition: skypoint.cpp:201
double cos() const
Get the cosine of this angle.
Definition: cachingdms.h:204
void setH(const double &x) override
Sets the angle in hours, supplied as a double.
Definition: cachingdms.h:91
dms angularDistanceTo(const SkyPoint *sp, double *const positionAngle=nullptr) const
Computes the angular distance between two SkyObjects.
Definition: skypoint.cpp:899
double cos() const
Compute the Angle's Cosine.
Definition: dms.h:290
double p2(int i1, int i2) const
Definition: ksnumbers.h:100
static const double altCrit
Critical height for atmospheric refraction corrections.
Definition: skypoint.h:718
void set(const dms &r, const dms &d)
Sets RA, Dec and RA0, Dec0 according to arguments.
Definition: skypoint.cpp:63
SkyPoint moveAway(const SkyPoint &from, double dist) const
Find the SkyPoint obtained by moving distance dist (arcseconds) away from the givenSkyPoint.
Definition: skypoint.cpp:371
void precess(const KSNumbers *num)
Precess this SkyPoint's catalog coordinates to the epoch described by the given KSNumbers object.
Definition: skypoint.cpp:223
SkyMapComposite * skyComposite()
Definition: kstarsdata.h:165
double physicalSize() const
Definition: ksplanetbase.h:192
An angle, stored as degrees, but expressible in many ways.
Definition: dms.h:37
virtual void setRadians(const double &Rad)
Set angle according to the argument, in radians.
Definition: dms.h:333
SkyPoint deprecess(const KSNumbers *num, long double epoch=J2000)
Obtain a Skypoint with RA0 and Dec0 set from the RA, Dec of this skypoint.
Definition: skypoint.cpp:257
static dms findAltitude(const SkyPoint *p, const KStarsDateTime &dt, const GeoLocation *geo, const double hour=0)
Compute the altitude of a given skypoint hour hours from the given date/time.
Definition: skypoint.cpp:1111
double vTopoToVGeo(double vtopo, double vsite[3])
Computes the radial velocity of a source referred to the center of the Earth from the radial velocity...
Definition: skypoint.cpp:1035
const CachingDms & ra() const
Definition: skypoint.h:263
double p1(int i1, int i2) const
Definition: ksnumbers.h:97
const CachingDms * obliquity() const
Definition: ksnumbers.h:56
double radians() const
Express the angle in radians.
Definition: dms.h:325
QDebug operator<<(QDebug d, const QCPVector2D &vec)
Definition: qcustomplot.h:446
void apparentCoord(long double jd0, long double jdf)
Computes the apparent coordinates for this SkyPoint for any epoch, accounting for the effects of prec...
Definition: skypoint.cpp:700
double vHeliocentric(double vlsr, long double jd)
Computes the radial velocity of a source referred to the solar system barycenter from the radial velo...
Definition: skypoint.cpp:972
const double & Degrees() const
Definition: dms.h:141
bool checkBendLight()
Check if this sky point is close enough to the sun for gravitational lensing to be significant.
Definition: skypoint.cpp:399
void beginStructure()
void setDec(dms d)
Sets Dec, the current Declination.
Definition: skypoint.h:169
void endStructure()
void HorizontalToEquatorial(const dms *LST, const dms *lat)
Determine the (RA, Dec) coordinates of the SkyPoint from its (Altitude, Azimuth) coordinates,...
Definition: skypoint.cpp:143
void setRA(dms &r)
Sets RA, the current Right Ascension.
Definition: skypoint.h:144
void reduceToRange(enum dms::AngleRanges range)
Reduce this angle to the given range.
Definition: dms.cpp:438
CachingDms earthPerihelionLongitude() const
Definition: ksnumbers.h:74
double sin() const
Compute the Angle's Sine.
Definition: dms.h:258
double vHelioToVlsr(double vhelio, long double jd)
Computes the radial velocity of a source referred to the Local Standard of Rest, also known as VLSR f...
Definition: skypoint.cpp:977
void setUsing_asin(const double &sine)
Sets the angle using asin()
Definition: cachingdms.cpp:93
void aberrate(const KSNumbers *num, bool reverse=false)
Determine the effects of aberration for this SkyPoint.
Definition: skypoint.cpp:454
void nutate(const KSNumbers *num, const bool reverse=false)
Apply the effects of nutation to this SkyPoint.
Definition: skypoint.cpp:275
dms constAberr() const
Definition: ksnumbers.h:59
double Hours() const
Definition: dms.h:168
SkyPoint catalogueCoord(long double jdf)
Computes the J2000.0 catalogue coordinates for this SkyPoint using the epoch removing aberration,...
Definition: skypoint.cpp:710
void addEterms(void)
Coordinates in the FK4 catalog include the effect of aberration due to the ellipticity of the orbit o...
Definition: skypoint.cpp:883
Relevant data about an observing location on Earth.
Definition: geolocation.h:27
This file is part of the KDE documentation.
Documentation copyright © 1996-2023 The KDE developers.
Generated on Fri Jun 9 2023 04:02:26 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.