kstars/kstars/ksasteroid.cpp

   1 /***************************************************************************
   2                           ksasteroid.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 "ksasteroid.h"
  21 #include "dms.h"
  22 #include "ksnumbers.h"
  23 #include "ksutils.h"
  24 #include "kstarsdata.h"
  25
  26 KSAsteroid::KSAsteroid( KStarsData *_kd, QString s, QString imfile,
  27                 long double _JD, double _a, double _e, dms _i, dms _w, dms _Node, dms _M, double _H )
  28  : KSPlanetBase(_kd, s, imfile), kd(_kd), JD(_JD), a(_a), e(_e), H(_H), i(_i), w(_w), M(_M), N(_Node) {
  29
  30         setType( 10 ); //Asteroid
  31         setMag( H );
  32         //Compute the orbital Period from Kepler's 3rd law:
  33         P = 365.2568984 * pow(a, 1.5); //period in days
  34 }
  35
  36 bool KSAsteroid::findGeocentricPosition( const KSNumbers *num, const KSPlanetBase *Earth ) {
  37         //Precess the longitude of the Ascending Node to the desired epoch:
  38         dms n = dms( double( N.Degrees() - 3.82394E-5 * ( num->julianDay() - J2000 )) ).reduce();
  39
  40         //determine the mean anomaly for the desired date.  This is the mean anomaly for the
  41         //ephemeis epoch, plus the number of days between the desired date and ephemeris epoch,
  42         //times the asteroid's mean daily motion (360/P):
  43         dms m = dms( double( M.Degrees() + ( num->julianDay() - JD ) * 360.0/P ) ).reduce();
  44         double sinm, cosm;
  45         m.SinCos( sinm, cosm );
  46
  47         //compute eccentric anomaly:
  48         double E = m.Degrees() + e*180.0/dms::PI * sinm * ( 1.0 + e*cosm );
  49
  50         if ( e > 0.05 ) { //need more accurate approximation, iterate...
  51                 double E0;
  52                 int iter(0);
  53                 do {
  54                         E0 = E;
  55                         iter++;
  56                         E = E0 - ( E0 - e*180.0/dms::PI *sin( E0*dms::DegToRad ) - m.Degrees() )/(1 - e*cos( E0*dms::DegToRad ) );
  57                 } while ( fabs( E - E0 ) > 0.001 && iter < 1000 );
  58         }
  59
  60         double sinE, cosE;
  61         dms E1( E );
  62         E1.SinCos( sinE, cosE );
  63
  64         double xv = a * ( cosE - e );
  65         double yv = a * sqrt( 1.0 - e*e ) * sinE;
  66
  67         //v is the true anomaly; r is the distance from the Sun
  68
  69         double v = atan( yv/xv ) / dms::DegToRad;
  70         //resolve atan ambiguity
  71         if ( xv < 0.0 ) v += 180.0;
  72
  73         double r = sqrt( xv*xv + yv*yv );
  74
  75         //vw is the sum of the true anomaly and the argument of perihelion
  76         dms vw( v + w.Degrees() );
  77         double sinN, cosN, sinvw, cosvw, sini, cosi;
  78
  79         N.SinCos( sinN, cosN );
  80         vw.SinCos( sinvw, cosvw );
  81         i.SinCos( sini, cosi );
  82
  83         //xh, yh, zh are the heliocentric cartesian coords with the ecliptic plane congruent with zh=0.
  84         double xh = r * ( cosN * cosvw - sinN * sinvw * cosi );
  85         double yh = r * ( sinN * cosvw + cosN * sinvw * cosi );
  86         double zh = r * ( sinvw * sini );
  87
  88         //xe, ye, ze are the Earth's heliocentric cartesian coords
  89         double cosBe, sinBe, cosLe, sinLe;
  90         Earth->ecLong()->SinCos( sinLe, cosLe );
  91         Earth->ecLat()->SinCos( sinBe, cosBe );
  92
  93         double xe = Earth->rsun() * cosBe * cosLe;
  94         double ye = Earth->rsun() * cosBe * sinLe;
  95         double ze = Earth->rsun() * sinBe;
  96
  97         //convert to geocentric ecliptic coordinates by subtracting Earth's coords:
  98         xh -= xe;
  99         yh -= ye;
 100         zh -= ze;
 101
 102         //Finally, the spherical ecliptic coordinates:
 103         double ELongRad = atan( yh/xh );
 104         //resolve atan ambiguity
 105         if ( xh < 0.0 ) ELongRad += dms::PI;
 106
 107         double rr = sqrt( xh*xh + yh*yh );
 108         double ELatRad = atan( zh/rr );  //(rr can't possibly be negative, so no atan ambiguity)
 109
 110         ep.longitude.setRadians( ELongRad );
 111         ep.latitude.setRadians( ELatRad );
 112         setRsun( r );
 113         setRearth( Earth );
 114
 115         EclipticToEquatorial( num->obliquity() );
 116         nutate( num );
 117         aberrate( num );
 118
 119         return true;
 120 }
 121
 122 //Unused virtual function from KSPlanetBase
 123 bool KSAsteroid::loadData() { return false; }
 124