simgear/ephemeris/celestialBody.cxx

   1 /**************************************************************************
   2  * celestialBody.cxx
   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).
   7  *
   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.
  12  *
  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.
  17  *
  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.
  21  *
  22  * $Id$
  23  **************************************************************************/
  24
  25 #include <simgear/debug/logstream.hxx>
  26
  27 #ifdef SG_MATH_EXCEPTION_CLASH
  28 #  define exception c_exception
  29 #endif
  30 #include <math.h>
  31
  32 #include "celestialBody.hxx"
  33 #include "star.hxx"
  34
  35
  36 /**************************************************************************
  37  * void CelestialBody::updatePosition(double mjd, Star *ourSun)
  38  *
  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
  43  * planet's magnitude.
  44  * The sun and moon have their own overloaded updatePosition member, as their
  45  * position is calculated an a slightly different manner.
  46  *
  47  * arguments:
  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.
  51  *
  52  * return value: none
  53  *
  54  *************************************************************************/
  55 void CelestialBody::updatePosition(double mjd, Star *ourSun)
  56 {
  57   double eccAnom, v, ecl, actTime,
  58     xv, yv, xh, yh, zh, xg, yg, zg, xe, ye, ze;
  59
  60   updateOrbElements(mjd);
  61   actTime = sgCalcActTime(mjd);
  62
  63   // calcualate the angle bewteen ecliptic and equatorial coordinate system
  64   ecl = SGD_DEGREES_TO_RADIANS * (23.4393 - 3.563E-7 *actTime);
  65
  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
  71
  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));
  76
  77   // calculate the ecliptic longitude and latitude
  78   xg = xh + ourSun->getxs();
  79   yg = yh + ourSun->getys();
  80   zg = zh;
  81
  82   lonEcl = atan2(yh, xh);
  83   latEcl = atan2(zh, sqrt(xh*xh+yh*yh));
  84
  85   xe = xg;
  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)" ); */
  92
  93   //calculate some variables specific to calculating the magnitude
  94   //of the planet
  95   R = sqrt (xg*xg + yg*yg + zg*zg);
  96   s = ourSun->getDistance();
  97
  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
 100   // checking.
 101
 102   double tmp = (r*r + R*R - s*s) / (2*r*R);
 103   if ( tmp > 1.0) {
 104       tmp = 1.0;
 105   } else if ( tmp < -1.0) {
 106       tmp = -1.0;
 107   }
 108
 109   FV = SGD_RADIANS_TO_DEGREES * acos( tmp );
 110 }
 111
 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.
 116  *
 117  * -Mean anomaly: the approximate angle between the perihelion and the current
 118  *  position. this angle increases uniformly with time.
 119  *
 120  * True anomaly: the actual angle between perihelion and current position.
 121  *
 122  * Eccentric anomaly: this is an auxilary angle, used in calculating the true
 123  * anomaly from the mean anomaly.
 124  *
 125  * -eccentricity. Indicates the amount in which the orbit deviates from a
 126  *  circle (0 = circle, 0-1, is ellipse, 1 = parabola, > 1 = hyperbola).
 127  *
 128  * This function is also known as solveKeplersEquation()
 129  *
 130  * arguments:
 131  * M: the mean anomaly
 132  * e: the eccentricity
 133  *
 134  * return value:
 135  * the eccentric anomaly
 136  *
 137  ****************************************************************************/
 138 double CelestialBody::sgCalcEccAnom(double M, double e)
 139 {
 140   double
 141     eccAnom, E0, E1, diff;
 142
 143   eccAnom = M + e * sin(M) * (1.0 + e * cos (M));
 144   // iterate to achieve a greater precision for larger eccentricities
 145   if (e > 0.05)
 146     {
 147       E0 = eccAnom;
 148       do
 149         {
 150           E1 = E0 - (E0 - e * sin(E0) - M) / (1 - e *cos(E0));
 151           diff = fabs(E0 - E1);
 152           E0 = E1;
 153         }
 154       while (diff > (SGD_DEGREES_TO_RADIANS * 0.001));
 155       return E0;
 156     }
 157   return eccAnom;
 158 }
 159
 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
 168  * e: eccenticity
 169  * M: mean anomaly
 170  * Each orbital element consists of a constant part and a variable part that
 171  * gradually changes over time.
 172  *
 173  * Argumetns:
 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)
 186 {
 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);
 194 }
 195
 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)
 202 {
 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;
 209 }
 210
 211 /****************************************************************************
 212  * inline void CelestialBody::updateOrbElements(double mjd)
 213  * given the current time, this private member calculates the actual
 214  * orbital elements
 215  *
 216  * Arguments: double mjd: the current modified julian date:
 217  *
 218  * return value: none
 219  ***************************************************************************/
 220 void CelestialBody::updateOrbElements(double mjd)
 221 {
 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);
 229 }
 230
 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).
 235  *
 236  * Argument: the current time
 237  *
 238  * return value: the (fractional) number of days until Jan 1, 2000.
 239  ****************************************************************************/
 240 double CelestialBody::sgCalcActTime(double mjd)
 241 {
 242   return (mjd - 36523.5);
 243 }
 244
 245 /*****************************************************************************
 246  * inline void CelestialBody::getPos(double* ra, double* dec)
 247  * gives public access to Right Ascension and declination
 248  *
 249  ****************************************************************************/
 250 void CelestialBody::getPos(double* ra, double* dec)
 251 {
 252   *ra  = rightAscension;
 253   *dec = declination;
 254 }
 255
 256 /*****************************************************************************
 257  * inline void CelestialBody::getPos(double* ra, double* dec, double* magnitude
 258  * gives public acces to the current Right ascension, declination, and
 259  * magnitude
 260  ****************************************************************************/
 261 void CelestialBody::getPos(double* ra, double* dec, double* magn)
 262 {
 263   *ra = rightAscension;
 264   *dec = declination;
 265   *magn = magnitude;
 266 }
 267
 268