File indexing completed on 2024-04-14 03:43:10 

0001 /*
0002     SPDX-FileCopyrightText: 2001 Jason Harris <jharris@30doradus.org>
0003 
0004     SPDX-License-Identifier: GPL-2.0-or-later
0005 */
0006 
0007 #include "ksplanet.h"
0008 
0009 #include "ksnumbers.h"
0010 #include "ksutils.h"
0011 #include "ksfilereader.h"
0012 
0013 #include <cmath>
0014 #include <typeinfo>
0015 
0016 #include "kstars_debug.h"
0017 
0018 // Qt version calming
0019 #include <qtskipemptyparts.h>
0020 
0021 KSPlanet::OrbitDataManager KSPlanet::odm;
0022 
0023 KSPlanet::OrbitDataManager::OrbitDataManager()
0024 {
0025     //EMPTY
0026 }
0027 
0028 bool KSPlanet::OrbitDataManager::readOrbitData(const QString &fname, QVector<OrbitData> *vector)
0029 {
0030     QFile f;
0031 
0032     if (KSUtils::openDataFile(f, fname))
0033     {
0034         KSFileReader fileReader(f); // close file is included
0035         QStringList fields;
0036         while (fileReader.hasMoreLines())
0037         {
0038             fields = fileReader.readLine().split(' ', Qt::SkipEmptyParts);
0039 
0040             if (fields.size() == 3)
0041             {
0042                 double A = fields[0].toDouble();
0043                 double B = fields[1].toDouble();
0044                 double C = fields[2].toDouble();
0045                 vector->append(OrbitData(A, B, C));
0046             }
0047         }
0048     }
0049     else
0050     {
0051         return false;
0052     }
0053     return true;
0054 }
0055 
0056 bool KSPlanet::OrbitDataManager::loadData(KSPlanet::OrbitDataColl &odc, const QString &n)
0057 {
0058     QString fname, snum;
0059     QFile f;
0060     int nCount = 0;
0061     QString nl = n.toLower();
0062 
0063     if (hash.contains(nl))
0064     {
0065         odc = hash[nl];
0066         return true; //orbit data already loaded
0067     }
0068 
0069     //Create a new OrbitDataColl
0070     OrbitDataColl ret;
0071 
0072     //Ecliptic Longitude
0073     for (int i = 0; i < 6; ++i)
0074     {
0075         snum.setNum(i);
0076         fname = nl + ".L" + snum + ".vsop";
0077         if (readOrbitData(fname, &ret.Lon[i]))
0078             nCount++;
0079     }
0080 
0081     if (nCount == 0)
0082         return false;
0083 
0084     //Ecliptic Latitude
0085     for (int i = 0; i < 6; ++i)
0086     {
0087         snum.setNum(i);
0088         fname = nl + ".B" + snum + ".vsop";
0089         if (readOrbitData(fname, &ret.Lat[i]))
0090             nCount++;
0091     }
0092 
0093     if (nCount == 0)
0094         return false;
0095 
0096     //Heliocentric Distance
0097     for (int i = 0; i < 6; ++i)
0098     {
0099         snum.setNum(i);
0100         fname = nl + ".R" + snum + ".vsop";
0101         if (readOrbitData(fname, &ret.Dst[i]))
0102             nCount++;
0103     }
0104 
0105     if (nCount == 0)
0106         return false;
0107 
0108     hash[nl] = ret;
0109     odc      = hash[nl];
0110 
0111     return true;
0112 }
0113 
0114 KSPlanet::KSPlanet(const QString &s, const QString &imfile, const QColor &c, double pSize)
0115     : KSPlanetBase(s, imfile, c, pSize)
0116 {
0117 }
0118 
0119 KSPlanet::KSPlanet(int n) : KSPlanetBase()
0120 {
0121     switch (n)
0122     {
0123         case MERCURY:
0124             KSPlanetBase::init(i18n("Mercury"), "mercury", KSPlanetBase::planetColor[KSPlanetBase::MERCURY], 4879.4);
0125             break;
0126         case VENUS:
0127             KSPlanetBase::init(i18n("Venus"), "venus", KSPlanetBase::planetColor[KSPlanetBase::VENUS], 12103.6);
0128             break;
0129         case MARS:
0130             KSPlanetBase::init(i18n("Mars"), "mars", KSPlanetBase::planetColor[KSPlanetBase::MARS], 6792.4);
0131             break;
0132         case JUPITER:
0133             KSPlanetBase::init(i18n("Jupiter"), "jupiter", KSPlanetBase::planetColor[KSPlanetBase::JUPITER], 142984.);
0134             break;
0135         case SATURN:
0136             KSPlanetBase::init(i18n("Saturn"), "saturn", KSPlanetBase::planetColor[KSPlanetBase::SATURN], 120536.);
0137             break;
0138         case URANUS:
0139             KSPlanetBase::init(i18n("Uranus"), "uranus", KSPlanetBase::planetColor[KSPlanetBase::URANUS], 51118.);
0140             break;
0141         case NEPTUNE:
0142             KSPlanetBase::init(i18n("Neptune"), "neptune", KSPlanetBase::planetColor[KSPlanetBase::NEPTUNE], 49572.);
0143             break;
0144         default:
0145             qDebug() << Q_FUNC_INFO << "Error: Illegal identifier in KSPlanet constructor: " << n;
0146             break;
0147     }
0148 }
0149 
0150 KSPlanet *KSPlanet::clone() const
0151 {
0152     Q_ASSERT(typeid(this) == typeid(static_cast<const KSPlanet *>(this))); // Ensure we are not slicing a derived class
0153     return new KSPlanet(*this);
0154 }
0155 
0156 // TODO: Get rid of this dirty hack post KDE 4.2 release
0157 QString KSPlanet::untranslatedName() const
0158 {
0159     if (name() == i18n("Mercury"))
0160         return "Mercury";
0161     else if (name() == i18n("Venus"))
0162         return "Venus";
0163     else if (name() == i18n("Mars"))
0164         return "Mars";
0165     else if (name() == i18n("Jupiter"))
0166         return "Jupiter";
0167     else if (name() == i18n("Saturn"))
0168         return "Saturn";
0169     else if (name() == i18n("Uranus"))
0170         return "Uranus";
0171     else if (name() == i18n("Neptune"))
0172         return "Neptune";
0173     else if (name() == i18n("Earth"))
0174         return "Earth";
0175     else if (name() == i18n("Earth Shadow"))
0176         return "Earth Shadow";
0177     else if (name() == i18n("Moon"))
0178         return "Moon";
0179     else if (name() == i18n("Sun"))
0180         return "Sun";
0181     else
0182         return name();
0183 }
0184 
0185 //we don't need the reference to the ODC, so just give it a junk variable
0186 bool KSPlanet::loadData()
0187 {
0188     OrbitDataColl odc;
0189     return odm.loadData(odc, untranslatedName());
0190 }
0191 
0192 void KSPlanet::calcEcliptic(double Tau, EclipticPosition &epret) const
0193 {
0194     double sum[6];
0195     OrbitDataColl odc;
0196     double Tpow[6];
0197 
0198     Tpow[0] = 1.0;
0199     for (int i = 1; i < 6; ++i)
0200     {
0201         Tpow[i] = Tpow[i - 1] * Tau;
0202     }
0203 
0204     if (!odm.loadData(odc, untranslatedName()))
0205     {
0206         epret.longitude = dms(0.0);
0207         epret.latitude  = dms(0.0);
0208         epret.radius    = 0.0;
0209         qCWarning(KSTARS) << "Could not get data for name:" << name() << "(" << untranslatedName() << ")";
0210         return;
0211     }
0212 
0213     //Ecliptic Longitude
0214     for (int i = 0; i < 6; ++i)
0215     {
0216         sum[i] = 0.0;
0217         for (int j = 0; j < odc.Lon[i].size(); ++j)
0218         {
0219             sum[i] += odc.Lon[i][j].A * cos(odc.Lon[i][j].B + odc.Lon[i][j].C * Tau);
0220             /*
0221             qDebug() << Q_FUNC_INFO << "sum[" << i <<"] =" << sum[i] <<
0222                 " A = " << odc.Lon[i][j].A << " B = " << odc.Lon[i][j].B <<
0223                 " C = " << odc.Lon[i][j].C;
0224                 */
0225         }
0226         sum[i] *= Tpow[i];
0227         //qDebug() << Q_FUNC_INFO << name() << " : sum[" << i << "] = " << sum[i];
0228     }
0229 
0230     epret.longitude.setRadians(sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5]);
0231     epret.longitude.setD(epret.longitude.reduce().Degrees());
0232 
0233     //Compute Ecliptic Latitude
0234     for (uint i = 0; i < 6; ++i)
0235     {
0236         sum[i] = 0.0;
0237         for (int j = 0; j < odc.Lat[i].size(); ++j)
0238         {
0239             sum[i] += odc.Lat[i][j].A * cos(odc.Lat[i][j].B + odc.Lat[i][j].C * Tau);
0240         }
0241         sum[i] *= Tpow[i];
0242     }
0243 
0244     epret.latitude.setRadians(sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5]);
0245 
0246     //Compute Heliocentric Distance
0247     for (uint i = 0; i < 6; ++i)
0248     {
0249         sum[i] = 0.0;
0250         for (int j = 0; j < odc.Dst[i].size(); ++j)
0251         {
0252             sum[i] += odc.Dst[i][j].A * cos(odc.Dst[i][j].B + odc.Dst[i][j].C * Tau);
0253         }
0254         sum[i] *= Tpow[i];
0255     }
0256 
0257     epret.radius = sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5];
0258 
0259     /*
0260     qDebug() << Q_FUNC_INFO << name() << " pre: Lat = " << epret.latitude.toDMSString() << " Long = " <<
0261         epret.longitude.toDMSString() << " Dist = " << epret.radius;
0262     */
0263 }
0264 
0265 bool KSPlanet::findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *Earth)
0266 {
0267     if (Earth != nullptr)
0268     {
0269         double sinL, sinL0, sinB, sinB0;
0270         double cosL, cosL0, cosB, cosB0;
0271         double x = 0.0, y = 0.0, z = 0.0;
0272 
0273         double olddst = -1000;
0274         double dst    = 0;
0275 
0276         EclipticPosition trialpos;
0277 
0278         double jm = num->julianMillenia();
0279 
0280         Earth->ecLong().SinCos(sinL0, cosL0);
0281         Earth->ecLat().SinCos(sinB0, cosB0);
0282 
0283         double eX = Earth->rsun() * cosB0 * cosL0;
0284         double eY = Earth->rsun() * cosB0 * sinL0;
0285         double eZ = Earth->rsun() * sinB0;
0286 
0287         bool once = true;
0288         while (fabs(dst - olddst) > .001)
0289         {
0290             calcEcliptic(jm, trialpos);
0291 
0292             // We store the heliocentric ecliptic coordinates the first time they are computed.
0293             if (once)
0294             {
0295                 helEcPos = trialpos;
0296                 once     = false;
0297             }
0298 
0299             olddst = dst;
0300 
0301             trialpos.longitude.SinCos(sinL, cosL);
0302             trialpos.latitude.SinCos(sinB, cosB);
0303 
0304             x = trialpos.radius * cosB * cosL - eX;
0305             y = trialpos.radius * cosB * sinL - eY;
0306             z = trialpos.radius * sinB - eZ;
0307 
0308             //distance from Earth
0309             dst = sqrt(x * x + y * y + z * z);
0310 
0311             //The light-travel time delay, in millenia
0312             //0.0057755183 is the inverse speed of light,
0313             //in days/AU
0314             double delay = (.0057755183 * dst) / 365250.0;
0315 
0316             jm = num->julianMillenia() - delay;
0317         }
0318 
0319         ep.longitude.setRadians(atan2(y, x));
0320         ep.longitude.reduce();
0321         ep.latitude.setRadians(atan2(z, sqrt(x * x + y * y)));
0322         setRsun(trialpos.radius);
0323         setRearth(dst);
0324 
0325         EclipticToEquatorial(num->obliquity());
0326 
0327         setRA0(ra());
0328         setDec0(dec());
0329         apparentCoord(J2000, lastPrecessJD);
0330 
0331         //nutate(num);
0332         //aberrate(num);
0333     }
0334     else
0335     {
0336         calcEcliptic(num->julianMillenia(), ep);
0337         helEcPos = ep;
0338     }
0339 
0340     //determine the position angle
0341     findPA(num);
0342 
0343     return true;
0344 }
0345 
0346 void KSPlanet::findMagnitude(const KSNumbers *num)
0347 {
0348     double cosDec, sinDec;
0349     dec().SinCos(cosDec, sinDec);
0350 
0351     /* Computation of the visual magnitude (V band) of the planet.
0352     * Algorithm provided by Pere Planesas (Observatorio Astronomico Nacional)
0353     * It has some similarity to J. Meeus algorithm in Astronomical Algorithms, Chapter 40.
0354     * */
0355 
0356     // Initialized to the faintest magnitude observable with the HST
0357     float magnitude = 30;
0358 
0359     double param = 5 * log10(rsun() * rearth());
0360     double phase = this->phase().Degrees();
0361     double f1    = phase / 100.;
0362 
0363     if (name() == i18n("Mercury"))
0364     {
0365         if (phase > 150.)
0366             f1 = 1.5;
0367         magnitude = -0.36 + param + 3.8 * f1 - 2.73 * f1 * f1 + 2 * f1 * f1 * f1;
0368     }
0369     else if (name() == i18n("Venus"))
0370     {
0371         magnitude = -4.29 + param + 0.09 * f1 + 2.39 * f1 * f1 - 0.65 * f1 * f1 * f1;
0372     }
0373     else if (name() == i18n("Mars"))
0374     {
0375         magnitude = -1.52 + param + 0.016 * phase;
0376     }
0377     else if (name() == i18n("Jupiter"))
0378     {
0379         magnitude = -9.25 + param + 0.005 * phase;
0380     }
0381     else if (name() == i18n("Saturn"))
0382     {
0383         double T     = num->julianCenturies();
0384         double a0    = (40.66 - 4.695 * T) * dms::PI / 180.;
0385         double d0    = (83.52 + 0.403 * T) * dms::PI / 180.;
0386         double sinx  = -cos(d0) * cosDec * cos(a0 - ra().radians());
0387         sinx         = fabs(sinx - sin(d0) * sinDec);
0388         double rings = -2.6 * sinx + 1.25 * sinx * sinx;
0389         magnitude    = -8.88 + param + 0.044 * phase + rings;
0390     }
0391     else if (name() == i18n("Uranus"))
0392     {
0393         magnitude = -7.19 + param + 0.0028 * phase;
0394     }
0395     else if (name() == i18n("Neptune"))
0396     {
0397         magnitude = -6.87 + param;
0398     }
0399     setMag(magnitude);
0400 }
0401 
0402 SkyObject::UID KSPlanet::getUID() const
0403 {
0404     SkyObject::UID n;
0405     if (name() == i18n("Mercury"))
0406     {
0407         n = 1;
0408     }
0409     else if (name() == i18n("Venus"))
0410     {
0411         n = 2;
0412     }
0413     else if (name() == i18n("Earth"))
0414     {
0415         n = 3;
0416     }
0417     else if (name() == i18n("Mars"))
0418     {
0419         n = 4;
0420     }
0421     else if (name() == i18n("Jupiter"))
0422     {
0423         n = 5;
0424     }
0425     else if (name() == i18n("Saturn"))
0426     {
0427         n = 6;
0428     }
0429     else if (name() == i18n("Uranus"))
0430     {
0431         n = 7;
0432     }
0433     else if (name() == i18n("Neptune"))
0434     {
0435         n = 8;
0436     }
0437     else
0438     {
0439         return SkyObject::invalidUID;
0440     }
0441     return solarsysUID(UID_SOL_BIGOBJ) | n;
0442 }