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>
29 #include "celestialBody.hxx"
33 /**************************************************************************
34 * void CelestialBody::updatePosition(double mjd, Star *ourSun)
36 * Basically, this member function provides a general interface for
37 * calculating the right ascension and declinaion. This function is
38 * used for calculating the planetary positions. For the planets, an
39 * overloaded member function is provided to additionally calculate the
41 * The sun and moon have their own overloaded updatePosition member, as their
42 * position is calculated an a slightly different manner.
45 * double mjd: provides the modified julian date.
46 * Star *ourSun: the sun's position is needed to convert heliocentric
47 * coordinates into geocentric coordinates.
51 *************************************************************************/
52 void CelestialBody::updatePosition(double mjd, Star *ourSun)
54 double eccAnom, v, ecl, actTime,
55 xv, yv, xh, yh, zh, xg, yg, zg, xe, ye, ze;
57 updateOrbElements(mjd);
58 actTime = sgCalcActTime(mjd);
60 // calcualate the angle bewteen ecliptic and equatorial coordinate system
61 ecl = SGD_DEGREES_TO_RADIANS * (23.4393 - 3.563E-7 *actTime);
63 eccAnom = sgCalcEccAnom(M, e); //calculate the eccentric anomaly
64 xv = a * (cos(eccAnom) - e);
65 yv = a * (sqrt (1.0 - e*e) * sin(eccAnom));
66 v = atan2(yv, xv); // the planet's true anomaly
67 r = sqrt (xv*xv + yv*yv); // the planet's distance
69 // calculate the planet's position in 3D space
70 xh = r * (cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(i));
71 yh = r * (sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(i));
72 zh = r * (sin(v+w) * sin(i));
74 // calculate the ecliptic longitude and latitude
75 xg = xh + ourSun->getxs();
76 yg = yh + ourSun->getys();
79 lonEcl = atan2(yh, xh);
80 latEcl = atan2(zh, sqrt(xh*xh+yh*yh));
83 ye = yg * cos(ecl) - zg * sin(ecl);
84 ze = yg * sin(ecl) + zg * cos(ecl);
85 rightAscension = atan2(ye, xe);
86 declination = atan2(ze, sqrt(xe*xe + ye*ye));
87 /* SG_LOG(SG_GENERAL, SG_INFO, "Planet found at : "
88 << rightAscension << " (ra), " << declination << " (dec)" ); */
90 //calculate some variables specific to calculating the magnitude
92 R = sqrt (xg*xg + yg*yg + zg*zg);
93 s = ourSun->getDistance();
95 // It is possible from these calculations for the argument to acos
96 // to exceed the valid range for acos(). So we do a little extra
99 double tmp = (r*r + R*R - s*s) / (2*r*R);
102 } else if ( tmp < -1.0) {
106 FV = SGD_RADIANS_TO_DEGREES * acos( tmp );
109 /****************************************************************************
110 * double CelestialBody::sgCalcEccAnom(double M, double e)
111 * this private member calculates the eccentric anomaly of a celestial body,
112 * given its mean anomaly and eccentricity.
114 * -Mean anomaly: the approximate angle between the perihelion and the current
115 * position. this angle increases uniformly with time.
117 * True anomaly: the actual angle between perihelion and current position.
119 * Eccentric anomaly: this is an auxilary angle, used in calculating the true
120 * anomaly from the mean anomaly.
122 * -eccentricity. Indicates the amount in which the orbit deviates from a
123 * circle (0 = circle, 0-1, is ellipse, 1 = parabola, > 1 = hyperbola).
125 * This function is also known as solveKeplersEquation()
128 * M: the mean anomaly
129 * e: the eccentricity
132 * the eccentric anomaly
134 ****************************************************************************/
135 double CelestialBody::sgCalcEccAnom(double M, double e)
138 eccAnom, E0, E1, diff;
140 eccAnom = M + e * sin(M) * (1.0 + e * cos (M));
141 // iterate to achieve a greater precision for larger eccentricities
147 E1 = E0 - (E0 - e * sin(E0) - M) / (1 - e *cos(E0));
148 diff = fabs(E0 - E1);
151 while (diff > (SGD_DEGREES_TO_RADIANS * 0.001));
157 /*****************************************************************************
158 * inline CelestialBody::CelestialBody
159 * public constructor for a generic celestialBody object.
160 * initializes the 6 primary orbital elements. The elements are:
161 * N: longitude of the ascending node
162 * i: inclination to the ecliptic
163 * w: argument of perihelion
164 * a: semi-major axis, or mean distance from the sun
167 * Each orbital element consists of a constant part and a variable part that
168 * gradually changes over time.
171 * the 13 arguments to the constructor constitute the first, constant
172 * ([NiwaeM]f) and the second variable ([NiwaeM]s) part of the orbital
173 * elements. The 13th argument is the current time. Note that the inclination
174 * is written with a capital (If, Is), because 'if' is a reserved word in the
175 * C/C++ programming language.
176 ***************************************************************************/
177 CelestialBody::CelestialBody(double Nf, double Ns,
178 double If, double Is,
179 double wf, double ws,
180 double af, double as,
181 double ef, double es,
182 double Mf, double Ms, double mjd)
184 NFirst = Nf; NSec = Ns;
185 iFirst = If; iSec = Is;
186 wFirst = wf; wSec = ws;
187 aFirst = af; aSec = as;
188 eFirst = ef; eSec = es;
189 MFirst = Mf; MSec = Ms;
190 updateOrbElements(mjd);
193 CelestialBody::CelestialBody(double Nf, double Ns,
194 double If, double Is,
195 double wf, double ws,
196 double af, double as,
197 double ef, double es,
198 double Mf, double Ms)
200 NFirst = Nf; NSec = Ns;
201 iFirst = If; iSec = Is;
202 wFirst = wf; wSec = ws;
203 aFirst = af; aSec = as;
204 eFirst = ef; eSec = es;
205 MFirst = Mf; MSec = Ms;
208 /****************************************************************************
209 * inline void CelestialBody::updateOrbElements(double mjd)
210 * given the current time, this private member calculates the actual
213 * Arguments: double mjd: the current modified julian date:
216 ***************************************************************************/
217 void CelestialBody::updateOrbElements(double mjd)
219 double actTime = sgCalcActTime(mjd);
220 M = SGD_DEGREES_TO_RADIANS * (MFirst + (MSec * actTime));
221 w = SGD_DEGREES_TO_RADIANS * (wFirst + (wSec * actTime));
222 N = SGD_DEGREES_TO_RADIANS * (NFirst + (NSec * actTime));
223 i = SGD_DEGREES_TO_RADIANS * (iFirst + (iSec * actTime));
224 e = eFirst + (eSec * actTime);
225 a = aFirst + (aSec * actTime);
228 /*****************************************************************************
229 * inline double CelestialBody::sgCalcActTime(double mjd)
230 * this private member function returns the offset in days from the epoch for
231 * wich the orbital elements are calculated (Jan, 1st, 2000).
233 * Argument: the current time
235 * return value: the (fractional) number of days until Jan 1, 2000.
236 ****************************************************************************/
237 double CelestialBody::sgCalcActTime(double mjd)
239 return (mjd - 36523.5);
242 /*****************************************************************************
243 * inline void CelestialBody::getPos(double* ra, double* dec)
244 * gives public access to Right Ascension and declination
246 ****************************************************************************/
247 void CelestialBody::getPos(double* ra, double* dec)
249 *ra = rightAscension;
253 /*****************************************************************************
254 * inline void CelestialBody::getPos(double* ra, double* dec, double* magnitude
255 * gives public acces to the current Right ascension, declination, and
257 ****************************************************************************/
258 void CelestialBody::getPos(double* ra, double* dec, double* magn)
260 *ra = rightAscension;