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