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 General Public License
19 * along with this program; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 **************************************************************************/
25 #include <simgear/debug/logstream.hxx>
27 #ifdef SG_MATH_EXCEPTION_CLASH
28 # define exception c_exception
32 #include "celestialBody.hxx"
36 /**************************************************************************
37 * void CelestialBody::updatePosition(double mjd, Star *ourSun)
39 * Basically, this member function provides a general interface for
40 * calculating the right ascension and declinaion. This function is
41 * used for calculating the planetary positions. For the planets, an
42 * overloaded member function is provided to additionally calculate the
44 * The sun and moon have their own overloaded updatePosition member, as their
45 * position is calculated an a slightly different manner.
48 * double mjd: provides the modified julian date.
49 * Star *ourSun: the sun's position is needed to convert heliocentric
50 * coordinates into geocentric coordinates.
54 *************************************************************************/
55 void CelestialBody::updatePosition(double mjd, Star *ourSun)
57 double eccAnom, v, ecl, actTime,
58 xv, yv, xh, yh, zh, xg, yg, zg, xe, ye, ze;
60 updateOrbElements(mjd);
61 actTime = sgCalcActTime(mjd);
63 // calcualate the angle bewteen ecliptic and equatorial coordinate system
64 ecl = SGD_DEGREES_TO_RADIANS * (23.4393 - 3.563E-7 *actTime);
66 eccAnom = sgCalcEccAnom(M, e); //calculate the eccentric anomaly
67 xv = a * (cos(eccAnom) - e);
68 yv = a * (sqrt (1.0 - e*e) * sin(eccAnom));
69 v = atan2(yv, xv); // the planet's true anomaly
70 r = sqrt (xv*xv + yv*yv); // the planet's distance
72 // calculate the planet's position in 3D space
73 xh = r * (cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(i));
74 yh = r * (sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(i));
75 zh = r * (sin(v+w) * sin(i));
77 // calculate the ecliptic longitude and latitude
78 xg = xh + ourSun->getxs();
79 yg = yh + ourSun->getys();
82 lonEcl = atan2(yh, xh);
83 latEcl = atan2(zh, sqrt(xh*xh+yh*yh));
86 ye = yg * cos(ecl) - zg * sin(ecl);
87 ze = yg * sin(ecl) + zg * cos(ecl);
88 rightAscension = atan2(ye, xe);
89 declination = atan2(ze, sqrt(xe*xe + ye*ye));
90 /* SG_LOG(SG_GENERAL, SG_INFO, "Planet found at : "
91 << rightAscension << " (ra), " << declination << " (dec)" ); */
93 //calculate some variables specific to calculating the magnitude
95 R = sqrt (xg*xg + yg*yg + zg*zg);
96 s = ourSun->getDistance();
98 // It is possible from these calculations for the argument to acos
99 // to exceed the valid range for acos(). So we do a little extra
102 double tmp = (r*r + R*R - s*s) / (2*r*R);
105 } else if ( tmp < -1.0) {
109 FV = SGD_RADIANS_TO_DEGREES * acos( tmp );
112 /****************************************************************************
113 * double CelestialBody::sgCalcEccAnom(double M, double e)
114 * this private member calculates the eccentric anomaly of a celestial body,
115 * given its mean anomaly and eccentricity.
117 * -Mean anomaly: the approximate angle between the perihelion and the current
118 * position. this angle increases uniformly with time.
120 * True anomaly: the actual angle between perihelion and current position.
122 * Eccentric anomaly: this is an auxilary angle, used in calculating the true
123 * anomaly from the mean anomaly.
125 * -eccentricity. Indicates the amount in which the orbit deviates from a
126 * circle (0 = circle, 0-1, is ellipse, 1 = parabola, > 1 = hyperbola).
128 * This function is also known as solveKeplersEquation()
131 * M: the mean anomaly
132 * e: the eccentricity
135 * the eccentric anomaly
137 ****************************************************************************/
138 double CelestialBody::sgCalcEccAnom(double M, double e)
141 eccAnom, E0, E1, diff;
143 eccAnom = M + e * sin(M) * (1.0 + e * cos (M));
144 // iterate to achieve a greater precision for larger eccentricities
150 E1 = E0 - (E0 - e * sin(E0) - M) / (1 - e *cos(E0));
151 diff = fabs(E0 - E1);
154 while (diff > (SGD_DEGREES_TO_RADIANS * 0.001));
160 /*****************************************************************************
161 * inline CelestialBody::CelestialBody
162 * public constructor for a generic celestialBody object.
163 * initializes the 6 primary orbital elements. The elements are:
164 * N: longitude of the ascending node
165 * i: inclination to the ecliptic
166 * w: argument of perihelion
167 * a: semi-major axis, or mean distance from the sun
170 * Each orbital element consists of a constant part and a variable part that
171 * gradually changes over time.
174 * the 13 arguments to the constructor constitute the first, constant
175 * ([NiwaeM]f) and the second variable ([NiwaeM]s) part of the orbital
176 * elements. The 13th argument is the current time. Note that the inclination
177 * is written with a capital (If, Is), because 'if' is a reserved word in the
178 * C/C++ programming language.
179 ***************************************************************************/
180 CelestialBody::CelestialBody(double Nf, double Ns,
181 double If, double Is,
182 double wf, double ws,
183 double af, double as,
184 double ef, double es,
185 double Mf, double Ms, double mjd)
187 NFirst = Nf; NSec = Ns;
188 iFirst = If; iSec = Is;
189 wFirst = wf; wSec = ws;
190 aFirst = af; aSec = as;
191 eFirst = ef; eSec = es;
192 MFirst = Mf; MSec = Ms;
193 updateOrbElements(mjd);
196 CelestialBody::CelestialBody(double Nf, double Ns,
197 double If, double Is,
198 double wf, double ws,
199 double af, double as,
200 double ef, double es,
201 double Mf, double Ms)
203 NFirst = Nf; NSec = Ns;
204 iFirst = If; iSec = Is;
205 wFirst = wf; wSec = ws;
206 aFirst = af; aSec = as;
207 eFirst = ef; eSec = es;
208 MFirst = Mf; MSec = Ms;
211 /****************************************************************************
212 * inline void CelestialBody::updateOrbElements(double mjd)
213 * given the current time, this private member calculates the actual
216 * Arguments: double mjd: the current modified julian date:
219 ***************************************************************************/
220 void CelestialBody::updateOrbElements(double mjd)
222 double actTime = sgCalcActTime(mjd);
223 M = SGD_DEGREES_TO_RADIANS * (MFirst + (MSec * actTime));
224 w = SGD_DEGREES_TO_RADIANS * (wFirst + (wSec * actTime));
225 N = SGD_DEGREES_TO_RADIANS * (NFirst + (NSec * actTime));
226 i = SGD_DEGREES_TO_RADIANS * (iFirst + (iSec * actTime));
227 e = eFirst + (eSec * actTime);
228 a = aFirst + (aSec * actTime);
231 /*****************************************************************************
232 * inline double CelestialBody::sgCalcActTime(double mjd)
233 * this private member function returns the offset in days from the epoch for
234 * wich the orbital elements are calculated (Jan, 1st, 2000).
236 * Argument: the current time
238 * return value: the (fractional) number of days until Jan 1, 2000.
239 ****************************************************************************/
240 double CelestialBody::sgCalcActTime(double mjd)
242 return (mjd - 36523.5);
245 /*****************************************************************************
246 * inline void CelestialBody::getPos(double* ra, double* dec)
247 * gives public access to Right Ascension and declination
249 ****************************************************************************/
250 void CelestialBody::getPos(double* ra, double* dec)
252 *ra = rightAscension;
256 /*****************************************************************************
257 * inline void CelestialBody::getPos(double* ra, double* dec, double* magnitude
258 * gives public acces to the current Right ascension, declination, and
260 ****************************************************************************/
261 void CelestialBody::getPos(double* ra, double* dec, double* magn)
263 *ra = rightAscension;