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 }