kstars/kstars/ksplanet.cpp

   1 /***************************************************************************
   2                           ksplanet.cpp  -  K Desktop Planetarium
   3                              -------------------
   4     begin                : Sun Jul 22 2001
   5     copyright            : (C) 2001 by Jason Harris
   6     email                : jharris@30doradus.org
   7  ***************************************************************************/
   8
   9 /***************************************************************************
  10  *                                                                         *
  11  *   This program is free software; you can redistribute it and/or modify  *
  12  *   it under the terms of the GNU General Public License as published by  *
  13  *   the Free Software Foundation; either version 2 of the License, or     *
  14  *   (at your option) any later version.                                   *
  15  *                                                                         *
  16  ***************************************************************************/
  17
  18 #include <math.h>
  19 #include <kdebug.h>
  20 #include <qfile.h>
  21
  22 #include "ksplanet.h"
  23 #include "ksnumbers.h"
  24 #include "ksutils.h"
  25 #include "ksfilereader.h"
  26
  27 KSPlanet::OrbitDataManager KSPlanet::odm;
  28
  29 KSPlanet::OrbitDataColl::OrbitDataColl() {
  30 // avoid memory leaks
  31         for (int i=0; i<6; i++) {
  32                 Lon[i].setAutoDelete(true);
  33                 Lat[i].setAutoDelete(true);
  34                 Dst[i].setAutoDelete(true);
  35         }
  36 }
  37
  38
  39 KSPlanet::OrbitDataManager::OrbitDataManager() : dict(31, true) {
  40         // delete all data automatically to avoid a leak
  41         dict.setAutoDelete(true);
  42 }
  43
  44 bool KSPlanet::OrbitDataManager::readOrbitData(QString fname,
  45                 QPtrVector<KSPlanet::OrbitData> *vector) {
  46         QString line;
  47         QFile f;
  48         double A, B, C;
  49
  50         QPtrList<OrbitData> DData;
  51
  52         if ( KSUtils::openDataFile( f, fname ) ) {
  53                 KSFileReader fileReader( f ); // close file is included
  54     while ( fileReader.hasMoreLines() ) {
  55       line = fileReader.readLine();
  56                         QTextIStream instream( &line );
  57                         instream >> A >> B >> C;
  58                         DData.append(new OrbitData(A, B, C));
  59
  60                 }
  61 /* old code
  62                 QTextStream stream( &f );
  63                 while ( !stream.eof() ) {
  64                         line = stream.readLine();
  65                         QTextIStream instream( &line );
  66                         instream >> A >> B >> C;
  67                         DData.append(new OrbitData(A, B, C));
  68
  69                 }
  70                 f.close();
  71 */
  72         } else {
  73                 return false;
  74         }
  75
  76         DData.toVector(vector);
  77
  78         return true;
  79 }
  80
  81 KSPlanet::OrbitDataColl *KSPlanet::OrbitDataManager::loadData(QString n) {
  82         QString fname, snum, line;
  83         QFile f;
  84         int nCount = 0;
  85         OrbitDataColl *ret;
  86
  87 //      kdDebug() << k_funcinfo << " Loading data named " << n << endl;
  88
  89         n = n.lower();
  90
  91         if ((ret = dict[n])) {
  92 //              kdDebug() << k_funcinfo << " already loaded - returning" << endl;
  93                 return ret;
  94         }
  95
  96         ret = new OrbitDataColl;
  97
  98         //Ecliptic Longitude
  99         for (int i=0; i<6; ++i) {
 100                 snum.setNum( i );
 101                 fname = n + ".L" + snum + ".vsop";
 102                 if (readOrbitData(fname, &(ret->Lon[i])))
 103                         nCount++;
 104         }
 105
 106         if ( nCount==0 ){ //No longitude data found!
 107                 delete ret;
 108                 return 0;
 109         }
 110
 111         //Ecliptic Latitude
 112         for (int i=0; i<6; ++i) {
 113                 snum.setNum( i );
 114                 fname = n + ".B" + snum + ".vsop";
 115                 if (readOrbitData(fname, &(ret->Lat[i])))
 116                         nCount++;
 117
 118         }
 119
 120         if (nCount==0){ //no latitude data found!
 121                 delete ret;
 122                 return 0;
 123         }
 124
 125         //Heliocentric Distance
 126         for (int i=0; i<6; ++i) {
 127                 snum.setNum( i );
 128                 fname = n + ".R" + snum + ".vsop";
 129                 if (readOrbitData(fname, &(ret->Dst[i])))
 130                         nCount++;
 131
 132         }
 133
 134         if (nCount==0){ //no distance data found!
 135                 delete ret;
 136                 return 0;
 137         }
 138
 139         dict.insert(n, ret);
 140
 141
 142 //      kdDebug() << k_funcinfo << " successful load" << endl;
 143
 144         return ret;
 145 }
 146
 147 KSPlanet::KSPlanet( KStarsData *kd, QString s, QString imfile, double pSize )
 148  : KSPlanetBase(kd, s, imfile, pSize ), data_loaded(false) {
 149 }
 150
 151 bool KSPlanet::loadData() {
 152         return (odm.loadData(name()) != 0);
 153 }
 154
 155 void KSPlanet::calcEcliptic(double Tau, EclipticPosition &epret) const {
 156         double sum[6];
 157         OrbitDataColl * odc;
 158         double Tpow[6];
 159
 160         Tpow[0] = 1.0;
 161         for (int i=1; i<6; ++i) {
 162                 Tpow[i] = Tpow[i-1] * Tau;
 163         }
 164
 165         if (!(odc =  odm.loadData(name()))) {
 166                 epret.longitude = 0.0;
 167                 epret.latitude = 0.0;
 168                 epret.radius = 0.0;
 169                 kdError() << "Could not get data for '" << name() << "'" << endl;
 170                 return;
 171         }
 172
 173         //Ecliptic Longitude
 174         for (int i=0; i<6; ++i) {
 175                 sum[i] = 0.0;
 176                 for (uint j = 0; j < odc->Lon[i].size(); ++j) {
 177                         sum[i] += odc->Lon[i][j]->A * cos( odc->Lon[i][j]->B + odc->Lon[i][j]->C*Tau );
 178                         /*
 179                         kdDebug() << "sum[" << i <<"] =" << sum[i] <<
 180                                 " A = " << odc->Lon[i][j]->A << " B = " << odc->Lon[i][j]->B <<
 181                                 " C = " << odc->Lon[i][j]->C << endl;
 182                                 */
 183                 }
 184                 sum[i] *= Tpow[i];
 185                 //kdDebug() << name() << " : sum[" << i << "] = " << sum[i] <<endl;
 186         }
 187
 188         epret.longitude.setRadians( sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5] );
 189         epret.longitude.setD( epret.longitude.reduce().Degrees() );
 190
 191         //Compute Ecliptic Latitude
 192         for (uint i=0; i<6; ++i) {
 193                 sum[i] = 0.0;
 194                 for (uint j = 0; j < odc->Lat[i].size(); ++j) {
 195                         sum[i] += odc->Lat[i][j]->A * cos( odc->Lat[i][j]->B + odc->Lat[i][j]->C*Tau );
 196                 }
 197                 sum[i] *= Tpow[i];
 198         }
 199
 200
 201         epret.latitude.setRadians( sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5] );
 202
 203         //Compute Heliocentric Distance
 204         for (uint i=0; i<6; ++i) {
 205                 sum[i] = 0.0;
 206                 for (uint j = 0; j < odc->Dst[i].size(); ++j) {
 207                         sum[i] += odc->Dst[i][j]->A * cos( odc->Dst[i][j]->B + odc->Dst[i][j]->C*Tau );
 208                 }
 209                 sum[i] *= Tpow[i];
 210         }
 211
 212         epret.radius = sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5];
 213
 214         /*
 215         kdDebug() << name() << " pre: Lat = " << epret.latitude.toDMSString() << " Long = " <<
 216                 epret.longitude.toDMSString() << " Dist = " << epret.radius << endl;
 217         */
 218
 219 }
 220
 221 bool KSPlanet::findGeocentricPosition( const KSNumbers *num, const KSPlanetBase *Earth ) {
 222
 223         if ( Earth != NULL ) {
 224                 double sinL, sinL0, sinB, sinB0;
 225                 double cosL, cosL0, cosB, cosB0;
 226                 double x = 0.0, y = 0.0, z = 0.0;
 227
 228                 double olddst = -1000;
 229                 double dst = 0;
 230
 231                 EclipticPosition trialpos;
 232
 233                 double jm = num->julianMillenia();
 234
 235                 Earth->ecLong()->SinCos( sinL0, cosL0 );
 236                 Earth->ecLat()->SinCos( sinB0, cosB0 );
 237
 238                 double eX = Earth->rsun()*cosB0*cosL0;
 239                 double eY = Earth->rsun()*cosB0*sinL0;
 240                 double eZ = Earth->rsun()*sinB0;
 241
 242                 bool once=true;
 243                 while (fabs(dst - olddst) > .001) {
 244                         calcEcliptic(jm, trialpos);
 245
 246                         // We store the heliocentric ecliptic coordinates the first time they are computed.
 247                         if(once){
 248                                 helEcPos = trialpos;
 249                                 once=false;
 250                         }
 251
 252                         olddst = dst;
 253
 254                         trialpos.longitude.SinCos( sinL, cosL );
 255                         trialpos.latitude.SinCos( sinB, cosB );
 256
 257                         x = trialpos.radius*cosB*cosL - eX;
 258                         y = trialpos.radius*cosB*sinL - eY;
 259                         z = trialpos.radius*sinB - eZ;
 260
 261                         //distance from Earth
 262                         dst = sqrt(x*x + y*y + z*z);
 263
 264                         double delay = (.0057755183 * dst) / 365250.0;
 265
 266                         jm = num->julianMillenia() - delay;
 267
 268                 }
 269
 270                 ep.longitude.setRadians( atan( y/x ) );
 271                 if (x<0) ep.longitude.setD( ep.longitude.Degrees() + 180.0 ); //resolve atan ambiguity
 272                 ep.latitude.setRadians( atan( z/( sqrt( x*x + y*y ) ) ) );
 273                 setRsun( trialpos.radius );
 274                 setRearth( dst );
 275
 276                 EclipticToEquatorial( num->obliquity() );
 277
 278                 nutate(num);
 279                 aberrate(num);
 280
 281         } else {
 282
 283                 calcEcliptic(num->julianMillenia(), ep);
 284                 helEcPos = ep;
 285         }
 286
 287         //determine the position angle
 288         findPA( num );
 289
 290         return true;
 291 }