1 /**************************************************************************
3 * Written by Durk Talsma. Originally started October 1997, for distribution
4 * with the FlightGear project. Version 2 was written in August and
5 * September 1998. This code is based upon algorithms and data kindly
6 * provided by Mr. Paul Schlyter. (pausch@saaf.se).
8 * This library is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU Library General Public
10 * License as published by the Free Software Foundation; either
11 * version 2 of the License, or (at your option) any later version.
13 * This library is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * Library General Public License for more details.
18 * You should have received a copy of the GNU Library General Public
19 * License along with this library; if not, write to the
20 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
21 * Boston, MA 02111-1307, USA.
24 **************************************************************************/
26 #include <simgear/debug/logstream.hxx>
28 #ifdef SG_MATH_EXCEPTION_CLASH
29 # define exception c_exception
33 #include "celestialBody.hxx"
37 /**************************************************************************
38 * void CelestialBody::updatePosition(double mjd, Star *ourSun)
40 * Basically, this member function provides a general interface for
41 * calculating the right ascension and declinaion. This function is
42 * used for calculating the planetary positions. For the planets, an
43 * overloaded member function is provided to additionally calculate the
45 * The sun and moon have their own overloaded updatePosition member, as their
46 * position is calculated an a slightly different manner.
49 * double mjd: provides the modified julian date.
50 * Star *ourSun: the sun's position is needed to convert heliocentric
51 * coordinates into geocentric coordinates.
55 *************************************************************************/
56 void CelestialBody::updatePosition(double mjd, Star *ourSun)
58 double eccAnom, v, ecl, actTime,
59 xv, yv, xh, yh, zh, xg, yg, zg, xe, ye, ze;
61 updateOrbElements(mjd);
62 actTime = sgCalcActTime(mjd);
64 // calcualate the angle bewteen ecliptic and equatorial coordinate system
65 ecl = SGD_DEGREES_TO_RADIANS * (23.4393 - 3.563E-7 *actTime);
67 eccAnom = sgCalcEccAnom(M, e); //calculate the eccentric anomaly
68 xv = a * (cos(eccAnom) - e);
69 yv = a * (sqrt (1.0 - e*e) * sin(eccAnom));
70 v = atan2(yv, xv); // the planet's true anomaly
71 r = sqrt (xv*xv + yv*yv); // the planet's distance
73 // calculate the planet's position in 3D space
74 xh = r * (cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(i));
75 yh = r * (sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(i));
76 zh = r * (sin(v+w) * sin(i));
78 // calculate the ecliptic longitude and latitude
79 xg = xh + ourSun->getxs();
80 yg = yh + ourSun->getys();
83 lonEcl = atan2(yh, xh);
84 latEcl = atan2(zh, sqrt(xh*xh+yh*yh));
87 ye = yg * cos(ecl) - zg * sin(ecl);
88 ze = yg * sin(ecl) + zg * cos(ecl);
89 rightAscension = atan2(ye, xe);
90 declination = atan2(ze, sqrt(xe*xe + ye*ye));
91 /* SG_LOG(SG_GENERAL, SG_INFO, "Planet found at : "
92 << rightAscension << " (ra), " << declination << " (dec)" ); */
94 //calculate some variables specific to calculating the magnitude
96 R = sqrt (xg*xg + yg*yg + zg*zg);
97 s = ourSun->getDistance();
99 // It is possible from these calculations for the argument to acos
100 // to exceed the valid range for acos(). So we do a little extra
103 double tmp = (r*r + R*R - s*s) / (2*r*R);
106 } else if ( tmp < -1.0) {
110 FV = SGD_RADIANS_TO_DEGREES * acos( tmp );
113 /****************************************************************************
114 * double CelestialBody::sgCalcEccAnom(double M, double e)
115 * this private member calculates the eccentric anomaly of a celestial body,
116 * given its mean anomaly and eccentricity.
118 * -Mean anomaly: the approximate angle between the perihelion and the current
119 * position. this angle increases uniformly with time.
121 * True anomaly: the actual angle between perihelion and current position.
123 * Eccentric anomaly: this is an auxilary angle, used in calculating the true
124 * anomaly from the mean anomaly.
126 * -eccentricity. Indicates the amount in which the orbit deviates from a
127 * circle (0 = circle, 0-1, is ellipse, 1 = parabola, > 1 = hyperbola).
129 * This function is also known as solveKeplersEquation()
132 * M: the mean anomaly
133 * e: the eccentricity
136 * the eccentric anomaly
138 ****************************************************************************/
139 double CelestialBody::sgCalcEccAnom(double M, double e)
142 eccAnom, E0, E1, diff;
144 eccAnom = M + e * sin(M) * (1.0 + e * cos (M));
145 // iterate to achieve a greater precision for larger eccentricities
151 E1 = E0 - (E0 - e * sin(E0) - M) / (1 - e *cos(E0));
152 diff = fabs(E0 - E1);
155 while (diff > (SGD_DEGREES_TO_RADIANS * 0.001));
161 /*****************************************************************************
162 * inline CelestialBody::CelestialBody
163 * public constructor for a generic celestialBody object.
164 * initializes the 6 primary orbital elements. The elements are:
165 * N: longitude of the ascending node
166 * i: inclination to the ecliptic
167 * w: argument of perihelion
168 * a: semi-major axis, or mean distance from the sun
171 * Each orbital element consists of a constant part and a variable part that
172 * gradually changes over time.
175 * the 13 arguments to the constructor constitute the first, constant
176 * ([NiwaeM]f) and the second variable ([NiwaeM]s) part of the orbital
177 * elements. The 13th argument is the current time. Note that the inclination
178 * is written with a capital (If, Is), because 'if' is a reserved word in the
179 * C/C++ programming language.
180 ***************************************************************************/
181 CelestialBody::CelestialBody(double Nf, double Ns,
182 double If, double Is,
183 double wf, double ws,
184 double af, double as,
185 double ef, double es,
186 double Mf, double Ms, double mjd)
188 NFirst = Nf; NSec = Ns;
189 iFirst = If; iSec = Is;
190 wFirst = wf; wSec = ws;
191 aFirst = af; aSec = as;
192 eFirst = ef; eSec = es;
193 MFirst = Mf; MSec = Ms;
194 updateOrbElements(mjd);
197 CelestialBody::CelestialBody(double Nf, double Ns,
198 double If, double Is,
199 double wf, double ws,
200 double af, double as,
201 double ef, double es,
202 double Mf, double Ms)
204 NFirst = Nf; NSec = Ns;
205 iFirst = If; iSec = Is;
206 wFirst = wf; wSec = ws;
207 aFirst = af; aSec = as;
208 eFirst = ef; eSec = es;
209 MFirst = Mf; MSec = Ms;
212 /****************************************************************************
213 * inline void CelestialBody::updateOrbElements(double mjd)
214 * given the current time, this private member calculates the actual
217 * Arguments: double mjd: the current modified julian date:
220 ***************************************************************************/
221 void CelestialBody::updateOrbElements(double mjd)
223 double actTime = sgCalcActTime(mjd);
224 M = SGD_DEGREES_TO_RADIANS * (MFirst + (MSec * actTime));
225 w = SGD_DEGREES_TO_RADIANS * (wFirst + (wSec * actTime));
226 N = SGD_DEGREES_TO_RADIANS * (NFirst + (NSec * actTime));
227 i = SGD_DEGREES_TO_RADIANS * (iFirst + (iSec * actTime));
228 e = eFirst + (eSec * actTime);
229 a = aFirst + (aSec * actTime);
232 /*****************************************************************************
233 * inline double CelestialBody::sgCalcActTime(double mjd)
234 * this private member function returns the offset in days from the epoch for
235 * wich the orbital elements are calculated (Jan, 1st, 2000).
237 * Argument: the current time
239 * return value: the (fractional) number of days until Jan 1, 2000.
240 ****************************************************************************/
241 double CelestialBody::sgCalcActTime(double mjd)
243 return (mjd - 36523.5);
246 /*****************************************************************************
247 * inline void CelestialBody::getPos(double* ra, double* dec)
248 * gives public access to Right Ascension and declination
250 ****************************************************************************/
251 void CelestialBody::getPos(double* ra, double* dec)
253 *ra = rightAscension;
257 /*****************************************************************************
258 * inline void CelestialBody::getPos(double* ra, double* dec, double* magnitude
259 * gives public acces to the current Right ascension, declination, and
261 ****************************************************************************/
262 void CelestialBody::getPos(double* ra, double* dec, double* magn)
264 *ra = rightAscension;