kstars

ksasteroid.cpp

00001 /***************************************************************************
00002                           ksasteroid.cpp  -  K Desktop Planetarium
00003                              -------------------
00004     begin                : Wed 19 Feb 2003
00005     copyright            : (C) 2001 by Jason Harris
00006     email                : jharris@30doradus.org
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
00015  *                                                                         *
00016  ***************************************************************************/
00017 
00018 #include <kdebug.h>
00019 
00020 #include "ksasteroid.h"
00021 #include "dms.h"
00022 #include "ksnumbers.h"
00023 #include "ksutils.h"
00024 #include "kstarsdata.h"
00025 
00026 KSAsteroid::KSAsteroid( KStarsData *_kd, QString s, QString imfile,
00027         long double _JD, double _a, double _e, dms _i, dms _w, dms _Node, dms _M, double _H )
00028  : KSPlanetBase(_kd, s, imfile), kd(_kd), JD(_JD), a(_a), e(_e), H(_H), i(_i), w(_w), M(_M), N(_Node) {
00029 
00030     setType( 10 ); //Asteroid
00031     setMag( H );
00032     //Compute the orbital Period from Kepler's 3rd law:
00033     P = 365.2568984 * pow(a, 1.5); //period in days
00034 }
00035 
00036 bool KSAsteroid::findGeocentricPosition( const KSNumbers *num, const KSPlanetBase *Earth ) {
00037     //Precess the longitude of the Ascending Node to the desired epoch:
00038     dms n = dms( double( N.Degrees() - 3.82394E-5 * ( num->julianDay() - J2000 )) ).reduce();
00039 
00040     //determine the mean anomaly for the desired date.  This is the mean anomaly for the
00041     //ephemeis epoch, plus the number of days between the desired date and ephemeris epoch,
00042     //times the asteroid's mean daily motion (360/P):
00043     dms m = dms( double( M.Degrees() + ( num->julianDay() - JD ) * 360.0/P ) ).reduce();
00044     double sinm, cosm;
00045     m.SinCos( sinm, cosm );
00046 
00047     //compute eccentric anomaly:
00048     double E = m.Degrees() + e*180.0/dms::PI * sinm * ( 1.0 + e*cosm );
00049 
00050     if ( e > 0.05 ) { //need more accurate approximation, iterate...
00051         double E0;
00052         int iter(0);
00053         do {
00054             E0 = E;
00055             iter++;
00056             E = E0 - ( E0 - e*180.0/dms::PI *sin( E0*dms::DegToRad ) - m.Degrees() )/(1 - e*cos( E0*dms::DegToRad ) );
00057         } while ( fabs( E - E0 ) > 0.001 && iter < 1000 );
00058     }
00059 
00060     double sinE, cosE;
00061     dms E1( E );
00062     E1.SinCos( sinE, cosE );
00063 
00064     double xv = a * ( cosE - e );
00065     double yv = a * sqrt( 1.0 - e*e ) * sinE;
00066 
00067     //v is the true anomaly; r is the distance from the Sun
00068 
00069     double v = atan( yv/xv ) / dms::DegToRad;
00070     //resolve atan ambiguity
00071     if ( xv < 0.0 ) v += 180.0;
00072 
00073     double r = sqrt( xv*xv + yv*yv );
00074 
00075     //vw is the sum of the true anomaly and the argument of perihelion
00076     dms vw( v + w.Degrees() );
00077     double sinN, cosN, sinvw, cosvw, sini, cosi;
00078 
00079     N.SinCos( sinN, cosN );
00080     vw.SinCos( sinvw, cosvw );
00081     i.SinCos( sini, cosi );
00082 
00083     //xh, yh, zh are the heliocentric cartesian coords with the ecliptic plane congruent with zh=0.
00084     double xh = r * ( cosN * cosvw - sinN * sinvw * cosi );
00085     double yh = r * ( sinN * cosvw + cosN * sinvw * cosi );
00086     double zh = r * ( sinvw * sini );
00087 
00088     //the spherical ecliptic coordinates:
00089     double ELongRad = atan( yh/xh );
00090     //resolve atan ambiguity
00091     if ( xh < 0.0 ) ELongRad += dms::PI;
00092     double ELatRad = atan( zh/r );  //(r can't possibly be negative, so no atan ambiguity)
00093 
00094     helEcPos.longitude.setRadians( ELongRad );
00095     helEcPos.latitude.setRadians( ELatRad );
00096     setRsun( r );
00097 
00098     if ( Earth ) {
00099         //xe, ye, ze are the Earth's heliocentric cartesian coords
00100         double cosBe, sinBe, cosLe, sinLe;
00101         Earth->ecLong()->SinCos( sinLe, cosLe );
00102         Earth->ecLat()->SinCos( sinBe, cosBe );
00103     
00104         double xe = Earth->rsun() * cosBe * cosLe;
00105         double ye = Earth->rsun() * cosBe * sinLe;
00106         double ze = Earth->rsun() * sinBe;
00107     
00108         //convert to geocentric ecliptic coordinates by subtracting Earth's coords:
00109         xh -= xe;
00110         yh -= ye;
00111         zh -= ze;
00112     }
00113 
00114     //the spherical geocentricecliptic coordinates:
00115     ELongRad = atan( yh/xh );
00116     //resolve atan ambiguity
00117     if ( xh < 0.0 ) ELongRad += dms::PI;
00118 
00119     double rr = sqrt( xh*xh + yh*yh + zh*zh );
00120     ELatRad = atan( zh/rr );  //(rr can't possibly be negative, so no atan ambiguity)
00121 
00122     ep.longitude.setRadians( ELongRad );
00123     ep.latitude.setRadians( ELatRad );
00124     if ( Earth ) setRearth( Earth );
00125     
00126     EclipticToEquatorial( num->obliquity() );
00127     nutate( num );
00128     aberrate( num );
00129 
00130     return true;
00131 }
00132 
00133 //Unused virtual function from KSPlanetBase
00134 bool KSAsteroid::loadData() { return false; }
00135
kstars

ksasteroid.cpp

kstars

API Reference
