simgear/ephemeris/moon.cxx

   1 /**************************************************************************
   2  * moon.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 program is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU General Public License as
  10  * published by the Free Software Foundation; either version 2 of the
  11  * License, or (at your option) any later version.
  12  *
  13  * This program is distributed in the hope that it will be useful, but
  14  * WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16  * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
  21  *
  22  * $Id$
  23  **************************************************************************/
  24
  25
  26 #include <string.h>
  27
  28 #include <simgear/debug/logstream.hxx>
  29 #include <simgear/misc/fgpath.hxx>
  30
  31 #include <Main/options.hxx>
  32 #include <Objects/texload.h>
  33
  34 #ifdef __BORLANDC__
  35 #  define exception c_exception
  36 #endif
  37 #include <math.h>
  38
  39 #include <FDM/flight.hxx>
  40
  41 #include "moon.hxx"
  42
  43
  44 /*************************************************************************
  45  * Moon::Moon(FGTime *t)
  46  * Public constructor for class Moon. Initializes the orbital elements and
  47  * sets up the moon texture.
  48  * Argument: The current time.
  49  * the hard coded orbital elements for Moon are passed to
  50  * CelestialBody::CelestialBody();
  51  ************************************************************************/
  52 Moon::Moon(FGTime *t) :
  53   CelestialBody(125.1228, -0.0529538083,
  54                 5.1454,    0.00000,
  55                 318.0634,  0.1643573223,
  56                 60.266600, 0.000000,
  57                 0.054900,  0.000000,
  58                 115.3654,  13.0649929509, t)
  59 {
  60 }
  61
  62 Moon::Moon() :
  63   CelestialBody(125.1228, -0.0529538083,
  64                 5.1454,    0.00000,
  65                 318.0634,  0.1643573223,
  66                 60.266600, 0.000000,
  67                 0.054900,  0.000000,
  68                 115.3654,  13.0649929509)
  69 {
  70 }
  71
  72
  73 Moon::~Moon()
  74 {
  75 }
  76
  77
  78 /*****************************************************************************
  79  * void Moon::updatePosition(FGTime *t, Star *ourSun)
  80  * this member function calculates the actual topocentric position (i.e.)
  81  * the position of the moon as seen from the current position on the surface
  82  * of the moon.
  83  ****************************************************************************/
  84 void Moon::updatePosition(FGTime *t, double lat, Star *ourSun)
  85 {
  86   double
  87     eccAnom, ecl, actTime,
  88     xv, yv, v, r, xh, yh, zh, xg, yg, zg, xe, ye, ze,
  89     Ls, Lm, D, F, mpar, gclat, rho, HA, g,
  90     geoRa, geoDec;
  91
  92   fgAIRCRAFT *air;
  93   FGInterface *f;
  94
  95   air = &current_aircraft;
  96   f = air->fdm_state;
  97
  98   updateOrbElements(t);
  99   actTime = fgCalcActTime(t);
 100
 101   // calculate the angle between ecliptic and equatorial coordinate system
 102   // in Radians
 103   ecl = ((DEG_TO_RAD * 23.4393) - (DEG_TO_RAD * 3.563E-7) * actTime);
 104   eccAnom = fgCalcEccAnom(M, e);  // Calculate the eccentric anomaly
 105   xv = a * (cos(eccAnom) - e);
 106   yv = a * (sqrt(1.0 - e*e) * sin(eccAnom));
 107   v = atan2(yv, xv);               // the moon's true anomaly
 108   r = sqrt (xv*xv + yv*yv);       // and its distance
 109
 110   // estimate the geocentric rectangular coordinates here
 111   xh = r * (cos(N) * cos (v+w) - sin (N) * sin(v+w) * cos(i));
 112   yh = r * (sin(N) * cos (v+w) + cos (N) * sin(v+w) * cos(i));
 113   zh = r * (sin(v+w) * sin(i));
 114
 115   // calculate the ecliptic latitude and longitude here
 116   lonEcl = atan2 (yh, xh);
 117   latEcl = atan2(zh, sqrt(xh*xh + yh*yh));
 118
 119   /* Calculate a number of perturbatioin, i.e. disturbances caused by the
 120    * gravitational infuence of the sun and the other major planets.
 121    * The largest of these even have a name */
 122   Ls = ourSun->getM() + ourSun->getw();
 123   Lm = M + w + N;
 124   D = Lm - Ls;
 125   F = Lm - N;
 126
 127   lonEcl += DEG_TO_RAD * (-1.274 * sin (M - 2*D)
 128                           +0.658 * sin (2*D)
 129                           -0.186 * sin(ourSun->getM())
 130                           -0.059 * sin(2*M - 2*D)
 131                           -0.057 * sin(M - 2*D + ourSun->getM())
 132                           +0.053 * sin(M + 2*D)
 133                           +0.046 * sin(2*D - ourSun->getM())
 134                           +0.041 * sin(M - ourSun->getM())
 135                           -0.035 * sin(D)
 136                           -0.031 * sin(M + ourSun->getM())
 137                           -0.015 * sin(2*F - 2*D)
 138                           +0.011 * sin(M - 4*D)
 139                           );
 140   latEcl += DEG_TO_RAD * (-0.173 * sin(F-2*D)
 141                           -0.055 * sin(M - F - 2*D)
 142                           -0.046 * sin(M + F - 2*D)
 143                           +0.033 * sin(F + 2*D)
 144                           +0.017 * sin(2*M + F)
 145                           );
 146   r += (-0.58 * cos(M - 2*D)
 147         -0.46 * cos(2*D)
 148         );
 149   // FG_LOG(FG_GENERAL, FG_INFO, "Running moon update");
 150   xg = r * cos(lonEcl) * cos(latEcl);
 151   yg = r * sin(lonEcl) * cos(latEcl);
 152   zg = r *               sin(latEcl);
 153
 154   xe = xg;
 155   ye = yg * cos(ecl) -zg * sin(ecl);
 156   ze = yg * sin(ecl) +zg * cos(ecl);
 157
 158   geoRa  = atan2(ye, xe);
 159   geoDec = atan2(ze, sqrt(xe*xe + ye*ye));
 160
 161   /* FG_LOG( FG_GENERAL, FG_INFO,
 162           "(geocentric) geoRa = (" << (RAD_TO_DEG * geoRa)
 163           << "), geoDec= (" << (RAD_TO_DEG * geoDec) << ")" ); */
 164
 165
 166   // Given the moon's geocentric ra and dec, calculate its
 167   // topocentric ra and dec. i.e. the position as seen from the
 168   // surface of the earth, instead of the center of the earth
 169
 170   // First calculate the moon's parrallax, that is, the apparent size of the
 171   // (equatorial) radius of the earth, as seen from the moon
 172   mpar = asin ( 1 / r);
 173   // FG_LOG( FG_GENERAL, FG_INFO, "r = " << r << " mpar = " << mpar );
 174   // FG_LOG( FG_GENERAL, FG_INFO, "lat = " << f->get_Latitude() );
 175
 176   gclat = lat - 0.003358 *
 177       sin (2 * DEG_TO_RAD * lat );
 178   // FG_LOG( FG_GENERAL, FG_INFO, "gclat = " << gclat );
 179
 180   rho = 0.99883 + 0.00167 * cos(2 * DEG_TO_RAD * lat);
 181   // FG_LOG( FG_GENERAL, FG_INFO, "rho = " << rho );
 182
 183   if (geoRa < 0)
 184     geoRa += (2*FG_PI);
 185
 186   HA = t->getLst() - (3.8197186 * geoRa);
 187   /* FG_LOG( FG_GENERAL, FG_INFO, "t->getLst() = " << t->getLst()
 188           << " HA = " << HA ); */
 189
 190   g = atan (tan(gclat) / cos ((HA / 3.8197186)));
 191   // FG_LOG( FG_GENERAL, FG_INFO, "g = " << g );
 192
 193   rightAscension = geoRa - mpar * rho * cos(gclat) * sin(HA) / cos (geoDec);
 194   declination = geoDec - mpar * rho * sin (gclat) * sin (g - geoDec) / sin(g);
 195
 196   /* FG_LOG( FG_GENERAL, FG_INFO,
 197           "Ra = (" << (RAD_TO_DEG *rightAscension)
 198           << "), Dec= (" << (RAD_TO_DEG *declination) << ")" ); */
 199 }