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