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

0001 /*
0002     SPDX-FileCopyrightText: 2001 Jason Harris <jharris@30doradus.org>
0003 
0004     SPDX-License-Identifier: GPL-2.0-or-later
0005 */
0006 
0007 #include "kssun.h"
0008 
0009 #include <typeinfo>
0010 
0011 #include <cmath>
0012 
0013 #include "ksnumbers.h"
0014 #include "kstarsdata.h"
0015 #include "kstarsdatetime.h"
0016 
0017 KSSun::KSSun() : KSPlanet(i18n("Sun"), "sun", Qt::yellow, 1392000. /*diameter in km*/)
0018 {
0019     setMag(-26.73);
0020 }
0021 
0022 KSSun *KSSun::clone() const
0023 {
0024     Q_ASSERT(typeid(this) == typeid(static_cast<const KSSun *>(this))); // Ensure we are not slicing a derived class
0025     return new KSSun(*this);
0026 }
0027 
0028 bool KSSun::loadData()
0029 {
0030     OrbitDataColl odc;
0031     return (odm.loadData(odc, "earth") != 0);
0032 }
0033 
0034 // We don't need to do anything here
0035 void KSSun::findMagnitude(const KSNumbers *)
0036 {
0037 }
0038 
0039 bool KSSun::findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *Earth)
0040 {
0041     if (Earth)
0042     {
0043         //
0044         // For the precision we need, the earth's orbit is circular.
0045         // So don't bother to iterate like KSPlanet does. Just subtract
0046         // The current delay and recompute (once).
0047         //
0048 
0049         //The light-travel time delay, in millenia
0050         //0.0057755183 is the inverse speed of light, in days/AU
0051         double delay = (.0057755183 * Earth->rsun()) / 365250.0;
0052         //
0053         // MHH 2002-02-04 I don't like this. But it avoids code duplication.
0054         // Maybe we can find a better way.
0055         //
0056         const KSPlanet *pEarth = static_cast<const KSPlanet *>(Earth);
0057         EclipticPosition trialpos;
0058         pEarth->calcEcliptic(num->julianMillenia() - delay, trialpos);
0059 
0060         setEcLong((trialpos.longitude + dms(180.0)).reduce());
0061         setEcLat(-trialpos.latitude);
0062 
0063         setRearth(Earth->rsun());
0064     }
0065     else
0066     {
0067         double sum[6];
0068         dms EarthLong, EarthLat; //heliocentric coords of Earth
0069         OrbitDataColl odc;
0070         double T = num->julianMillenia(); //Julian millenia since J2000
0071         double Tpow[6];
0072 
0073         Tpow[0] = 1.0;
0074         for (int i = 1; i < 6; ++i)
0075         {
0076             Tpow[i] = Tpow[i - 1] * T;
0077         }
0078         //First, find heliocentric coordinates
0079 
0080         if (!odm.loadData(odc, "earth"))
0081             return false;
0082 
0083         //Ecliptic Longitude
0084         for (int i = 0; i < 6; ++i)
0085         {
0086             sum[i] = 0.0;
0087             for (int j = 0; j < odc.Lon[i].size(); ++j)
0088             {
0089                 sum[i] += odc.Lon[i][j].A * cos(odc.Lon[i][j].B + odc.Lon[i][j].C * T);
0090             }
0091             sum[i] *= Tpow[i];
0092             //qDebug() << Q_FUNC_INFO << name() << " : sum[" << i << "] = " << sum[i];
0093         }
0094 
0095         EarthLong.setRadians(sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5]);
0096         EarthLong = EarthLong.reduce();
0097 
0098         //Compute Ecliptic Latitude
0099         for (int i = 0; i < 6; ++i)
0100         {
0101             sum[i] = 0.0;
0102             for (int j = 0; j < odc.Lat[i].size(); ++j)
0103             {
0104                 sum[i] += odc.Lat[i][j].A * cos(odc.Lat[i][j].B + odc.Lat[i][j].C * T);
0105             }
0106             sum[i] *= Tpow[i];
0107         }
0108 
0109         EarthLat.setRadians(sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5]);
0110 
0111         //Compute Heliocentric Distance
0112         for (int i = 0; i < 6; ++i)
0113         {
0114             sum[i] = 0.0;
0115             for (int j = 0; j < odc.Dst[i].size(); ++j)
0116             {
0117                 sum[i] += odc.Dst[i][j].A * cos(odc.Dst[i][j].B + odc.Dst[i][j].C * T);
0118             }
0119             sum[i] *= Tpow[i];
0120         }
0121 
0122         ep.radius = sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5];
0123         setRearth(ep.radius);
0124 
0125         setEcLong((EarthLong + dms(180.0)).reduce());
0126         setEcLat(-EarthLat);
0127     }
0128 
0129     //Finally, convert Ecliptic coords to Ra, Dec.  Ecliptic latitude is zero, by definition
0130     EclipticToEquatorial(num->obliquity());
0131 
0132     nutate(num);
0133 
0134     // Store in RA0 and Dec0, the unaberrated coordinates
0135     setRA0(ra());
0136     setDec0(dec());
0137 
0138     precess(num);
0139 
0140     aberrate(num);
0141 
0142     // We obtain the apparent geocentric ecliptic coordinates. That is, after
0143     // nutation and aberration have been applied.
0144     EquatorialToEcliptic(num->obliquity());
0145 
0146     //Determine the position angle
0147     findPA(num);
0148 
0149     return true;
0150 }
0151 
0152 SkyObject::UID KSSun::getUID() const
0153 {
0154     return solarsysUID(UID_SOL_BIGOBJ) | 0;
0155 }