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