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