File indexing completed on 2024-03-24 15:17:56

0001 /*
0002     SPDX-FileCopyrightText: 2001-2005 Jason Harris <jharris@30doradus.org>
0003     SPDX-FileCopyrightText: 2004-2005 Pablo de Vicente <p.devicente@wanadoo.es>
0004 
0005     SPDX-License-Identifier: GPL-2.0-or-later
0006 */
0007 
0008 #include "skypoint.h"
0009 
0010 #include "dms.h"
0011 #include "ksnumbers.h"
0012 #include "kstarsdatetime.h"
0013 #include "kssun.h"
0014 #include "kstarsdata.h"
0015 #include "Options.h"
0016 #include "skyobject.h"
0017 #include "skycomponents/skymapcomposite.h"
0018 
0019 #include "config-kstars.h"
0020 
0021 #include <KLocalizedString>
0022 
0023 #include <QDebug>
0024 
0025 #include <cmath>
0026 
0027 
0028 
0029 // N.B. To use libnova depending on the availability of the package,
0030 // uncomment the following:
0031 /*
0032 #ifdef HAVE_LIBNOVA
0033 #define SKYPOINT_USE_LIBNOVA 1
0034 #endif
0035 */
0036 
0037 #ifdef SKYPOINT_USE_LIBNOVA
0038 #include <libnova/libnova.h>
0039 bool SkyPoint::implementationIsLibnova = true;
0040 #else
0041 bool SkyPoint::implementationIsLibnova = false;
0042 #endif
0043 
0044 #ifdef PROFILE_COORDINATE_CONVERSION
0045 #include <ctime> // For profiling, remove if not profiling.
0046 long unsigned SkyPoint::eqToHzCalls = 0;
0047 double SkyPoint::cpuTime_EqToHz     = 0.;
0048 #endif
0049 
0050 KSSun *SkyPoint::m_Sun         = nullptr;
0051 const double SkyPoint::altCrit = -1.0;
0052 
0053 SkyPoint::SkyPoint()
0054 {
0055     // Default constructor. Set nonsense values
0056     RA0.setD(-1);   // RA >= 0 always :-)
0057     Dec0.setD(180); // Dec is between -90 and 90 Degrees :-)
0058     RA            = RA0;
0059     Dec           = Dec0;
0060     lastPrecessJD = J2000; // By convention, we use J2000 coordinates
0061 }
0062 
0063 void SkyPoint::set(const dms &r, const dms &d)
0064 {
0065     RA0 = RA = r;
0066     Dec0 = Dec    = d;
0067     lastPrecessJD = J2000; // By convention, we use J2000 coordinates
0068 }
0069 
0070 void SkyPoint::EquatorialToHorizontal(const dms *LST, const dms *lat)
0071 {
0072     //    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!";
0073     CachingDms _LST(*LST), _lat(*lat);
0074     EquatorialToHorizontal(&_LST, &_lat);
0075 }
0076 
0077 void SkyPoint::EquatorialToHorizontal(const CachingDms *LST, const CachingDms *lat)
0078 {
0079 #ifdef PROFILE_COORDINATE_CONVERSION
0080     std::clock_t start = std::clock();
0081 #endif
0082     //Uncomment for spherical trig version
0083     double AltRad, AzRad;
0084     double sindec, cosdec, sinlat, coslat, sinHA, cosHA;
0085     double sinAlt, cosAlt;
0086 
0087     CachingDms HourAngle =
0088         (*LST) - ra(); // Using CachingDms subtraction operator to find cos/sin of HourAngle without calling sincos()
0089 
0090     lat->SinCos(sinlat, coslat);
0091     dec().SinCos(sindec, cosdec);
0092     HourAngle.SinCos(sinHA, cosHA);
0093 
0094     sinAlt = sindec * sinlat + cosdec * coslat * cosHA;
0095     AltRad = asin(sinAlt);
0096 
0097     cosAlt = sqrt(
0098                  1 -
0099                  sinAlt *
0100                  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.
0101     if (cosAlt == 0.)
0102         cosAlt = cos(AltRad);
0103 
0104     double arg = (sindec - sinlat * sinAlt) / (coslat * cosAlt);
0105     if (arg <= -1.0)
0106         AzRad = dms::PI;
0107     else if (arg >= 1.0)
0108         AzRad = 0.0;
0109     else
0110         AzRad = acos(arg);
0111 
0112     if (sinHA > 0.0 && AzRad != 0.0)
0113         AzRad = 2.0 * dms::PI - AzRad; // resolve acos() ambiguity
0114 
0115     Alt.setRadians(AltRad);
0116     Az.setRadians(AzRad);
0117 #ifdef PROFILE_COORDINATE_CONVERSION
0118     std::clock_t stop = std::clock();
0119     cpuTime_EqToHz += double(stop - start) / double(CLOCKS_PER_SEC); // Accumulate time in seconds
0120     ++eqToHzCalls;
0121 #endif
0122 
0123     // //Uncomment for XYZ version
0124     //      double xr, yr, zr, xr1, zr1, sa, ca;
0125     //  //Z-axis rotation by -LST
0126     //  dms a = dms( -1.0*LST->Degrees() );
0127     //  a.SinCos( sa, ca );
0128     //  xr1 = m_X*ca + m_Y*sa;
0129     //  yr  = -1.0*m_X*sa + m_Y*ca;
0130     //  zr1 = m_Z;
0131     //
0132     //  //Y-axis rotation by lat - 90.
0133     //  a = dms( lat->Degrees() - 90.0 );
0134     //  a.SinCos( sa, ca );
0135     //  xr = xr1*ca - zr1*sa;
0136     //  zr = xr1*sa + zr1*ca;
0137     //
0138     //  //FIXME: eventually, we will work with XYZ directly
0139     //  Alt.setRadians( asin( zr ) );
0140     //  Az.setRadians( atan2( yr, xr ) );
0141 }
0142 
0143 void SkyPoint::HorizontalToEquatorial(const dms *LST, const dms *lat)
0144 {
0145     double HARad, DecRad;
0146     double sinlat, coslat, sinAlt, cosAlt, sinAz, cosAz;
0147     double sinDec, cosDec;
0148 
0149     lat->SinCos(sinlat, coslat);
0150     alt().SinCos(sinAlt, cosAlt);
0151     Az.SinCos(sinAz, cosAz);
0152 
0153     sinDec = sinAlt * sinlat + cosAlt * coslat * cosAz;
0154     DecRad = asin(sinDec);
0155     cosDec = cos(DecRad);
0156     Dec.setRadians(DecRad);
0157 
0158     double x = (sinAlt - sinlat * sinDec) / (coslat * cosDec);
0159 
0160     //Under certain circumstances, x can be very slightly less than -1.0000, or slightly
0161     //greater than 1.0000, leading to a crash on acos(x).  However, the value isn't
0162     //*really* out of range; it's a kind of roundoff error.
0163     if (x < -1.0 && x > -1.000001)
0164         HARad = dms::PI;
0165     else if (x > 1.0 && x < 1.000001)
0166         HARad = 0.0;
0167     else if (x < -1.0)
0168     {
0169         //qWarning() << "Coordinate out of range.";
0170         HARad = dms::PI;
0171     }
0172     else if (x > 1.0)
0173     {
0174         //qWarning() << "Coordinate out of range.";
0175         HARad = 0.0;
0176     }
0177     else
0178         HARad = acos(x);
0179 
0180     if (sinAz > 0.0)
0181         HARad = 2.0 * dms::PI - HARad; // resolve acos() ambiguity
0182 
0183     RA.setRadians(LST->radians() - HARad);
0184     RA.reduceToRange(dms::ZERO_TO_2PI);
0185 }
0186 
0187 void SkyPoint::findEcliptic(const CachingDms *Obliquity, dms &EcLong, dms &EcLat)
0188 {
0189     double sinRA, cosRA, sinOb, cosOb, sinDec, cosDec;
0190     ra().SinCos(sinRA, cosRA);
0191     dec().SinCos(sinDec, cosDec);
0192     Obliquity->SinCos(sinOb, cosOb);
0193 
0194     double ycosDec        = sinRA * cosOb * cosDec + sinDec * sinOb;
0195     double ELongRad = atan2(ycosDec, cosDec * cosRA);
0196     EcLong.setRadians(ELongRad);
0197     EcLong.reduceToRange(dms::ZERO_TO_2PI);
0198     EcLat.setRadians(asin(sinDec * cosOb - cosDec * sinOb * sinRA)); // FIXME: Haversine
0199 }
0200 
0201 void SkyPoint::setFromEcliptic(const CachingDms *Obliquity, const dms &EcLong, const dms &EcLat)
0202 {
0203     double sinLong, cosLong, sinLat, cosLat, sinObliq, cosObliq;
0204     EcLong.SinCos(sinLong, cosLong);
0205     EcLat.SinCos(sinLat, cosLat);
0206     Obliquity->SinCos(sinObliq, cosObliq);
0207 
0208     // double sinDec = sinLat * cosObliq + cosLat * sinObliq * sinLong;
0209 
0210     double ycosLat = sinLong * cosObliq * cosLat - sinLat * sinObliq;
0211     //    double RARad =  atan2( y, cosLong );
0212     RA.setUsing_atan2(ycosLat, cosLong * cosLat);
0213     RA.reduceToRange(dms::ZERO_TO_2PI);
0214     // Dec.setUsing_asin(sinDec);
0215 
0216     // Use Haversine to set declination
0217     Dec.setRadians(dms::PI / 2.0 - 2.0 * asin(sqrt(0.5 * (
0218                        1.0 - sinLat * cosObliq
0219                        - cosLat * sinObliq * sinLong
0220                    ))));
0221 }
0222 
0223 void SkyPoint::precess(const KSNumbers *num)
0224 {
0225     double cosRA0, sinRA0, cosDec0, sinDec0;
0226     const Eigen::Matrix3d &precessionMatrix = num->p2();
0227     Eigen::Vector3d v, s;
0228 
0229     RA0.SinCos(sinRA0, cosRA0);
0230     Dec0.SinCos(sinDec0, cosDec0);
0231 
0232     s[0] = cosRA0 * cosDec0;
0233     s[1] = sinRA0 * cosDec0;
0234     s[2] = sinDec0;
0235 
0236     // NOTE: Rotation matrices are the fastest way to do rotations on
0237     // a vector. Quaternions need more multiplications. The rotation
0238     // matrix compensates in some sense by having more 'precomputed'
0239     // multiplications. The matrix elements seem to cache nicely, so
0240     // there isn't much overhead in accessing them.
0241 
0242     //Multiply P2 and s to get v, the vector representing the new coords.
0243     // for ( unsigned int i=0; i<3; ++i ) {
0244     //     v[i] = 0.0;
0245     //     for (uint j=0; j< 3; ++j) {
0246     //         v[i] += num->p2( j, i )*s[j];
0247     //     }
0248     // }
0249     v.noalias() = precessionMatrix * s;
0250 
0251     //Extract RA, Dec from the vector:
0252     RA.setUsing_atan2(v[1], v[0]);
0253     RA.reduceToRange(dms::ZERO_TO_2PI);
0254     Dec.setUsing_asin(v[2]);
0255 }
0256 
0257 SkyPoint SkyPoint::deprecess(const KSNumbers *num, long double epoch)
0258 {
0259     SkyPoint p1(RA, Dec);
0260     long double now = num->julianDay();
0261     p1.precessFromAnyEpoch(now, epoch);
0262     if ((std::isnan(RA0.Degrees()) || std::isnan(Dec0.Degrees())) ||
0263             (!std::isnan(Dec0.Degrees()) && fabs(Dec0.Degrees()) > 90.0))
0264     {
0265         // We have invalid RA0 and Dec0, so set them if epoch = J2000. Otherwise, do not touch.
0266         if (epoch == J2000L)
0267         {
0268             RA0  = p1.ra();
0269             Dec0 = p1.dec();
0270         }
0271     }
0272     return p1;
0273 }
0274 
0275 void SkyPoint::nutate(const KSNumbers *num, const bool reverse)
0276 {
0277     //Step 2: Nutation
0278     if (fabs(Dec.Degrees()) < 80.0) //approximate method
0279     {
0280         double cosRA, sinRA, cosDec, sinDec, tanDec;
0281         double cosOb, sinOb;
0282         double dRA, dDec;
0283 
0284         RA.SinCos(sinRA, cosRA);
0285         Dec.SinCos(sinDec, cosDec);
0286 
0287         tanDec = sinDec / cosDec;
0288 
0289         // Equ 23.1 in Jean Meeus' book
0290         // nut_ecliptic / num->obliquity() is called epsilon in Meeus
0291         // nut.longitude / num->dEcLong() is called delta psi in Meeus
0292         // nut.obliquity / num->dObliq() is called delta epsilon in Meeus
0293         // Meeus notes that these expressions are invalid if the star is
0294         // close to one of the celestial poles
0295 
0296         // N.B. These expressions are valid for FK5 coordinates
0297         // (presumably also valid for ICRS by extension), but not for
0298         // FK4. See the "Important Remark" on Page 152 of Jean Meeus'
0299         // book.
0300 
0301 #ifdef SKYPOINT_USE_LIBNOVA
0302         // code lifted from libnova ln_get_equ_nut, tailored to our needs
0303         // with the option to add or remove nutation
0304         struct ln_nutation nut;
0305         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
0306 
0307         double nut_ecliptic = ln_deg_to_rad(nut.ecliptic + nut.obliquity);
0308         sinOb = sin(nut_ecliptic);
0309         cosOb = cos(nut_ecliptic);
0310 
0311         dRA = nut.longitude * (cosOb + sinOb * sinRA * tanDec) - nut.obliquity * cosRA * tanDec;
0312         dDec = nut.longitude * (sinOb * cosRA) + nut.obliquity * sinRA;
0313 #else
0314         num->obliquity()->SinCos(sinOb, cosOb);
0315 
0316         // N.B. num->dEcLong() and num->dObliq() methods return in
0317         // degrees, whereby the corrections dRA and dDec are also in
0318         // degrees.
0319         dRA  = num->dEcLong() * (cosOb + sinOb * sinRA * tanDec) - num->dObliq() * cosRA * tanDec;
0320         dDec = num->dEcLong() * (sinOb * cosRA) + num->dObliq() * sinRA;
0321 #endif
0322         // the sign changed to remove nutation
0323         if (reverse)
0324         {
0325             dRA = -dRA;
0326             dDec = -dDec;
0327         }
0328         RA.setD(RA.Degrees() + dRA);
0329         Dec.setD(Dec.Degrees() + dDec);
0330     }
0331     else //exact method
0332     {
0333         // NOTE: Meeus declares that you must add Δψ to the ecliptic
0334         // longitude of the body to get a more accurate precession
0335         // result, but fails to emphasize that the NCP of the two
0336         // coordinates systems is different. The (RA, Dec) without
0337         // nutation computed, i.e. the mean place of the sky point is
0338         // referenced to the un-nutated geocentric frame (without the
0339         // obliquity correction), whereas the (RA, Dec) after nutation
0340         // is applied is referenced to the nutated geocentric
0341         // frame. This is more clearly explained in the "Explanatory
0342         // Supplement to the Astronomical Almanac" by
0343         // K. P. Seidelmann, which can be borrowed on the internet
0344         // archive as of this writing, see page 114:
0345         // https://archive.org/details/explanatorysuppl00pken
0346         //
0347         // The rotation matrix formulation in (3.222-3) and the figure
0348         // 3.222.1 make this clear
0349 
0350         // TODO apply reverse test above 80 degrees
0351         dms EcLong, EcLat;
0352         CachingDms obliquityWithoutNutation(*num->obliquity() - dms(num->dObliq()));
0353 
0354         if (reverse)
0355         {
0356             //Subtract dEcLong from the Ecliptic Longitude
0357             findEcliptic(num->obliquity(), EcLong, EcLat);
0358             dms newLong(EcLong.Degrees() - num->dEcLong());
0359             setFromEcliptic(&obliquityWithoutNutation, newLong, EcLat); // FIXME: Check
0360         }
0361         else
0362         {
0363             //Add dEcLong to the Ecliptic Longitude
0364             findEcliptic(&obliquityWithoutNutation, EcLong, EcLat);
0365             dms newLong(EcLong.Degrees() + num->dEcLong());
0366             setFromEcliptic(num->obliquity(), newLong, EcLat);
0367         }
0368     }
0369 }
0370 
0371 SkyPoint SkyPoint::moveAway(const SkyPoint &from, double dist) const
0372 {
0373     CachingDms lat1, dtheta;
0374 
0375     if (dist == 0.0)
0376     {
0377         qDebug() << Q_FUNC_INFO << "moveAway called with zero distance!";
0378         return *this;
0379     }
0380 
0381     double dst = fabs(dist * dms::DegToRad / 3600.0); // In radian
0382 
0383     // Compute the bearing angle w.r.t. the RA axis ("latitude")
0384     CachingDms dRA(ra() - from.ra());
0385     CachingDms dDec(dec() - from.dec());
0386     double bearing = atan2(dRA.sin() / dRA.cos(), dDec.sin()); // Do not use dRA = PI / 2!!
0387     //double bearing = atan2( dDec.radians() , dRA.radians() );
0388 
0389     //    double dir0 = (dist >= 0 ) ? bearing : bearing + dms::PI; // in radian
0390     double dir0   = bearing + std::signbit(dist) * dms::PI; // might be faster?
0391     double sinDst = sin(dst), cosDst = cos(dst);
0392 
0393     lat1.setUsing_asin(dec().sin() * cosDst + dec().cos() * sinDst * cos(dir0));
0394     dtheta.setUsing_atan2(sin(dir0) * sinDst * dec().cos(), cosDst - dec().sin() * lat1.sin());
0395 
0396     return SkyPoint(ra() + dtheta, lat1);
0397 }
0398 
0399 bool SkyPoint::checkBendLight()
0400 {
0401     // First see if we are close enough to the sun to bother about the
0402     // gravitational lensing effect. We correct for the effect at
0403     // least till b = 10 solar radii, where the effect is only about
0404     // 0.06".  Assuming min. sun-earth distance is 200 solar radii.
0405     static const dms maxAngle(1.75 * (30.0 / 200.0) / dms::DegToRad);
0406 
0407     if (!m_Sun)
0408     {
0409         SkyComposite *skycomopsite = KStarsData::Instance()->skyComposite();
0410 
0411         if (skycomopsite == nullptr)
0412             return false;
0413 
0414         m_Sun = dynamic_cast<KSSun *>(skycomopsite->findByName(i18n("Sun")));
0415 
0416         if (m_Sun == nullptr)
0417             return false;
0418     }
0419 
0420     // TODO: This can be optimized further. We only need a ballpark estimate of the distance to the sun to start with.
0421     return (fabs(angularDistanceTo(static_cast<const SkyPoint *>(m_Sun)).Degrees()) <=
0422             maxAngle.Degrees()); // NOTE: dynamic_cast is slow and not important here.
0423 }
0424 
0425 bool SkyPoint::bendlight()
0426 {
0427     // NOTE: This should be applied before aberration
0428     // NOTE: One must call checkBendLight() before unnecessarily calling this.
0429     // We correct for GR effects
0430 
0431     // NOTE: This code is buggy. The sun needs to be initialized to
0432     // the current epoch -- but we are not certain that this is the
0433     // case. We have, as of now, no way of telling if the sun is
0434     // initialized or not. If we initialize the sun here, we will be
0435     // slowing down the program rather substantially and potentially
0436     // introducing bugs. Therefore, we just ignore this problem, and
0437     // hope that whenever the user is interested in seeing the effects
0438     // of GR, we have the sun initialized correctly. This is usually
0439     // the case. When the sun is not correctly initialized, rearth()
0440     // is not computed, so we just assume it is nominally equal to 1
0441     // AU to get a reasonable estimate.
0442     Q_ASSERT(m_Sun);
0443     double corr_sec = 1.75 * m_Sun->physicalSize() /
0444                       ((std::isfinite(m_Sun->rearth()) ? m_Sun->rearth() : 1) * AU_KM *
0445                        angularDistanceTo(static_cast<const SkyPoint *>(m_Sun)).sin());
0446     Q_ASSERT(corr_sec > 0);
0447 
0448     SkyPoint sp = moveAway(*m_Sun, corr_sec);
0449     setRA(sp.ra());
0450     setDec(sp.dec());
0451     return true;
0452 }
0453 
0454 void SkyPoint::aberrate(const KSNumbers *num, bool reverse)
0455 {
0456 #ifdef SKYPOINT_USE_LIBNOVA
0457     ln_equ_posn pos { RA.Degrees(), Dec.Degrees() };
0458     ln_equ_posn abPos { 0, 0 };
0459     ln_get_equ_aber(&pos, num->julianDay(), &abPos);
0460     if (reverse)
0461     {
0462         RA.setD(RA.Degrees() * 2 - abPos.ra);
0463         Dec.setD(Dec.Degrees() * 2 - abPos.dec);
0464     }
0465     else
0466     {
0467         RA.setD(abPos.ra);
0468         Dec.setD(abPos.dec);
0469     }
0470 
0471 #else
0472     // N.B. These expressions are valid for FK5 coordinates
0473     // (presumably also valid for ICRS by extension), but not for
0474     // FK4. See the "Important Remark" on Page 152 of Jean Meeus'
0475     // book.
0476 
0477     // N.B. Even though Meeus does not say this explicitly, these
0478     // equations must not apply near the pole. Therefore, we fall-back
0479     // to the expressions provided by Meeus in ecliptic coordinates
0480     // (Equ 23.2) when we are near the pole.
0481 
0482     double K = num->constAberr().Degrees(); //constant of aberration
0483     double e = num->earthEccentricity(); // eccentricity of Earth's orbit
0484 
0485     if (fabs(Dec.Degrees()) < 80.0)
0486     {
0487 
0488         double cosRA, sinRA, cosDec, sinDec;
0489         double cosL, sinL, cosP, sinP;
0490         double cosOb, sinOb;
0491 
0492 
0493         RA.SinCos(sinRA, cosRA);
0494         Dec.SinCos(sinDec, cosDec);
0495 
0496         num->obliquity()->SinCos(sinOb, cosOb);
0497         // double tanOb = sinOb/cosOb;
0498 
0499         num->sunTrueLongitude().SinCos(sinL, cosL);
0500         num->earthPerihelionLongitude().SinCos(sinP, cosP);
0501 
0502         //Step 3: Aberration
0503         // These are the expressions given in Jean Meeus, Equation (23.3)
0504 
0505         // double dRA = -1.0 * K * ( cosRA * cosL * cosOb + sinRA * sinL )/cosDec
0506         //               + e * K * ( cosRA * cosP * cosOb + sinRA * sinP )/cosDec;
0507 
0508         // double dDec = -1.0 * K * ( cosL * cosOb * ( tanOb * cosDec - sinRA * sinDec ) + cosRA * sinDec * sinL )
0509         //                + e * K * ( cosP * cosOb * ( tanOb * cosDec - sinRA * sinDec ) + cosRA * sinDec * sinP );
0510 
0511         // However, we have factorized the expressions below by pulling
0512         // out common factors to make it more optimized, in case the
0513         // compiler fails to spot these optimizations.
0514 
0515         // N.B. I had failed to factor out the expressions correctly,
0516         // making mistakes, in c5e709bd91, which should now be
0517         // fixed. --asimha
0518 
0519         // FIXME: Check if the unit tests have sufficient coverage to
0520         // check this expression
0521 
0522         double dRA = (K / cosDec) * (
0523                          cosRA * cosOb * (e * cosP - cosL)
0524                          + sinRA * (e * sinP - sinL)
0525                      );
0526         double dDec = K * (
0527                           (sinOb * cosDec - cosOb * sinDec * sinRA) * (e * cosP - cosL)
0528                           + cosRA * sinDec * (e * sinP - sinL)
0529                       );
0530 
0531         // N.B. Meeus points out that the result has the same units as
0532         // K, so the corrections are in degrees.
0533 
0534         if (reverse)
0535         {
0536             dRA = -dRA;
0537             dDec = -dDec;
0538         }
0539         RA.setD(RA.Degrees() + dRA);
0540         Dec.setD(Dec.Degrees() + dDec);
0541     }
0542     else
0543     {
0544         dms EcLong, EcLat;
0545         double sinEcLat, cosEcLat;
0546         const auto &L = num->sunTrueLongitude();
0547         const auto &P = num->earthPerihelionLongitude();
0548 
0549         findEcliptic(num->obliquity(), EcLong, EcLat);
0550         EcLat.SinCos(sinEcLat, cosEcLat);
0551 
0552         double sin_L_minus_EcLong, cos_L_minus_EcLong;
0553         double sin_P_minus_EcLong, cos_P_minus_EcLong;
0554         (L - EcLong).SinCos(sin_L_minus_EcLong, cos_L_minus_EcLong);
0555         (P - EcLong).SinCos(sin_P_minus_EcLong, cos_P_minus_EcLong);
0556 
0557         // Equation (23.2) in Meeus
0558         // N.B. dEcLong, dEcLat are in degrees, because K is expressed in degrees.
0559         double dEcLong = (K / cosEcLat) * (e * cos_P_minus_EcLong - cos_L_minus_EcLong);
0560         double dEcLat = K * sinEcLat * (e * sin_P_minus_EcLong - sin_L_minus_EcLong);
0561 
0562         // Note: These are approximate corrections, so it is
0563         // appropriate to change their sign to reverse them to first
0564         // order in the corrections. As a result, the forward and
0565         // reverse calculations will not be exact inverses, but only
0566         // approximate inverses.
0567         if (reverse)
0568         {
0569             dEcLong = -dEcLong;
0570             dEcLat = -dEcLat;
0571         }
0572 
0573         // Update the ecliptic coordinates to their new values
0574         EcLong.setD(EcLong.Degrees() + dEcLong);
0575         EcLat.setD(EcLat.Degrees() + dEcLat);
0576         setFromEcliptic(num->obliquity(), EcLong, EcLat);
0577     }
0578 #endif
0579 }
0580 
0581 // Note: This method is one of the major rate determining factors in how fast the map pans / zooms in or out
0582 void SkyPoint::updateCoords(const KSNumbers *num, bool /*includePlanets*/, const CachingDms *lat, const CachingDms *LST,
0583                             bool forceRecompute)
0584 {
0585     //Correct the catalog coordinates for the time-dependent effects
0586     //of precession, nutation and aberration
0587     bool recompute, lens;
0588 
0589     // NOTE: The same short-circuiting checks are also implemented in
0590     // StarObject::JITUpdate(), even before calling
0591     // updateCoords(). While this is code-duplication, these bits of
0592     // code need to be really optimized, at least for stars. For
0593     // optimization purposes, the code is left duplicated in two
0594     // places. Please be wary of changing one without changing the
0595     // other.
0596 
0597     Q_ASSERT(std::isfinite(lastPrecessJD));
0598 
0599     if (Options::useRelativistic() && checkBendLight())
0600     {
0601         recompute = true;
0602         lens      = true;
0603     }
0604     else
0605     {
0606         recompute = (Options::alwaysRecomputeCoordinates() || forceRecompute ||
0607                      std::abs(lastPrecessJD - num->getJD()) >= 0.00069444); // Update once per solar minute
0608         lens      = false;
0609     }
0610     if (recompute)
0611     {
0612         precess(num);
0613         nutate(num);
0614         if (lens)
0615             bendlight(); // FIXME: Shouldn't we apply this on the horizontal coordinates?
0616         aberrate(num);
0617         lastPrecessJD = num->getJD();
0618         Q_ASSERT(std::isfinite(RA.Degrees()) && std::isfinite(Dec.Degrees()));
0619     }
0620 
0621     if (lat || LST)
0622         qWarning() << i18n("lat and LST parameters should only be used in KSPlanetBase objects.");
0623 }
0624 
0625 void SkyPoint::precessFromAnyEpoch(long double jd0, long double jdf)
0626 {
0627     double cosRA, sinRA, cosDec, sinDec;
0628     double v[3], s[3];
0629 
0630     RA  = RA0;
0631     Dec = Dec0; // Is this necessary?
0632 
0633     if (jd0 == jdf)
0634         return;
0635 
0636     RA.SinCos(sinRA, cosRA);
0637     Dec.SinCos(sinDec, cosDec);
0638 
0639     if (jd0 == B1950L)
0640     {
0641         B1950ToJ2000();
0642         jd0 = J2000L;
0643         RA.SinCos(sinRA, cosRA);
0644         Dec.SinCos(sinDec, cosDec);
0645     }
0646 
0647     if (jd0 != jdf)
0648     {
0649         // The original coordinate is referred to the FK5 system and
0650         // is NOT J2000.
0651         if (jd0 != J2000L)
0652         {
0653             //v is a column vector representing input coordinates.
0654             v[0] = cosRA * cosDec;
0655             v[1] = sinRA * cosDec;
0656             v[2] = sinDec;
0657 
0658             //Need to first precess to J2000.0 coords
0659             //s is the product of P1 and v; s represents the
0660             //coordinates precessed to J2000
0661             KSNumbers num(jd0);
0662             for (unsigned int i = 0; i < 3; ++i)
0663             {
0664                 s[i] = num.p1(0, i) * v[0] + num.p1(1, i) * v[1] + num.p1(2, i) * v[2];
0665             }
0666 
0667             //Input coords already in J2000, set s accordingly.
0668         }
0669         else
0670         {
0671             s[0] = cosRA * cosDec;
0672             s[1] = sinRA * cosDec;
0673             s[2] = sinDec;
0674         }
0675 
0676         if (jdf == B1950L)
0677         {
0678             RA.setRadians(atan2(s[1], s[0]));
0679             Dec.setRadians(asin(s[2]));
0680             J2000ToB1950();
0681 
0682             return;
0683         }
0684 
0685         KSNumbers num(jdf);
0686         for (unsigned int i = 0; i < 3; ++i)
0687         {
0688             v[i] = num.p2(0, i) * s[0] + num.p2(1, i) * s[1] + num.p2(2, i) * s[2];
0689         }
0690 
0691         RA.setUsing_atan2(v[1], v[0]);
0692         Dec.setUsing_asin(v[2]);
0693 
0694         RA.reduceToRange(dms::ZERO_TO_2PI);
0695 
0696         return;
0697     }
0698 }
0699 
0700 void SkyPoint::apparentCoord(long double jd0, long double jdf)
0701 {
0702     precessFromAnyEpoch(jd0, jdf);
0703     KSNumbers num(jdf);
0704     nutate(&num);
0705     if (Options::useRelativistic() && checkBendLight())
0706         bendlight();
0707     aberrate(&num);
0708 }
0709 
0710 SkyPoint SkyPoint::catalogueCoord(long double jdf)
0711 {
0712     KSNumbers num(jdf);
0713 
0714     // remove abberation
0715     aberrate(&num, true);
0716 
0717     // remove nutation
0718     nutate(&num, true);
0719 
0720     // remove precession
0721     // the start position needs to be in RA0,Dec0
0722     RA0 = RA;
0723     Dec0 = Dec;
0724     // from now to J2000
0725     precessFromAnyEpoch(jdf, static_cast<long double>(J2000));
0726     // the J2000 position is in RA,Dec, move to RA0, Dec0
0727     RA0 = RA;
0728     Dec0 = Dec;
0729     lastPrecessJD = J2000;
0730 
0731     SkyPoint sp(RA0, Dec0);
0732     return sp;
0733 }
0734 
0735 void SkyPoint::Equatorial1950ToGalactic(dms &galLong, dms &galLat)
0736 {
0737     double a = 192.25;
0738     double sinb, cosb, sina_RA, cosa_RA, sinDEC, cosDEC, tanDEC;
0739 
0740     dms c(303.0);
0741     dms b(27.4);
0742     tanDEC = tan(Dec.radians());
0743 
0744     b.SinCos(sinb, cosb);
0745     dms(a - RA.Degrees()).SinCos(sina_RA, cosa_RA);
0746     Dec.SinCos(sinDEC, cosDEC);
0747 
0748     galLong.setRadians(c.radians() - atan2(sina_RA, cosa_RA * sinb - tanDEC * cosb));
0749     galLong.reduceToRange(dms::ZERO_TO_2PI);
0750 
0751     galLat.setRadians(asin(sinDEC * sinb + cosDEC * cosb * cosa_RA));
0752 }
0753 
0754 void SkyPoint::GalacticToEquatorial1950(const dms *galLong, const dms *galLat)
0755 {
0756     double a = 123.0;
0757     double sinb, cosb, singLat, cosgLat, tangLat, singLong_a, cosgLong_a;
0758 
0759     dms c(12.25);
0760     dms b(27.4);
0761     tangLat = tan(galLat->radians());
0762     galLat->SinCos(singLat, cosgLat);
0763 
0764     dms(galLong->Degrees() - a).SinCos(singLong_a, cosgLong_a);
0765     b.SinCos(sinb, cosb);
0766 
0767     RA.setRadians(c.radians() + atan2(singLong_a, cosgLong_a * sinb - tangLat * cosb));
0768     RA.reduceToRange(dms::ZERO_TO_2PI);
0769 
0770     Dec.setRadians(asin(singLat * sinb + cosgLat * cosb * cosgLong_a));
0771 }
0772 
0773 void SkyPoint::B1950ToJ2000(void)
0774 {
0775     double cosRA, sinRA, cosDec, sinDec;
0776     //  double cosRA0, sinRA0, cosDec0, sinDec0;
0777     double v[3], s[3];
0778 
0779     // 1984 January 1 0h
0780     KSNumbers num(2445700.5L);
0781 
0782     // Eterms due to aberration
0783     addEterms();
0784     RA.SinCos(sinRA, cosRA);
0785     Dec.SinCos(sinDec, cosDec);
0786 
0787     // Precession from B1950 to J1984
0788     s[0] = cosRA * cosDec;
0789     s[1] = sinRA * cosDec;
0790     s[2] = sinDec;
0791     for (unsigned int i = 0; i < 3; ++i)
0792     {
0793         v[i] = num.p2b(0, i) * s[0] + num.p2b(1, i) * s[1] + num.p2b(2, i) * s[2];
0794     }
0795 
0796     // RA zero-point correction at 1984 day 1, 0h.
0797     RA.setRadians(atan2(v[1], v[0]));
0798     Dec.setRadians(asin(v[2]));
0799 
0800     RA.setH(RA.Hours() + 0.06390 / 3600.);
0801     RA.SinCos(sinRA, cosRA);
0802     Dec.SinCos(sinDec, cosDec);
0803 
0804     s[0] = cosRA * cosDec;
0805     s[1] = sinRA * cosDec;
0806     s[2] = sinDec;
0807 
0808     // Precession from 1984 to J2000.
0809 
0810     for (unsigned int i = 0; i < 3; ++i)
0811     {
0812         v[i] = num.p1(0, i) * s[0] + num.p1(1, i) * s[1] + num.p1(2, i) * s[2];
0813     }
0814 
0815     RA.setRadians(atan2(v[1], v[0]));
0816     Dec.setRadians(asin(v[2]));
0817 }
0818 
0819 void SkyPoint::J2000ToB1950(void)
0820 {
0821     double cosRA, sinRA, cosDec, sinDec;
0822     //  double cosRA0, sinRA0, cosDec0, sinDec0;
0823     double v[3], s[3];
0824 
0825     // 1984 January 1 0h
0826     KSNumbers num(2445700.5L);
0827 
0828     RA.SinCos(sinRA, cosRA);
0829     Dec.SinCos(sinDec, cosDec);
0830 
0831     s[0] = cosRA * cosDec;
0832     s[1] = sinRA * cosDec;
0833     s[2] = sinDec;
0834 
0835     // Precession from J2000 to 1984 day, 0h.
0836 
0837     for (unsigned int i = 0; i < 3; ++i)
0838     {
0839         v[i] = num.p2(0, i) * s[0] + num.p2(1, i) * s[1] + num.p2(2, i) * s[2];
0840     }
0841 
0842     RA.setRadians(atan2(v[1], v[0]));
0843     Dec.setRadians(asin(v[2]));
0844 
0845     // RA zero-point correction at 1984 day 1, 0h.
0846 
0847     RA.setH(RA.Hours() - 0.06390 / 3600.);
0848     RA.SinCos(sinRA, cosRA);
0849     Dec.SinCos(sinDec, cosDec);
0850 
0851     // Precession from B1950 to J1984
0852 
0853     s[0] = cosRA * cosDec;
0854     s[1] = sinRA * cosDec;
0855     s[2] = sinDec;
0856     for (unsigned int i = 0; i < 3; ++i)
0857     {
0858         v[i] = num.p1b(0, i) * s[0] + num.p1b(1, i) * s[1] + num.p1b(2, i) * s[2];
0859     }
0860 
0861     RA.setRadians(atan2(v[1], v[0]));
0862     Dec.setRadians(asin(v[2]));
0863 
0864     // Eterms due to aberration
0865     subtractEterms();
0866 }
0867 
0868 SkyPoint SkyPoint::Eterms(void)
0869 {
0870     double sd, cd, sinEterm, cosEterm;
0871     dms raTemp, raDelta, decDelta;
0872 
0873     Dec.SinCos(sd, cd);
0874     raTemp.setH(RA.Hours() + 11.25);
0875     raTemp.SinCos(sinEterm, cosEterm);
0876 
0877     raDelta.setH(0.0227 * sinEterm / (3600. * cd));
0878     decDelta.setD(0.341 * cosEterm * sd / 3600. + 0.029 * cd / 3600.);
0879 
0880     return SkyPoint(raDelta, decDelta);
0881 }
0882 
0883 void SkyPoint::addEterms(void)
0884 {
0885     SkyPoint spd = Eterms();
0886 
0887     RA  = RA + spd.ra();
0888     Dec = Dec + spd.dec();
0889 }
0890 
0891 void SkyPoint::subtractEterms(void)
0892 {
0893     SkyPoint spd = Eterms();
0894 
0895     RA  = RA - spd.ra();
0896     Dec = Dec - spd.dec();
0897 }
0898 
0899 dms SkyPoint::angularDistanceTo(const SkyPoint *sp, double *const positionAngle) const
0900 {
0901     // double dalpha = sp->ra().radians() - ra().radians() ;
0902     // double ddelta = sp->dec().radians() - dec().radians();
0903     CachingDms dalpha = sp->ra() - ra();
0904     CachingDms ddelta = sp->dec() - dec();
0905 
0906     // double sa = sin(dalpha/2.);
0907     // double sd = sin(ddelta/2.);
0908 
0909     // double hava = sa*sa;
0910     // double havd = sd*sd;
0911 
0912     // Compute the haversin directly:
0913     double hava = (1 - dalpha.cos()) / 2.;
0914     double havd = (1 - ddelta.cos()) / 2.;
0915 
0916     // Haversine law
0917     double aux = havd + (sp->dec().cos()) * dec().cos() * hava;
0918 
0919     dms angDist;
0920     angDist.setRadians(2. * fabs(asin(sqrt(aux))));
0921 
0922     if (positionAngle)
0923     {
0924         // Also compute the position angle of the line from this SkyPoint to sp
0925         //*positionAngle = acos( tan(-ddelta)/tan( angDist.radians() ) ); // FIXME: Might fail for large ddelta / zero angDist
0926         //if( -dalpha < 0 )
0927         //            *positionAngle = 2*M_PI - *positionAngle;
0928         *positionAngle =
0929             atan2f(dalpha.sin(), (dec().cos()) * tan(sp->dec().radians()) - (dec().sin()) * dalpha.cos()) * 180 / M_PI;
0930     }
0931     return angDist;
0932 }
0933 
0934 double SkyPoint::vRSun(long double jd0)
0935 {
0936     double ca, sa, cd, sd, vsun;
0937     double cosRA, sinRA, cosDec, sinDec;
0938 
0939     /* Sun apex (where the sun goes) coordinates */
0940 
0941     dms asun(270.9592); // Right ascention: 18h 3m 50.2s [J2000]
0942     dms dsun(30.00467); // Declination: 30^o 0' 16.8'' [J2000]
0943     vsun = 20.;         // [km/s]
0944 
0945     asun.SinCos(sa, ca);
0946     dsun.SinCos(sd, cd);
0947 
0948     /* We need an auxiliary SkyPoint since we need the
0949     * source referred to the J2000 equinox and we do not want to overwrite
0950     * the current values
0951     */
0952 
0953     SkyPoint aux;
0954     aux.set(RA0, Dec0);
0955 
0956     aux.precessFromAnyEpoch(jd0, J2000L);
0957 
0958     aux.ra().SinCos(sinRA, cosRA);
0959     aux.dec().SinCos(sinDec, cosDec);
0960 
0961     /* Computation is done performing the scalar product of a unitary vector
0962     in the direction of the source with the vector velocity of Sun, both being in the
0963     LSR reference system:
0964     Vlsr    = Vhel + Vsun.u_radial =>
0965     Vlsr    = Vhel + vsun(cos D cos A,cos D sen A,sen D).(cos d cos a,cos d sen a,sen d)
0966     Vhel    = Vlsr - Vsun.u_radial
0967     */
0968 
0969     return vsun * (cd * cosDec * (cosRA * ca + sa * sinRA) + sd * sinDec);
0970 }
0971 
0972 double SkyPoint::vHeliocentric(double vlsr, long double jd0)
0973 {
0974     return vlsr - vRSun(jd0);
0975 }
0976 
0977 double SkyPoint::vHelioToVlsr(double vhelio, long double jd0)
0978 {
0979     return vhelio + vRSun(jd0);
0980 }
0981 
0982 double SkyPoint::vREarth(long double jd0)
0983 {
0984     double sinRA, sinDec, cosRA, cosDec;
0985 
0986     /* u_radial = unitary vector in the direction of the source
0987         Vlsr    = Vhel + Vsun.u_radial
0988             = Vgeo + VEarth.u_radial + Vsun.u_radial  =>
0989 
0990         Vgeo    = (Vlsr -Vsun.u_radial) - VEarth.u_radial
0991             =  Vhel - VEarth.u_radial
0992             =  Vhel - (vx, vy, vz).(cos d cos a,cos d sen a,sen d)
0993     */
0994 
0995     /* We need an auxiliary SkyPoint since we need the
0996     * source referred to the J2000 equinox and we do not want to overwrite
0997     * the current values
0998     */
0999 
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 
1050 dms SkyPoint::altRefracted() const
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 
1121 SkyPoint SkyPoint::timeTransformed(const SkyPoint *p, const KStarsDateTime &dt, const GeoLocation *geo,
1122                                    const double hour)
1123 {
1124     Q_ASSERT(p);
1125     if (!p)
1126         return SkyPoint(NaN::d, NaN::d);
1127 
1128     // Jasem 2015-08-24 Using correct procedure to find altitude
1129     SkyPoint sp                   = *p; // make a copy
1130     KStarsDateTime targetDateTime = dt.addSecs(hour * 3600.0);
1131     dms LST                       = geo->GSTtoLST(targetDateTime.gst());
1132     sp.EquatorialToHorizontal(&LST, geo->lat());
1133     return sp;
1134 }
1135 
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 
1152 dms SkyPoint::parallacticAngle(const CachingDms &LST, const CachingDms &lat)
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