00001 /************************************************************************** 00002 * celestialBody.cxx 00003 * Written by Durk Talsma. Originally started October 1997, for distribution 00004 * with the FlightGear project. Version 2 was written in August and 00005 * September 1998. This code is based upon algorithms and data kindly 00006 * provided by Mr. Paul Schlyter. (pausch@saaf.se). 00007 * 00008 * This library is free software; you can redistribute it and/or 00009 * modify it under the terms of the GNU Library General Public 00010 * License as published by the Free Software Foundation; either 00011 * version 2 of the License, or (at your option) any later version. 00012 * 00013 * This library is distributed in the hope that it will be useful, 00014 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00015 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00016 * Library General Public License for more details. 00017 * 00018 * You should have received a copy of the GNU General Public License 00019 * along with this program; if not, write to the Free Software 00020 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 00021 * 00022 * $Id: celestialBody_8cxx_source.html,v 1.3 2010/02/23 22:10:07 curt Exp $ 00023 **************************************************************************/ 00024 00025 #include <simgear/debug/logstream.hxx> 00026 00027 #include <math.h> 00028 00029 #include "celestialBody.hxx" 00030 #include "star.hxx" 00031 00032 00033 /************************************************************************** 00034 * void CelestialBody::updatePosition(double mjd, Star *ourSun) 00035 * 00036 * Basically, this member function provides a general interface for 00037 * calculating the right ascension and declinaion. This function is 00038 * used for calculating the planetary positions. For the planets, an 00039 * overloaded member function is provided to additionally calculate the 00040 * planet's magnitude. 00041 * The sun and moon have their own overloaded updatePosition member, as their 00042 * position is calculated an a slightly different manner. 00043 * 00044 * arguments: 00045 * double mjd: provides the modified julian date. 00046 * Star *ourSun: the sun's position is needed to convert heliocentric 00047 * coordinates into geocentric coordinates. 00048 * 00049 * return value: none 00050 * 00051 *************************************************************************/ 00052 void CelestialBody::updatePosition(double mjd, Star *ourSun) 00053 { 00054 double eccAnom, v, ecl, actTime, 00055 xv, yv, xh, yh, zh, xg, yg, zg, xe, ye, ze; 00056 00057 updateOrbElements(mjd); 00058 actTime = sgCalcActTime(mjd); 00059 00060 // calcualate the angle bewteen ecliptic and equatorial coordinate system 00061 ecl = SGD_DEGREES_TO_RADIANS * (23.4393 - 3.563E-7 *actTime); 00062 00063 eccAnom = sgCalcEccAnom(M, e); //calculate the eccentric anomaly 00064 xv = a * (cos(eccAnom) - e); 00065 yv = a * (sqrt (1.0 - e*e) * sin(eccAnom)); 00066 v = atan2(yv, xv); // the planet's true anomaly 00067 r = sqrt (xv*xv + yv*yv); // the planet's distance 00068 00069 // calculate the planet's position in 3D space 00070 xh = r * (cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(i)); 00071 yh = r * (sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(i)); 00072 zh = r * (sin(v+w) * sin(i)); 00073 00074 // calculate the ecliptic longitude and latitude 00075 xg = xh + ourSun->getxs(); 00076 yg = yh + ourSun->getys(); 00077 zg = zh; 00078 00079 lonEcl = atan2(yh, xh); 00080 latEcl = atan2(zh, sqrt(xh*xh+yh*yh)); 00081 00082 xe = xg; 00083 ye = yg * cos(ecl) - zg * sin(ecl); 00084 ze = yg * sin(ecl) + zg * cos(ecl); 00085 rightAscension = atan2(ye, xe); 00086 declination = atan2(ze, sqrt(xe*xe + ye*ye)); 00087 /* SG_LOG(SG_GENERAL, SG_INFO, "Planet found at : " 00088 << rightAscension << " (ra), " << declination << " (dec)" ); */ 00089 00090 //calculate some variables specific to calculating the magnitude 00091 //of the planet 00092 R = sqrt (xg*xg + yg*yg + zg*zg); 00093 s = ourSun->getDistance(); 00094 00095 // It is possible from these calculations for the argument to acos 00096 // to exceed the valid range for acos(). So we do a little extra 00097 // checking. 00098 00099 double tmp = (r*r + R*R - s*s) / (2*r*R); 00100 if ( tmp > 1.0) { 00101 tmp = 1.0; 00102 } else if ( tmp < -1.0) { 00103 tmp = -1.0; 00104 } 00105 00106 FV = SGD_RADIANS_TO_DEGREES * acos( tmp ); 00107 } 00108 00109 /**************************************************************************** 00110 * double CelestialBody::sgCalcEccAnom(double M, double e) 00111 * this private member calculates the eccentric anomaly of a celestial body, 00112 * given its mean anomaly and eccentricity. 00113 * 00114 * -Mean anomaly: the approximate angle between the perihelion and the current 00115 * position. this angle increases uniformly with time. 00116 * 00117 * True anomaly: the actual angle between perihelion and current position. 00118 * 00119 * Eccentric anomaly: this is an auxilary angle, used in calculating the true 00120 * anomaly from the mean anomaly. 00121 * 00122 * -eccentricity. Indicates the amount in which the orbit deviates from a 00123 * circle (0 = circle, 0-1, is ellipse, 1 = parabola, > 1 = hyperbola). 00124 * 00125 * This function is also known as solveKeplersEquation() 00126 * 00127 * arguments: 00128 * M: the mean anomaly 00129 * e: the eccentricity 00130 * 00131 * return value: 00132 * the eccentric anomaly 00133 * 00134 ****************************************************************************/ 00135 double CelestialBody::sgCalcEccAnom(double M, double e) 00136 { 00137 double 00138 eccAnom, E0, E1, diff; 00139 00140 eccAnom = M + e * sin(M) * (1.0 + e * cos (M)); 00141 // iterate to achieve a greater precision for larger eccentricities 00142 if (e > 0.05) 00143 { 00144 E0 = eccAnom; 00145 do 00146 { 00147 E1 = E0 - (E0 - e * sin(E0) - M) / (1 - e *cos(E0)); 00148 diff = fabs(E0 - E1); 00149 E0 = E1; 00150 } 00151 while (diff > (SGD_DEGREES_TO_RADIANS * 0.001)); 00152 return E0; 00153 } 00154 return eccAnom; 00155 } 00156 00157 /***************************************************************************** 00158 * inline CelestialBody::CelestialBody 00159 * public constructor for a generic celestialBody object. 00160 * initializes the 6 primary orbital elements. The elements are: 00161 * N: longitude of the ascending node 00162 * i: inclination to the ecliptic 00163 * w: argument of perihelion 00164 * a: semi-major axis, or mean distance from the sun 00165 * e: eccenticity 00166 * M: mean anomaly 00167 * Each orbital element consists of a constant part and a variable part that 00168 * gradually changes over time. 00169 * 00170 * Argumetns: 00171 * the 13 arguments to the constructor constitute the first, constant 00172 * ([NiwaeM]f) and the second variable ([NiwaeM]s) part of the orbital 00173 * elements. The 13th argument is the current time. Note that the inclination 00174 * is written with a capital (If, Is), because 'if' is a reserved word in the 00175 * C/C++ programming language. 00176 ***************************************************************************/ 00177 CelestialBody::CelestialBody(double Nf, double Ns, 00178 double If, double Is, 00179 double wf, double ws, 00180 double af, double as, 00181 double ef, double es, 00182 double Mf, double Ms, double mjd) 00183 { 00184 NFirst = Nf; NSec = Ns; 00185 iFirst = If; iSec = Is; 00186 wFirst = wf; wSec = ws; 00187 aFirst = af; aSec = as; 00188 eFirst = ef; eSec = es; 00189 MFirst = Mf; MSec = Ms; 00190 updateOrbElements(mjd); 00191 } 00192 00193 CelestialBody::CelestialBody(double Nf, double Ns, 00194 double If, double Is, 00195 double wf, double ws, 00196 double af, double as, 00197 double ef, double es, 00198 double Mf, double Ms) 00199 { 00200 NFirst = Nf; NSec = Ns; 00201 iFirst = If; iSec = Is; 00202 wFirst = wf; wSec = ws; 00203 aFirst = af; aSec = as; 00204 eFirst = ef; eSec = es; 00205 MFirst = Mf; MSec = Ms; 00206 } 00207 00208 /**************************************************************************** 00209 * inline void CelestialBody::updateOrbElements(double mjd) 00210 * given the current time, this private member calculates the actual 00211 * orbital elements 00212 * 00213 * Arguments: double mjd: the current modified julian date: 00214 * 00215 * return value: none 00216 ***************************************************************************/ 00217 void CelestialBody::updateOrbElements(double mjd) 00218 { 00219 double actTime = sgCalcActTime(mjd); 00220 M = SGD_DEGREES_TO_RADIANS * (MFirst + (MSec * actTime)); 00221 w = SGD_DEGREES_TO_RADIANS * (wFirst + (wSec * actTime)); 00222 N = SGD_DEGREES_TO_RADIANS * (NFirst + (NSec * actTime)); 00223 i = SGD_DEGREES_TO_RADIANS * (iFirst + (iSec * actTime)); 00224 e = eFirst + (eSec * actTime); 00225 a = aFirst + (aSec * actTime); 00226 } 00227 00228 /***************************************************************************** 00229 * inline double CelestialBody::sgCalcActTime(double mjd) 00230 * this private member function returns the offset in days from the epoch for 00231 * wich the orbital elements are calculated (Jan, 1st, 2000). 00232 * 00233 * Argument: the current time 00234 * 00235 * return value: the (fractional) number of days until Jan 1, 2000. 00236 ****************************************************************************/ 00237 double CelestialBody::sgCalcActTime(double mjd) 00238 { 00239 return (mjd - 36523.5); 00240 } 00241 00242 /***************************************************************************** 00243 * inline void CelestialBody::getPos(double* ra, double* dec) 00244 * gives public access to Right Ascension and declination 00245 * 00246 ****************************************************************************/ 00247 void CelestialBody::getPos(double* ra, double* dec) 00248 { 00249 *ra = rightAscension; 00250 *dec = declination; 00251 } 00252 00253 /***************************************************************************** 00254 * inline void CelestialBody::getPos(double* ra, double* dec, double* magnitude 00255 * gives public acces to the current Right ascension, declination, and 00256 * magnitude 00257 ****************************************************************************/ 00258 void CelestialBody::getPos(double* ra, double* dec, double* magn) 00259 { 00260 *ra = rightAscension; 00261 *dec = declination; 00262 *magn = magnitude; 00263 } 00264 00265
1.6.1