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 #include <math.h>
  28
  29 #include "celestialBody.hxx"
  30 #include "star.hxx"
  31
  32
  33 /**************************************************************************
  34  * void CelestialBody::updatePosition(double mjd, Star *ourSun)
  35  *
  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
  40  * planet's magnitude.
  41  * The sun and moon have their own overloaded updatePosition member, as their
  42  * position is calculated an a slightly different manner.
  43  *
  44  * arguments:
  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.
  48  *
  49  * return value: none
  50  *
  51  *************************************************************************/
  52 void CelestialBody::updatePosition(double mjd, Star *ourSun)
  53 {
  54   double eccAnom, v, ecl, actTime,
  55     xv, yv, xh, yh, zh, xg, yg, zg, xe, ye, ze;
  56
  57   updateOrbElements(mjd);
  58   actTime = sgCalcActTime(mjd);
  59
  60   // calcualate the angle bewteen ecliptic and equatorial coordinate system
  61   ecl = SGD_DEGREES_TO_RADIANS * (23.4393 - 3.563E-7 *actTime);
  62
  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
  68
  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));
  73
  74   // calculate the ecliptic longitude and latitude
  75   xg = xh + ourSun->getxs();
  76   yg = yh + ourSun->getys();
  77   zg = zh;
  78
  79   lonEcl = atan2(yh, xh);
  80   latEcl = atan2(zh, sqrt(xh*xh+yh*yh));
  81
  82   xe = xg;
  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)" ); */
  89
  90   //calculate some variables specific to calculating the magnitude
  91   //of the planet
  92   R = sqrt (xg*xg + yg*yg + zg*zg);
  93   s = ourSun->getDistance();
  94
  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
  97   // checking.
  98
  99   double tmp = (r*r + R*R - s*s) / (2*r*R);
 100   if ( tmp > 1.0) {
 101       tmp = 1.0;
 102   } else if ( tmp < -1.0) {
 103       tmp = -1.0;
 104   }
 105
 106   FV = SGD_RADIANS_TO_DEGREES * acos( tmp );
 107 }
 108
 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.
 113  *
 114  * -Mean anomaly: the approximate angle between the perihelion and the current
 115  *  position. this angle increases uniformly with time.
 116  *
 117  * True anomaly: the actual angle between perihelion and current position.
 118  *
 119  * Eccentric anomaly: this is an auxilary angle, used in calculating the true
 120  * anomaly from the mean anomaly.
 121  *
 122  * -eccentricity. Indicates the amount in which the orbit deviates from a
 123  *  circle (0 = circle, 0-1, is ellipse, 1 = parabola, > 1 = hyperbola).
 124  *
 125  * This function is also known as solveKeplersEquation()
 126  *
 127  * arguments:
 128  * M: the mean anomaly
 129  * e: the eccentricity
 130  *
 131  * return value:
 132  * the eccentric anomaly
 133  *
 134  ****************************************************************************/
 135 double CelestialBody::sgCalcEccAnom(double M, double e)
 136 {
 137   double
 138     eccAnom, E0, E1, diff;
 139
 140   eccAnom = M + e * sin(M) * (1.0 + e * cos (M));
 141   // iterate to achieve a greater precision for larger eccentricities
 142   if (e > 0.05)
 143     {
 144       E0 = eccAnom;
 145       do
 146         {
 147           E1 = E0 - (E0 - e * sin(E0) - M) / (1 - e *cos(E0));
 148           diff = fabs(E0 - E1);
 149           E0 = E1;
 150         }
 151       while (diff > (SGD_DEGREES_TO_RADIANS * 0.001));
 152       return E0;
 153     }
 154   return eccAnom;
 155 }
 156
 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
 165  * e: eccenticity
 166  * M: mean anomaly
 167  * Each orbital element consists of a constant part and a variable part that
 168  * gradually changes over time.
 169  *
 170  * Argumetns:
 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)
 183 {
 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);
 191 }
 192
 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)
 199 {
 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;
 206 }
 207
 208 /****************************************************************************
 209  * inline void CelestialBody::updateOrbElements(double mjd)
 210  * given the current time, this private member calculates the actual
 211  * orbital elements
 212  *
 213  * Arguments: double mjd: the current modified julian date:
 214  *
 215  * return value: none
 216  ***************************************************************************/
 217 void CelestialBody::updateOrbElements(double mjd)
 218 {
 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);
 226 }
 227
 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).
 232  *
 233  * Argument: the current time
 234  *
 235  * return value: the (fractional) number of days until Jan 1, 2000.
 236  ****************************************************************************/
 237 double CelestialBody::sgCalcActTime(double mjd)
 238 {
 239   return (mjd - 36523.5);
 240 }
 241
 242 /*****************************************************************************
 243  * inline void CelestialBody::getPos(double* ra, double* dec)
 244  * gives public access to Right Ascension and declination
 245  *
 246  ****************************************************************************/
 247 void CelestialBody::getPos(double* ra, double* dec)
 248 {
 249   *ra  = rightAscension;
 250   *dec = declination;
 251 }
 252
 253 /*****************************************************************************
 254  * inline void CelestialBody::getPos(double* ra, double* dec, double* magnitude
 255  * gives public acces to the current Right ascension, declination, and
 256  * magnitude
 257  ****************************************************************************/
 258 void CelestialBody::getPos(double* ra, double* dec, double* magn)
 259 {
 260   *ra = rightAscension;
 261   *dec = declination;
 262   *magn = magnitude;
 263 }
 264
 265