1 /***************************************************************************
2 ksasteroid.cpp - K Desktop Planetarium
4 begin : Wed 19 Feb 2003
5 copyright : (C) 2001 by Jason Harris
6 email : jharris@30doradus.org
7 ***************************************************************************/
9 /***************************************************************************
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. *
16 ***************************************************************************/
20 #include "ksasteroid.h"
22 #include "ksnumbers.h"
24 #include "kstarsdata.h"
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
) {
30 setType( 10 ); //Asteroid
32 //Compute the orbital Period from Kepler's 3rd law:
33 P
= 365.2568984 * pow(a
, 1.5); //period in days
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();
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();
45 m
.SinCos( sinm
, cosm
);
47 //compute eccentric anomaly:
48 double E
= m
.Degrees() + e
*180.0/dms::PI
* sinm
* ( 1.0 + e
*cosm
);
50 if ( e
> 0.05 ) { //need more accurate approximation, iterate...
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 );
62 E1
.SinCos( sinE
, cosE
);
64 double xv
= a
* ( cosE
- e
);
65 double yv
= a
* sqrt( 1.0 - e
*e
) * sinE
;
67 //v is the true anomaly; r is the distance from the Sun
69 double v
= atan( yv
/xv
) / dms::DegToRad
;
70 //resolve atan ambiguity
71 if ( xv
< 0.0 ) v
+= 180.0;
73 double r
= sqrt( xv
*xv
+ yv
*yv
);
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
;
79 N
.SinCos( sinN
, cosN
);
80 vw
.SinCos( sinvw
, cosvw
);
81 i
.SinCos( sini
, cosi
);
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
);
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
);
93 double xe
= Earth
->rsun() * cosBe
* cosLe
;
94 double ye
= Earth
->rsun() * cosBe
* sinLe
;
95 double ze
= Earth
->rsun() * sinBe
;
97 //convert to geocentric ecliptic coordinates by subtracting Earth's coords:
102 //Finally, the spherical ecliptic coordinates:
103 double ELongRad
= atan( yh
/xh
);
104 //resolve atan ambiguity
105 if ( xh
< 0.0 ) ELongRad
+= dms::PI
;
107 double rr
= sqrt( xh
*xh
+ yh
*yh
);
108 double ELatRad
= atan( zh
/rr
); //(rr can't possibly be negative, so no atan ambiguity)
110 ep
.longitude
.setRadians( ELongRad
);
111 ep
.latitude
.setRadians( ELatRad
);
115 EclipticToEquatorial( num
->obliquity() );
122 //Unused virtual function from KSPlanetBase
123 bool KSAsteroid::loadData() { return false; }