kstars/kstars/kscomet.cpp

   1 /***************************************************************************
   2                           kscomet.cpp  -  K Desktop Planetarium
   3                              -------------------
   4     begin                : Wed 19 Feb 2003
   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 <kdebug.h>
  19
  20 #include "kstarsdata.h"
  21 #include "kstarsdatetime.h"
  22 #include "ksnumbers.h"
  23 #include "dms.h"
  24 #include "kscomet.h"
  25
  26
  27 KSComet::KSComet( KStarsData *_kd, QString _s, QString imfile,
  28                 long double _JD, double _q, double _e, dms _i, dms _w, dms _Node, double Tp )
  29  : KSPlanetBase(_kd, _s, imfile), kd(_kd), JD(_JD), q(_q), e(_e), i(_i), w(_w), N(_Node) {
  30
  31         setType( 9 ); //Comet
  32
  33         //Find the Julian Day of Perihelion from Tp
  34         //Tp is a double which encodes a date like: YYYYMMDD.DDDDD (e.g., 19730521.33333
  35         int year = int( Tp/10000.0 );
  36         int month = int( (int(Tp) % 10000)/100.0 );
  37         int day = int( int(Tp) % 100 );
  38         double Hour = 24.0 * ( Tp - int(Tp) );
  39         int h = int( Hour );
  40         int m = int( 60.0 * ( Hour - h ) );
  41         int s = int( 60.0 * ( 60.0 * ( Hour - h) - m ) );
  42
  43         JDp = KStarsDateTime( ExtDate( year, month, day ), QTime( h, m, s ) ).djd();
  44
  45         //compute the semi-major axis, a:
  46         a = q/(1.0-e);
  47
  48         //Compute the orbital Period from Kepler's 3rd law:
  49         P = 365.2568984 * pow(a, 1.5); //period in days
  50
  51         //If the name contains a "/", make this name2 and make name a truncated version without the leading "P/" or "C/"
  52         if ( name().contains( "/" ) ) {
  53                 setLongName( name() );
  54                 setName( name().mid( name().find("/") + 1 ) );
  55         }
  56 }
  57
  58 bool KSComet::findGeocentricPosition( const KSNumbers *num, const KSPlanetBase *Earth ) {
  59         double v(0.0), r(0.0);
  60
  61         //Precess the longitude of the Ascending Node to the desired epoch:
  62         dms n = dms( double(N.Degrees() - 3.82394E-5 * ( num->julianDay() - J2000 )) ).reduce();
  63
  64         if ( e > 0.98 ) {
  65                 //Use near-parabolic approximation
  66                 double k = 0.01720209895; //Gauss gravitational constant
  67                 double a = 0.75 * ( num->julianDay() - JDp ) * k * sqrt( (1+e)/(q*q*q) );
  68                 double b = sqrt( 1.0 + a*a );
  69                 double W = pow((b+a),1.0/3.0) - pow((b-a),1.0/3.0);
  70                 double c = 1.0 + 1.0/(W*W);
  71                 double f = (1.0-e)/(1.0+e);
  72                 double g = f/(c*c);
  73
  74                 double a1 = (2.0/3.0) + (2.0*W*W/5.0);
  75                 double a2 = (7.0/5.0) + (33.0*W*W/35.0) + (37.0*W*W*W*W/175.0);
  76                 double a3 = W*W*( (432.0/175.0) + (956.0*W*W/1125.0) + (84.0*W*W*W*W/1575.0) );
  77                 double w = W*(1.0 + g*c*( a1 + a2*g + a3*g*g ));
  78
  79                 v = 2.0*atan(w) / dms::DegToRad;
  80                 r = q*( 1.0 + w*w )/( 1.0 + w*w*f );
  81         } else {
  82                 //Use normal ellipse method
  83                 //Determine Mean anomaly for desired date:
  84                 dms m = dms( double(360.0*( num->julianDay() - JDp )/P) ).reduce();
  85                 double sinm, cosm;
  86                 m.SinCos( sinm, cosm );
  87
  88                 //compute eccentric anomaly:
  89                 double E = m.Degrees() + e*180.0/dms::PI * sinm * ( 1.0 + e*cosm );
  90
  91                 if ( e > 0.05 ) { //need more accurate approximation, iterate...
  92                         double E0;
  93                         int iter(0);
  94                         do {
  95                                 E0 = E;
  96                                 iter++;
  97                                 E = E0 - ( E0 - e*180.0/dms::PI *sin( E0*dms::DegToRad ) - m.Degrees() )/(1 - e*cos( E0*dms::DegToRad ) );
  98                         } while ( fabs( E - E0 ) > 0.001 && iter < 1000 );
  99                 }
 100
 101                 double sinE, cosE;
 102                 dms E1( E );
 103                 E1.SinCos( sinE, cosE );
 104
 105                 double xv = a * ( cosE - e );
 106                 double yv = a * sqrt( 1.0 - e*e ) * sinE;
 107
 108                 //v is the true anomaly; r is the distance from the Sun
 109
 110                 v = atan( yv/xv ) / dms::DegToRad;
 111                 //resolve atan ambiguity
 112                 if ( xv < 0.0 ) v += 180.0;
 113
 114                 r = sqrt( xv*xv + yv*yv );
 115         }
 116
 117         //vw is the sum of the true anomaly and the argument of perihelion
 118         dms vw( v + w.Degrees() );
 119         double sinN, cosN, sinvw, cosvw, sini, cosi;
 120
 121         n.SinCos( sinN, cosN );
 122         vw.SinCos( sinvw, cosvw );
 123         i.SinCos( sini, cosi );
 124
 125         //xh, yh, zh are the heliocentric cartesian coords with the ecliptic plane congruent with zh=0.
 126         double xh = r * ( cosN * cosvw - sinN * sinvw * cosi );
 127         double yh = r * ( sinN * cosvw + cosN * sinvw * cosi );
 128         double zh = r * ( sinvw * sini );
 129
 130         //xe, ye, ze are the Earth's heliocentric cartesian coords
 131         double cosBe, sinBe, cosLe, sinLe;
 132         Earth->ecLong()->SinCos( sinLe, cosLe );
 133         Earth->ecLat()->SinCos( sinBe, cosBe );
 134
 135         double xe = Earth->rsun() * cosBe * cosLe;
 136         double ye = Earth->rsun() * cosBe * sinLe;
 137         double ze = Earth->rsun() * sinBe;
 138
 139         //convert to geocentric ecliptic coordinates by subtracting Earth's coords:
 140         xh -= xe;
 141         yh -= ye;
 142         zh -= ze;
 143
 144         //Finally, the spherical ecliptic coordinates:
 145         double ELongRad = atan( yh/xh );
 146         //resolve atan ambiguity
 147         if ( xh < 0.0 ) ELongRad += dms::PI;
 148
 149         double rr = sqrt( xh*xh + yh*yh );
 150         double ELatRad = atan( zh/rr );  //(rr can't possibly be negative, so no atan ambiguity)
 151
 152         ep.longitude.setRadians( ELongRad );
 153         ep.latitude.setRadians( ELatRad );
 154         setRsun( r );
 155         setRearth( Earth );
 156
 157         EclipticToEquatorial( num->obliquity() );
 158         nutate( num );
 159         aberrate( num );
 160
 161         return true;
 162 }
 163
 164 //Unused virtual function from KSPlanetBase
 165 bool KSComet::loadData() { return false; }