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 #if 0
  61   int width, height;
  62
  63   FG_LOG( FG_GENERAL, FG_INFO, "Initializing Moon Texture");
  64 #ifdef GL_VERSION_1_1
  65   xglGenTextures(1, &moon_texid);
  66   xglBindTexture(GL_TEXTURE_2D, moon_texid);
  67 #elif GL_EXT_texture_object
  68   xglGenTexturesEXT(1, &moon_texid);
  69   xglBindTextureEXT(GL_TEXTURE_2D, moon_texid);
  70 #else
  71 #  error port me
  72 #endif
  73
  74   glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
  75   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
  76   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
  77
  78   // load in the texture data
  79   FGPath tpath( current_options.get_fg_root() );
  80   tpath.append( "Textures" );
  81   tpath.append( "moon.rgb" );
  82
  83   if ( (moon_texbuf = read_rgb_texture(tpath.c_str(), &width, &height))
  84        == NULL )
  85   {
  86     // Try compressed
  87     FGPath fg_tpath = tpath;
  88     fg_tpath.append( ".gz" );
  89     if ( (moon_texbuf = read_rgb_texture(fg_tpath.c_str(), &width, &height))
  90          == NULL )
  91     {
  92         FG_LOG( FG_GENERAL, FG_ALERT,
  93                 "Error in loading moon texture " << tpath.str() );
  94         exit(-1);
  95     }
  96   }
  97
  98   glTexImage2D( GL_TEXTURE_2D,
  99                 0,
 100                 GL_RGB,
 101                 256, 256,
 102                 0,
 103                 GL_RGB, GL_UNSIGNED_BYTE,
 104                 moon_texbuf);
 105
 106   // setup the halo texture
 107   FG_LOG( FG_GENERAL, FG_INFO, "Initializing Moon Texture");
 108 #ifdef GL_VERSION_1_1
 109   xglGenTextures(1, &moon_halotexid);
 110   xglBindTexture(GL_TEXTURE_2D, moon_halotexid);
 111 #elif GL_EXT_texture_object
 112   xglGenTexturesEXT(1, &moon_halotexid);
 113   xglBindTextureEXT(GL_TEXTURE_2D, moon_halotexid);
 114 #else
 115 #  error port me
 116 #endif
 117
 118   glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
 119   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
 120   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
 121   setHalo();
 122   glTexImage2D( GL_TEXTURE_2D,
 123                 0,
 124                 GL_RGBA,
 125                 256, 256,
 126                 0,
 127                 GL_RGBA, GL_UNSIGNED_BYTE,
 128                 moon_halotexbuf);
 129   moonObject = gluNewQuadric();
 130 #endif
 131 }
 132
 133 Moon::Moon() :
 134   CelestialBody(125.1228, -0.0529538083,
 135                 5.1454,    0.00000,
 136                 318.0634,  0.1643573223,
 137                 60.266600, 0.000000,
 138                 0.054900,  0.000000,
 139                 115.3654,  13.0649929509)
 140 {
 141 }
 142
 143
 144 Moon::~Moon()
 145 {
 146     //delete moonObject;
 147     // delete moon_texbuf;
 148     // delete moon_halotexbuf;
 149 }
 150
 151
 152 #if 0
 153 static int texWidth = 256;      /* 64x64 is plenty */
 154
 155 void Moon::setHalo()
 156 {
 157   int texSize;
 158   //void *textureBuf;
 159   GLubyte *p;
 160   int i,j;
 161   double radius;
 162
 163   texSize = texWidth*texWidth;
 164
 165   moon_halotexbuf = new GLubyte[texSize*4];
 166   if (!moon_halotexbuf)
 167     return;  // Ugly!
 168
 169   p = moon_halotexbuf;
 170
 171   radius = (double)(texWidth / 2);
 172
 173   for (i=0; i < texWidth; i++) {
 174     for (j=0; j < texWidth; j++) {
 175       double x, y, d;
 176
 177       *p = 0xff;
 178       *(p+1) = 0xff;
 179       *(p+2) = 0xff;
 180
 181       x = fabs((double)(i - (texWidth / 2)));
 182       y = fabs((double)(j - (texWidth / 2)));
 183
 184       d = sqrt((x * x) + (y * y));
 185       if (d < radius) {
 186           double t = 1.0 - (d / radius); // t is 1.0 at center, 0.0 at edge */
 187           // inverse square looks nice
 188           *(p+3) = (int)((double) 0x20 * (t*t));
 189       } else {
 190           *(p+3) = 0x00;
 191       }
 192       p += 4;
 193     }
 194   }
 195   //gluBuild2DMipmaps(GL_TEXTURE_2D, 1, texWidth, texWidth,
 196   //        GL_LUMINANCE,
 197   //        GL_UNSIGNED_BYTE, textureBuf);
 198   //free(textureBuf);
 199 }
 200 #endif
 201
 202
 203 /*****************************************************************************
 204  * void Moon::updatePosition(FGTime *t, Star *ourSun)
 205  * this member function calculates the actual topocentric position (i.e.)
 206  * the position of the moon as seen from the current position on the surface
 207  * of the moon.
 208  ****************************************************************************/
 209 void Moon::updatePosition(FGTime *t, double lat, Star *ourSun)
 210 {
 211   double
 212     eccAnom, ecl, actTime,
 213     xv, yv, v, r, xh, yh, zh, xg, yg, zg, xe, ye, ze,
 214     Ls, Lm, D, F, mpar, gclat, rho, HA, g,
 215     geoRa, geoDec;
 216
 217   fgAIRCRAFT *air;
 218   FGInterface *f;
 219
 220   air = &current_aircraft;
 221   f = air->fdm_state;
 222
 223   updateOrbElements(t);
 224   actTime = fgCalcActTime(t);
 225
 226   // calculate the angle between ecliptic and equatorial coordinate system
 227   // in Radians
 228   ecl = ((DEG_TO_RAD * 23.4393) - (DEG_TO_RAD * 3.563E-7) * actTime);
 229   eccAnom = fgCalcEccAnom(M, e);  // Calculate the eccentric anomaly
 230   xv = a * (cos(eccAnom) - e);
 231   yv = a * (sqrt(1.0 - e*e) * sin(eccAnom));
 232   v = atan2(yv, xv);               // the moon's true anomaly
 233   r = sqrt (xv*xv + yv*yv);       // and its distance
 234
 235   // estimate the geocentric rectangular coordinates here
 236   xh = r * (cos(N) * cos (v+w) - sin (N) * sin(v+w) * cos(i));
 237   yh = r * (sin(N) * cos (v+w) + cos (N) * sin(v+w) * cos(i));
 238   zh = r * (sin(v+w) * sin(i));
 239
 240   // calculate the ecliptic latitude and longitude here
 241   lonEcl = atan2 (yh, xh);
 242   latEcl = atan2(zh, sqrt(xh*xh + yh*yh));
 243
 244   /* Calculate a number of perturbatioin, i.e. disturbances caused by the
 245    * gravitational infuence of the sun and the other major planets.
 246    * The largest of these even have a name */
 247   Ls = ourSun->getM() + ourSun->getw();
 248   Lm = M + w + N;
 249   D = Lm - Ls;
 250   F = Lm - N;
 251
 252   lonEcl += DEG_TO_RAD * (-1.274 * sin (M - 2*D)
 253                           +0.658 * sin (2*D)
 254                           -0.186 * sin(ourSun->getM())
 255                           -0.059 * sin(2*M - 2*D)
 256                           -0.057 * sin(M - 2*D + ourSun->getM())
 257                           +0.053 * sin(M + 2*D)
 258                           +0.046 * sin(2*D - ourSun->getM())
 259                           +0.041 * sin(M - ourSun->getM())
 260                           -0.035 * sin(D)
 261                           -0.031 * sin(M + ourSun->getM())
 262                           -0.015 * sin(2*F - 2*D)
 263                           +0.011 * sin(M - 4*D)
 264                           );
 265   latEcl += DEG_TO_RAD * (-0.173 * sin(F-2*D)
 266                           -0.055 * sin(M - F - 2*D)
 267                           -0.046 * sin(M + F - 2*D)
 268                           +0.033 * sin(F + 2*D)
 269                           +0.017 * sin(2*M + F)
 270                           );
 271   r += (-0.58 * cos(M - 2*D)
 272         -0.46 * cos(2*D)
 273         );
 274   // FG_LOG(FG_GENERAL, FG_INFO, "Running moon update");
 275   xg = r * cos(lonEcl) * cos(latEcl);
 276   yg = r * sin(lonEcl) * cos(latEcl);
 277   zg = r *               sin(latEcl);
 278
 279   xe = xg;
 280   ye = yg * cos(ecl) -zg * sin(ecl);
 281   ze = yg * sin(ecl) +zg * cos(ecl);
 282
 283   geoRa  = atan2(ye, xe);
 284   geoDec = atan2(ze, sqrt(xe*xe + ye*ye));
 285
 286   /* FG_LOG( FG_GENERAL, FG_INFO,
 287           "(geocentric) geoRa = (" << (RAD_TO_DEG * geoRa)
 288           << "), geoDec= (" << (RAD_TO_DEG * geoDec) << ")" ); */
 289
 290
 291   // Given the moon's geocentric ra and dec, calculate its
 292   // topocentric ra and dec. i.e. the position as seen from the
 293   // surface of the earth, instead of the center of the earth
 294
 295   // First calculate the moon's parrallax, that is, the apparent size of the
 296   // (equatorial) radius of the earth, as seen from the moon
 297   mpar = asin ( 1 / r);
 298   // FG_LOG( FG_GENERAL, FG_INFO, "r = " << r << " mpar = " << mpar );
 299   // FG_LOG( FG_GENERAL, FG_INFO, "lat = " << f->get_Latitude() );
 300
 301   gclat = lat - 0.003358 *
 302       sin (2 * DEG_TO_RAD * lat );
 303   // FG_LOG( FG_GENERAL, FG_INFO, "gclat = " << gclat );
 304
 305   rho = 0.99883 + 0.00167 * cos(2 * DEG_TO_RAD * lat);
 306   // FG_LOG( FG_GENERAL, FG_INFO, "rho = " << rho );
 307
 308   if (geoRa < 0)
 309     geoRa += (2*FG_PI);
 310
 311   HA = t->getLst() - (3.8197186 * geoRa);
 312   /* FG_LOG( FG_GENERAL, FG_INFO, "t->getLst() = " << t->getLst()
 313           << " HA = " << HA ); */
 314
 315   g = atan (tan(gclat) / cos ((HA / 3.8197186)));
 316   // FG_LOG( FG_GENERAL, FG_INFO, "g = " << g );
 317
 318   rightAscension = geoRa - mpar * rho * cos(gclat) * sin(HA) / cos (geoDec);
 319   declination = geoDec - mpar * rho * sin (gclat) * sin (g - geoDec) / sin(g);
 320
 321   /* FG_LOG( FG_GENERAL, FG_INFO,
 322           "Ra = (" << (RAD_TO_DEG *rightAscension)
 323           << "), Dec= (" << (RAD_TO_DEG *declination) << ")" ); */
 324 }
 325
 326
 327 #if 0
 328 /************************************************************************
 329  * void Moon::newImage()
 330  *
 331  * This function regenerates a new visual image of the moon, which is added to
 332  * solarSystem display list.
 333  *
 334  * Arguments: Right Ascension and declination
 335  *
 336  * return value: none
 337  **************************************************************************/
 338 void Moon::newImage()
 339 {
 340   fgLIGHT *l = &cur_light_params;
 341   float moon_angle = l->moon_angle;
 342
 343   /*double x_2, x_4, x_8, x_10;
 344   GLfloat ambient;
 345   GLfloat amb[4];*/
 346   int moonSize = 550;
 347
 348   GLfloat moonColor[4] = {0.85, 0.75, 0.35, 1.0};
 349   GLfloat black[4] = {0.0, 0.0, 0.0, 1.0};
 350   GLfloat white[4] = {1.0, 1.0, 1.0, 0.0};
 351
 352   if( moon_angle*RAD_TO_DEG < 100 )
 353     {
 354       FG_LOG( FG_ASTRO, FG_INFO, "Generating Moon Image" );
 355
 356       xglPushMatrix();
 357       {
 358         xglRotatef(((RAD_TO_DEG * rightAscension)- 90.0), 0.0, 0.0, 1.0);
 359         xglRotatef((RAD_TO_DEG * declination), 1.0, 0.0, 0.0);
 360
 361         FG_LOG( FG_GENERAL, FG_INFO,
 362                 "Ra = (" << (RAD_TO_DEG *rightAscension)
 363                 << "), Dec= (" << (RAD_TO_DEG *declination) << ")" );
 364         xglTranslatef(0.0, 60000.0, 0.0);
 365         glEnable(GL_BLEND);     // BLEND ENABLED
 366
 367         // Draw the halo...
 368         if (current_options.get_textures())
 369           {
 370             // glBlendFunc(GL_SRC_ALPHA, GL_ONE);
 371             glBlendFunc ( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA ) ;
 372             glEnable(GL_TEXTURE_2D); // TEXTURE ENABLED
 373             glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
 374             glBindTexture(GL_TEXTURE_2D, moon_halotexid);
 375
 376             glBegin(GL_QUADS);
 377             glTexCoord2f(0.0f, 0.0f); glVertex3f(-5000, 0.0, -5000);
 378             glTexCoord2f(1.0f, 0.0f); glVertex3f( 5000, 0.0, -5000);
 379             glTexCoord2f(1.0f, 1.0f); glVertex3f( 5000, 0.0,  5000);
 380             glTexCoord2f(0.0f, 1.0f); glVertex3f(-5000, 0.0,  5000);
 381             glEnd();
 382           }
 383
 384         xglEnable(GL_LIGHTING); // LIGHTING ENABLED
 385         xglEnable( GL_LIGHT0 );
 386         // set lighting parameters
 387         xglLightfv(GL_LIGHT0, GL_AMBIENT, white );
 388         xglLightfv(GL_LIGHT0, GL_DIFFUSE, white );
 389         // Enable( GL_CULL_FACE );
 390         xglMaterialfv(GL_FRONT, GL_AMBIENT, black);
 391         xglMaterialfv(GL_FRONT, GL_DIFFUSE, moonColor);
 392         //glEnable(GL_TEXTURE_2D);
 393
 394         glBlendFunc(GL_ONE, GL_ONE);
 395
 396         //glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
 397         // Draw the moon-proper
 398
 399         if (current_options.get_textures())
 400           {
 401             glBindTexture(GL_TEXTURE_2D, moon_texid);
 402             gluQuadricTexture(moonObject, GL_TRUE );
 403           }
 404         gluSphere(moonObject,  moonSize, 12, 12 );
 405         glDisable(GL_TEXTURE_2D); // TEXTURE DISABLED
 406         glDisable(GL_BLEND);    // BLEND DISABLED
 407       }
 408       xglPopMatrix();
 409       glDisable(GL_LIGHTING);   // Lighting Disabled.
 410
 411     }
 412   else
 413     {
 414     }
 415 }
 416 #endif