kstars
ksplanet.cpp
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include <math.h>
00019 #include <kdebug.h>
00020 #include <qfile.h>
00021
00022 #include "ksplanet.h"
00023 #include "ksnumbers.h"
00024 #include "ksutils.h"
00025 #include "ksfilereader.h"
00026
00027 KSPlanet::OrbitDataManager KSPlanet::odm;
00028
00029 KSPlanet::OrbitDataColl::OrbitDataColl() {
00030
00031 for (int i=0; i<6; i++) {
00032 Lon[i].setAutoDelete(true);
00033 Lat[i].setAutoDelete(true);
00034 Dst[i].setAutoDelete(true);
00035 }
00036 }
00037
00038
00039 KSPlanet::OrbitDataManager::OrbitDataManager() : dict(31, true) {
00040
00041 dict.setAutoDelete(true);
00042 }
00043
00044 bool KSPlanet::OrbitDataManager::readOrbitData(QString fname,
00045 QPtrVector<KSPlanet::OrbitData> *vector) {
00046 QString line;
00047 QFile f;
00048 double A, B, C;
00049
00050 QPtrList<OrbitData> DData;
00051
00052 if ( KSUtils::openDataFile( f, fname ) ) {
00053 KSFileReader fileReader( f );
00054 while ( fileReader.hasMoreLines() ) {
00055 line = fileReader.readLine();
00056 QTextIStream instream( &line );
00057 instream >> A >> B >> C;
00058 DData.append(new OrbitData(A, B, C));
00059
00060 }
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072 } else {
00073 return false;
00074 }
00075
00076 DData.toVector(vector);
00077
00078 return true;
00079 }
00080
00081 KSPlanet::OrbitDataColl *KSPlanet::OrbitDataManager::loadData(QString n) {
00082 QString fname, snum, line;
00083 QFile f;
00084 int nCount = 0;
00085 OrbitDataColl *ret;
00086
00087
00088
00089 n = n.lower();
00090
00091 if ((ret = dict[n])) {
00092
00093 return ret;
00094 }
00095
00096 ret = new OrbitDataColl;
00097
00098
00099 for (int i=0; i<6; ++i) {
00100 snum.setNum( i );
00101 fname = n + ".L" + snum + ".vsop";
00102 if (readOrbitData(fname, &(ret->Lon[i])))
00103 nCount++;
00104 }
00105
00106 if ( nCount==0 ){
00107 delete ret;
00108 return 0;
00109 }
00110
00111
00112 for (int i=0; i<6; ++i) {
00113 snum.setNum( i );
00114 fname = n + ".B" + snum + ".vsop";
00115 if (readOrbitData(fname, &(ret->Lat[i])))
00116 nCount++;
00117
00118 }
00119
00120 if (nCount==0){
00121 delete ret;
00122 return 0;
00123 }
00124
00125
00126 for (int i=0; i<6; ++i) {
00127 snum.setNum( i );
00128 fname = n + ".R" + snum + ".vsop";
00129 if (readOrbitData(fname, &(ret->Dst[i])))
00130 nCount++;
00131
00132 }
00133
00134 if (nCount==0){
00135 delete ret;
00136 return 0;
00137 }
00138
00139 dict.insert(n, ret);
00140
00141
00142
00143
00144 return ret;
00145 }
00146
00147 KSPlanet::KSPlanet( KStarsData *kd, QString s, QString imfile, double pSize )
00148 : KSPlanetBase(kd, s, imfile, pSize ), data_loaded(false) {
00149 }
00150
00151 bool KSPlanet::loadData() {
00152 return (odm.loadData(name()) != 0);
00153 }
00154
00155 void KSPlanet::calcEcliptic(double Tau, EclipticPosition &epret) const {
00156 double sum[6];
00157 OrbitDataColl * odc;
00158 double Tpow[6];
00159
00160 Tpow[0] = 1.0;
00161 for (int i=1; i<6; ++i) {
00162 Tpow[i] = Tpow[i-1] * Tau;
00163 }
00164
00165 if (!(odc = odm.loadData(name()))) {
00166 epret.longitude = 0.0;
00167 epret.latitude = 0.0;
00168 epret.radius = 0.0;
00169 kdError() << "Could not get data for '" << name() << "'" << endl;
00170 return;
00171 }
00172
00173
00174 for (int i=0; i<6; ++i) {
00175 sum[i] = 0.0;
00176 for (uint j = 0; j < odc->Lon[i].size(); ++j) {
00177 sum[i] += odc->Lon[i][j]->A * cos( odc->Lon[i][j]->B + odc->Lon[i][j]->C*Tau );
00178
00179
00180
00181
00182
00183 }
00184 sum[i] *= Tpow[i];
00185
00186 }
00187
00188 epret.longitude.setRadians( sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5] );
00189 epret.longitude.setD( epret.longitude.reduce().Degrees() );
00190
00191
00192 for (uint i=0; i<6; ++i) {
00193 sum[i] = 0.0;
00194 for (uint j = 0; j < odc->Lat[i].size(); ++j) {
00195 sum[i] += odc->Lat[i][j]->A * cos( odc->Lat[i][j]->B + odc->Lat[i][j]->C*Tau );
00196 }
00197 sum[i] *= Tpow[i];
00198 }
00199
00200
00201 epret.latitude.setRadians( sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5] );
00202
00203
00204 for (uint i=0; i<6; ++i) {
00205 sum[i] = 0.0;
00206 for (uint j = 0; j < odc->Dst[i].size(); ++j) {
00207 sum[i] += odc->Dst[i][j]->A * cos( odc->Dst[i][j]->B + odc->Dst[i][j]->C*Tau );
00208 }
00209 sum[i] *= Tpow[i];
00210 }
00211
00212 epret.radius = sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5];
00213
00214
00215
00216
00217
00218
00219 }
00220
00221 bool KSPlanet::findGeocentricPosition( const KSNumbers *num, const KSPlanetBase *Earth ) {
00222
00223 if ( Earth != NULL ) {
00224 double sinL, sinL0, sinB, sinB0;
00225 double cosL, cosL0, cosB, cosB0;
00226 double x = 0.0, y = 0.0, z = 0.0;
00227
00228 double olddst = -1000;
00229 double dst = 0;
00230
00231 EclipticPosition trialpos;
00232
00233 double jm = num->julianMillenia();
00234
00235 Earth->ecLong()->SinCos( sinL0, cosL0 );
00236 Earth->ecLat()->SinCos( sinB0, cosB0 );
00237
00238 double eX = Earth->rsun()*cosB0*cosL0;
00239 double eY = Earth->rsun()*cosB0*sinL0;
00240 double eZ = Earth->rsun()*sinB0;
00241
00242 bool once=true;
00243 while (fabs(dst - olddst) > .001) {
00244 calcEcliptic(jm, trialpos);
00245
00246
00247 if(once){
00248 helEcPos = trialpos;
00249 once=false;
00250 }
00251
00252 olddst = dst;
00253
00254 trialpos.longitude.SinCos( sinL, cosL );
00255 trialpos.latitude.SinCos( sinB, cosB );
00256
00257 x = trialpos.radius*cosB*cosL - eX;
00258 y = trialpos.radius*cosB*sinL - eY;
00259 z = trialpos.radius*sinB - eZ;
00260
00261
00262 dst = sqrt(x*x + y*y + z*z);
00263
00264 double delay = (.0057755183 * dst) / 365250.0;
00265
00266 jm = num->julianMillenia() - delay;
00267
00268 }
00269
00270 ep.longitude.setRadians( atan( y/x ) );
00271 if (x<0) ep.longitude.setD( ep.longitude.Degrees() + 180.0 );
00272 ep.latitude.setRadians( atan( z/( sqrt( x*x + y*y ) ) ) );
00273 setRsun( trialpos.radius );
00274 setRearth( dst );
00275
00276 EclipticToEquatorial( num->obliquity() );
00277
00278 nutate(num);
00279 aberrate(num);
00280
00281 } else {
00282
00283 calcEcliptic(num->julianMillenia(), ep);
00284 helEcPos = ep;
00285 }
00286
00287
00288 findPA( num );
00289
00290 return true;
00291 }