simgear/ephemeris/moonpos.cxx

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