kstars/kstars/kspluto.cpp

   1 /***************************************************************************
   2                           kspluto.cpp  -  K Desktop Planetarium
   3                              -------------------
   4     begin                : Mon Sep 24 2001
   5     copyright            : (C) 2001 by Jason Harris
   6     email                : kstars@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
  20 #include <qfile.h>
  21 #include <kdebug.h>
  22
  23 #include "kspluto.h"
  24 #include "ksutils.h"
  25 #include "ksnumbers.h"
  26
  27 #ifdef B0
  28 // There are systems that #define B0 as a special termios flag (for baud rate)
  29 #undef B0
  30 #endif
  31
  32 int KSPluto::DATAARRAYSIZE = 106;
  33 bool KSPluto::data_loaded = false;
  34 double *KSPluto::freq = 0;
  35 KSPluto::XYZData *KSPluto::xdata = 0;
  36 KSPluto::XYZData *KSPluto::ydata = 0;
  37 KSPluto::XYZData *KSPluto::zdata = 0;
  38 int KSPluto::objects = 0;
  39
  40 KSPluto::KSPluto(KStarsData *kd, QString fn, double pSize ) : KSPlanetBase( kd, I18N_NOOP( "Pluto" ), fn, pSize ) {
  41         objects++;
  42 }
  43
  44 KSPluto::~KSPluto() {
  45         objects--;
  46         // delete arrays if all other objects are closed
  47         if (!objects) {
  48                 if (freq) {
  49                         delete [] freq;
  50                         freq = 0;
  51                 }
  52                 if (xdata) {
  53                         delete [] xdata;
  54                         xdata = 0;
  55                 }
  56                 if (ydata) {
  57                         delete [] ydata;
  58                         ydata = 0;
  59                 }
  60                 if (zdata) {
  61                         delete [] zdata;
  62                         zdata = 0;
  63                 }
  64         }
  65 }
  66
  67 bool KSPluto::loadData() {
  68         return loadData("pluto");
  69 }
  70
  71 bool KSPluto::loadData(QString fn) {
  72         QString fname, snum, line;
  73         QFile f;
  74
  75         if (data_loaded) return true;
  76
  77 //      kdDebug() << "inside pluto loadData" << endl;
  78         freq = new double[DATAARRAYSIZE];
  79         xdata = new XYZData[DATAARRAYSIZE];
  80         ydata = new XYZData[DATAARRAYSIZE];
  81         zdata = new XYZData[DATAARRAYSIZE];
  82
  83         //read in the periodic frequencies
  84         int n = 0;
  85         if ( KSUtils::openDataFile( f, fn.lower() + ".freq" ) ) {
  86                 QTextStream stream( &f );
  87                 while ( !stream.eof() ) {
  88                         line = stream.readLine();
  89                         QTextIStream instream( &line );
  90                         instream >> freq[n++];
  91
  92                         if (n > DATAARRAYSIZE) {
  93                                 kdDebug() << "We tried to read more than " << DATAARRAYSIZE <<
  94                                         "elements into KSPluto::freq" << endl;
  95                                 return false;
  96                         }
  97                 }
  98                 f.close();
  99         }
 100
 101         if (n==0) return false;
 102
 103         n=0;
 104         if ( KSUtils::openDataFile( f, fn.lower() + ".x" ) ) {
 105                 QTextStream stream( &f );
 106                 while ( !stream.eof() ) {
 107                         double ac,as;
 108                         line = stream.readLine();
 109                         QTextIStream instream( &line );
 110                         instream >> ac >> as;
 111                         xdata[n++] = XYZData(ac, as);
 112
 113                         if (n > DATAARRAYSIZE) {
 114                                 kdDebug() << "We tried to read more than " << DATAARRAYSIZE <<
 115                                         "elements into KSPluto::xdata" << endl;
 116                                 return false;
 117                         }
 118                 }
 119                 f.close();
 120         }
 121
 122         if (n==0) return false;
 123
 124         n=0;
 125         if ( KSUtils::openDataFile( f, fn.lower() + ".y" ) ) {
 126                 QTextStream stream( &f );
 127                 while ( !stream.eof() ) {
 128                         double ac,as;
 129                         line = stream.readLine();
 130                         QTextIStream instream( &line );
 131                         instream >> ac >> as;
 132                         ydata[n++] = XYZData(ac, as);
 133
 134                         if (n > DATAARRAYSIZE) {
 135                                 kdDebug() << "We tried to read more than " << DATAARRAYSIZE <<
 136                                         "elements into KSPluto::ydata" << endl;
 137                                 return false;
 138                         }
 139                 }
 140                 f.close();
 141         }
 142
 143         if (n==0) return false;
 144
 145         n=0;
 146         if ( KSUtils::openDataFile( f, fn.lower() + ".z" ) ) {
 147                 QTextStream stream( &f );
 148                 while ( !stream.eof() ) {
 149                         double ac,as;
 150                         line = stream.readLine();
 151                         QTextIStream instream( &line );
 152                         instream >> ac >> as;
 153                         zdata[n++] = XYZData(ac, as);
 154
 155                         if (n > DATAARRAYSIZE) {
 156                                 kdDebug() << "We tried to read more than " << DATAARRAYSIZE <<
 157                                         "elements into KSPluto::zdata" << endl;
 158                                 return false;
 159                         }
 160                 }
 161                 f.close();
 162         }
 163
 164         if (n==0) return false;
 165
 166         data_loaded=true;
 167         return true;
 168 }
 169
 170
 171 KSPluto::XYZpos KSPluto::calcRectCoords(double jd)  {
 172         double X, Y, Z, U;
 173
 174         loadData(name());
 175
 176         double T = 2.0*(jd -  2341972.5)/146120.0 - 1.0;
 177
 178         U =  T*73060.0;
 179
 180         //compute secular terms
 181         X = 98083308510. - 1465718392.*T + 11528487809.*T*T + 55397965917.*T*T*T;
 182         Y = 101846243715. + 57789.*T - 5487929294.*T*T + 8520205290.*T*T*T;
 183         Z = 2183700004. + 433209785.*T - 4911803413.*T*T - 14029741184.*T*T*T;
 184
 185
 186 // compute periodic X terms
 187         for(int n=0; n < DATAARRAYSIZE; ++n) {
 188                 double tempdouble = xdata[n].ac*cos( freq[n]*U ) + xdata[n].as*sin( freq[n]*U );
 189                 if (n > 81) tempdouble *= T;
 190                 if (n > 100) tempdouble *= T;
 191
 192                 X += tempdouble;
 193         }
 194
 195 // compute periodic Y terms
 196         for(int n=0; n < DATAARRAYSIZE; ++n) {
 197                 double tempdouble = ydata[n].ac*cos( freq[n]*U ) + ydata[n].as*sin( freq[n]*U );
 198                 if (n > 81) tempdouble *= T;
 199                 if (n > 100) tempdouble *= T;
 200
 201                 Y += tempdouble;
 202         }
 203
 204 // compute periodic Z terms
 205         for(int n=0; n < DATAARRAYSIZE; ++n) {
 206                 double tempdouble = zdata[n].ac*cos( freq[n]*U ) + zdata[n].as*sin( freq[n]*U );
 207                 if (n > 81) tempdouble *= T;
 208                 if (n > 100) tempdouble *= T;
 209
 210                 Z += tempdouble;
 211         }
 212
 213         return XYZpos(X,Y,Z);
 214 }
 215
 216 bool KSPluto::findGeocentricPosition( const KSNumbers *num, const KSPlanetBase *Earth ) {
 217         double X0, Y0, Z0, RARad;
 218         dms L0, B0; //geocentric ecliptic coords of Sun
 219         dms EarthLong, EarthLat; //heliocentric ecliptic coords of Earth
 220         double sinL0, cosL0, sinB0, cosB0;
 221
 222         XYZpos pos;
 223
 224         double olddst = -1000;
 225         double dst = 0;
 226
 227         double jd = num->julianDay();
 228
 229         Earth->ecLong()->SinCos( sinL0, cosL0 );
 230         Earth->ecLat()->SinCos( sinB0, cosB0 );
 231
 232         double eX = Earth->rsun()*cosB0*cosL0;
 233         double eY = Earth->rsun()*cosB0*sinL0;
 234         double eZ = Earth->rsun()*sinB0;
 235
 236         while (fabs(dst-olddst) > .001) {
 237                 pos = calcRectCoords(jd);
 238                 //convert X, Y, Z to AU (given in 1.0E10 AU)
 239                 pos.X /= 10000000000.0; pos.Y /= 10000000000.0; pos.Z /= 10000000000.0;
 240
 241                 /*
 242                 kdDebug() << "pluto: X = " << pos.X << " Y = " << pos.Y << " Z = " << pos.Z << endl;
 243                 */
 244
 245                 double dX = pos.X - eX;
 246                 double dY = pos.Y - eY;
 247                 double dZ = pos.Z - eZ;
 248
 249                 olddst = dst;
 250                 dst = sqrt( dX * dX + dY * dY + dZ * dZ );
 251                 double delay = .0057755183 * dst;
 252
 253                 jd = num->julianDay() - delay;
 254
 255         }
 256
 257         //Compute the heliocentric ecliptic coords
 258         setRsun( sqrt( pos.X*pos.X + pos.Y*pos.Y + pos.Z*pos.Z ) );
 259         L0.setRadians( atan( pos.Y / pos.X ) );
 260         if ( pos.X < 0.0 ) L0.setD( L0.Degrees() + 180.0 );
 261
 262 //      if ( pos.X > 0.0 && pos.Y < 0.0 ) L += 360.0;
 263         setHelEcLong( L0 );
 264
 265         B0.setRadians( asin( pos.Z / rsun() ) );
 266         setHelEcLat( B0 );
 267
 268         //L0, B0 are Sun's Ecliptic coords (L0 = Learth + 180; B0 = -Bearth)
 269         L0.setD( Earth->ecLong()->Degrees() + 180.0 );
 270         L0.setD( L0.reduce().Degrees() );
 271         B0.setD( -1.0*Earth->ecLat()->Degrees() );
 272
 273         L0.SinCos( sinL0, cosL0 );
 274         B0.SinCos( sinB0, cosB0 );
 275
 276         double cosOb, sinOb;
 277         num->obliquity()->SinCos( sinOb, cosOb );
 278
 279         X0 = Earth->rsun()*cosB0*cosL0;
 280         Y0 = Earth->rsun()*( cosB0*sinL0*cosOb - sinB0*sinOb );
 281         Z0 = Earth->rsun()*( cosB0*sinL0*sinOb + sinB0*cosOb );
 282
 283         //transform to geocentric rectangular coordinates by adding Sun's values
 284         pos.X += X0; pos.Y += Y0; pos.Z += Z0;
 285
 286         //Use Meeus's Eq. 32.10 to find Rsun, RA and Dec:
 287         RARad = atan( pos.Y / pos.X );
 288         if ( pos.X<0 ) RARad += dms::PI;
 289         if ( pos.X>0 && pos.Y<0 ) RARad += 2.0*dms::PI;
 290         dms newRA; newRA.setRadians( RARad );
 291         dms newDec; newDec.setRadians( asin( pos.Z/rsun() ) );
 292         setRA( newRA );
 293         setDec( newDec );
 294
 295         //compute Ecliptic coordinates
 296         EquatorialToEcliptic( num->obliquity() );
 297
 298         //determine the position angle
 299         findPA( num );
 300
 301         setRearth( Earth );
 302
 303         return true;
 304 }