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