Kstars

skypoint.cpp
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>
39bool SkyPoint::implementationIsLibnova = true;
40#else
41bool SkyPoint::implementationIsLibnova = false;
42#endif
43
44#ifdef PROFILE_COORDINATE_CONVERSION
45#include <ctime> // For profiling, remove if not profiling.
46long unsigned SkyPoint::eqToHzCalls = 0;
47double SkyPoint::cpuTime_EqToHz = 0.;
48#endif
49
50KSSun *SkyPoint::m_Sun = nullptr;
51const 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
63void 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
70void 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
143void 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
187void 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
201void 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
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
257SkyPoint 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
275void 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
371SkyPoint 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
454void 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
582void 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
625void 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
700void 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
754void 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
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
899dms 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
934double 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
972double SkyPoint::vHeliocentric(double vlsr, long double jd0)
973{
974 return vlsr - vRSun(jd0);
975}
976
977double SkyPoint::vHelioToVlsr(double vhelio, long double jd0)
978{
979 return vhelio + vRSun(jd0);
980}
981
982double 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
1015double SkyPoint::vGeocentric(double vhelio, long double jd0)
1016{
1017 return vhelio - vREarth(jd0);
1018}
1019
1020double SkyPoint::vGeoToVHelio(double vgeo, long double jd0)
1021{
1022 return vgeo + vREarth(jd0);
1023}
1024
1025double 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
1035double SkyPoint::vTopoToVGeo(double vtopo, double vsite[3])
1036{
1037 return vtopo + vRSite(vsite);
1038}
1039
1040double SkyPoint::vTopocentric(double vgeo, double vsite[3])
1041{
1042 return vgeo - vRSite(vsite);
1043}
1044
1045bool 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
1056{
1057 setAlt(unrefract(alt_apparent, Options::useRefraction()));
1058}
1059
1060void SkyPoint::setAltRefracted(double alt_apparent)
1061{
1062 setAlt(unrefract(alt_apparent, Options::useRefraction()));
1063}
1064
1066{
1067 return 1.02 / tan(dms::DegToRad * (alt + 10.3 / (alt + 5.11))) / 60;
1068}
1069
1070double 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.
1091double 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
1111dms 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
1136double 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
1144double 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
1174QDBusArgument &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
1182const 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
a dms subclass that caches its sine and cosine values every time the angle is changed.
Definition cachingdms.h:19
double cos() const
Get the cosine of this angle.
Definition cachingdms.h:204
void setUsing_atan2(const double &y, const double &x)
Sets the angle using atan2()
void setRadians(const double &a) override
Sets the angle in radians.
Definition cachingdms.h:136
double sin() const
Get the sine of this angle.
Definition cachingdms.h:190
void SinCos(double &s, double &c) const
Get the sine and cosine together.
Definition cachingdms.h:175
void setUsing_asin(const double &sine)
Sets the angle using asin()
void setD(const double &x) override
Sets the angle in degrees supplied as a double.
Definition cachingdms.h:59
void setH(const double &x) override
Sets the angle in hours, supplied as a double.
Definition cachingdms.h:91
Contains all relevant information for specifying a location on Earth: City Name, State/Province name,...
Definition geolocation.h:28
There are several time-dependent values used in position calculations, that are not specific to an ob...
Definition ksnumbers.h:43
double p1(int i1, int i2) const
Definition ksnumbers.h:97
dms constAberr() const
Definition ksnumbers.h:59
double p1b(int i1, int i2) const
Definition ksnumbers.h:103
CachingDms earthPerihelionLongitude() const
Definition ksnumbers.h:74
CachingDms sunTrueLongitude() const
Definition ksnumbers.h:71
const CachingDms * obliquity() const
Definition ksnumbers.h:56
long double getJD() const
Definition ksnumbers.h:128
double dEcLong() const
Definition ksnumbers.h:85
double dObliq() const
Definition ksnumbers.h:81
double p2(int i1, int i2) const
Definition ksnumbers.h:100
double earthEccentricity() const
Definition ksnumbers.h:77
long double julianDay() const
Definition ksnumbers.h:91
double p2b(int i1, int i2) const
Definition ksnumbers.h:106
double physicalSize() const
double rearth() const
Child class of KSPlanetBase; encapsulates information about the Sun.
Definition kssun.h:24
SkyMapComposite * skyComposite()
Definition kstarsdata.h:166
Extension of QDateTime for KStars KStarsDateTime can represent the date/time as a Julian Day,...
KStarsDateTime addSecs(double s) const
SkyComposite is a kind of container class for SkyComponent objects.
SkyObject * findByName(const QString &name, bool exact=true) override
Search the children of this SkyComposite for a SkyObject whose name matches the argument.
The sky coordinates of a point in the sky.
Definition skypoint.h:45
double maxAlt(const dms &lat) const
Return the object's altitude at the upper culmination for the given latitude.
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
static double refract(const double alt, bool conditional=true)
Apply refraction correction to altitude, depending on conditional.
static const double altCrit
Critical height for atmospheric refraction corrections.
Definition skypoint.h:727
void aberrate(const KSNumbers *num, bool reverse=false)
Determine the effects of aberration for this SkyPoint.
Definition skypoint.cpp:454
void precessFromAnyEpoch(long double jd0, long double jdf)
General case of precession.
Definition skypoint.cpp:625
static double refractionCorr(double alt)
Calculate refraction correction.
const CachingDms & dec() const
Definition skypoint.h:269
void setAltRefracted(dms alt_apparent)
Sets the apparent altitude, checking whether refraction corrections are enabled.
double minAlt(const dms &lat) const
Return the object's altitude at the lower culmination for the given latitude.
double vREarth(long double jd0)
Computes the velocity of any object projected on the direction of the source.
Definition skypoint.cpp:982
void setDec(dms d)
Sets Dec, the current Declination.
Definition skypoint.h:169
bool checkCircumpolar(const dms *gLat) const
Check if this point is circumpolar at the given geographic latitude.
void setRA(dms &r)
Sets RA, the current Right Ascension.
Definition skypoint.h:144
const CachingDms & ra() const
Definition skypoint.h:263
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.
dms altRefracted() const
void precess(const KSNumbers *num)
Precess this SkyPoint's catalog coordinates to the epoch described by the given KSNumbers object.
Definition skypoint.cpp:223
static SkyPoint timeTransformed(const SkyPoint *p, const KStarsDateTime &dt, const GeoLocation *geo, const double hour=0)
returns a time-transformed SkyPoint.
void findEcliptic(const CachingDms *Obliquity, dms &EcLong, dms &EcLat)
Determine the Ecliptic coordinates of the SkyPoint, given the Julian Date.
Definition skypoint.cpp:187
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
double vGeoToVHelio(double vgeo, long double jd)
Computes the radial velocity of a source referred to the solar system barycenter from the velocity re...
void nutate(const KSNumbers *num, const bool reverse=false)
Apply the effects of nutation to this SkyPoint.
Definition skypoint.cpp:275
dms angularDistanceTo(const SkyPoint *sp, double *const positionAngle=nullptr) const
Computes the angular distance between two SkyObjects.
Definition skypoint.cpp:899
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...
void B1950ToJ2000(void)
Exact precession from Besselian epoch 1950 to epoch J2000.
Definition skypoint.cpp:773
bool bendlight()
Correct for the effect of "bending" of light around the sun for positions near the sun.
Definition skypoint.cpp:425
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
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 EquatorialToHorizontal(const CachingDms *LST, const CachingDms *lat)
Determine the (Altitude, Azimuth) coordinates of the SkyPoint from its (RA, Dec) coordinates,...
Definition skypoint.cpp:77
double vRSite(double vsite[3])
Computes the velocity of any object (observer's site) projected on the direction of the source.
static double unrefract(const double alt, bool conditional=true)
Remove refraction correction, depending on conditional.
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
void set(const dms &r, const dms &d)
Sets RA, Dec and RA0, Dec0 according to arguments.
Definition skypoint.cpp:63
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...
void setAlt(dms alt)
Sets Alt, the Altitude.
Definition skypoint.h:194
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 dms & alt() const
Definition skypoint.h:281
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
bool checkBendLight()
Check if this sky point is close enough to the sun for gravitational lensing to be significant.
Definition skypoint.cpp:399
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
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
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...
void J2000ToB1950(void)
Exact precession from epoch J2000 Besselian epoch 1950.
Definition skypoint.cpp:819
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
SkyPoint()
Default constructor.
Definition skypoint.cpp:53
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
SkyPoint catalogueCoord(long double jdf)
Computes the J2000.0 catalogue coordinates for this SkyPoint using the epoch removing aberration,...
Definition skypoint.cpp:710
void Equatorial1950ToGalactic(dms &galLong, dms &galLat)
Computes galactic coordinates from equatorial coordinates referred to epoch 1950.
Definition skypoint.cpp:735
dms parallacticAngle(const CachingDms &LST, const CachingDms &lat)
Return the Parallactic Angle.
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
An angle, stored as degrees, but expressible in many ways.
Definition dms.h:38
double cos() const
Compute the Angle's Cosine.
Definition dms.h:290
double Hours() const
Definition dms.h:168
void reduceToRange(enum dms::AngleRanges range)
Reduce this angle to the given range.
Definition dms.cpp:446
virtual void setH(const double &x)
Sets floating-point value of angle, in hours.
Definition dms.h:210
void SinCos(double &s, double &c) const
Compute Sine and Cosine of the angle simultaneously.
Definition dms.h:447
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
double radians() const
Express the angle in radians.
Definition dms.h:325
virtual void setRadians(const double &Rad)
Set angle according to the argument, in radians.
Definition dms.h:333
virtual void setD(const double &x)
Sets floating-point value of angle, in degrees.
Definition dms.h:179
const double & Degrees() const
Definition dms.h:141
double sin() const
Compute the Angle's Sine.
Definition dms.h:258
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
QString i18n(const char *text, const TYPE &arg...)
KCALENDARCORE_EXPORT QDataStream & operator>>(QDataStream &in, const KCalendarCore::Alarm::Ptr &)
void beginStructure()
void endStructure()
QTextStream & dec(QTextStream &stream)
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Mon Nov 4 2024 16:38:44 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.