Astro/moon.c

   1 /**************************************************************************
   2  * moon.c
   3  * Written by Durk Talsma. Started October 1997, for the flight gear project.
   4  *
   5  * This program is free software; you can redistribute it and/or
   6  * modify it under the terms of the GNU General Public License as
   7  * published by the Free Software Foundation; either version 2 of the
   8  * License, or (at your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18  *
  19  * $Id$
  20  * (Log is kept at end of this file)
  21  **************************************************************************/
  22
  23
  24 #include <math.h>
  25 #include <GL/glut.h>
  26 #include <XGL/xgl.h>
  27
  28 #include <Astro/orbits.h>
  29 #include <Astro/moon.h>
  30
  31 #include <Aircraft/aircraft.h>
  32 #include <Include/fg_constants.h>
  33 #include <Include/general.h>
  34 #include <Main/views.h>
  35 #include <Time/fg_time.h>
  36 #include <Main/fg_debug.h>
  37
  38 struct CelestialCoord moonPos;
  39
  40 static float xMoon, yMoon, zMoon;
  41 static GLint moon;
  42
  43 /*
  44 static GLfloat vdata[12][3] =
  45 {
  46    {-X, 0.0, Z }, { X, 0.0, Z }, {-X, 0.0, -Z}, {X, 0.0, -Z },
  47    { 0.0, Z, X }, { 0.0, Z, -X}, {0.0, -Z, -X}, {0.0, -Z, -X},
  48    { Z, X, 0.0 }, { -Z, X, 0.0}, {Z, -X, 0.0 }, {-Z, -X, 0.0}
  49 };
  50
  51 static GLuint tindices[20][3] =
  52 {
  53    {0,4,1}, {0,9,4}, {9,5,4}, {4,5,8}, {4,8,1},
  54    {8,10,1}, {8,3,10}, {5,3,8}, {5,2,3}, {2,7,3},
  55    {7,10,3}, {7,6,10}, {7,11,6}, {11,0,6}, {0,1,6},
  56    {6,1,10}, {9,0,11}, {9,11,2}, {9,2,5}, {7,2,11}
  57 };*/
  58
  59 /* -------------------------------------------------------------
  60       This section contains the code that generates a yellow
  61       Icosahedron. It's under development... (of Course)
  62 ______________________________________________________________*/
  63 /*
  64 void NormalizeVector(float v[3])
  65 {
  66    GLfloat d = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
  67    if (d == 0.0)
  68    {
  69       printf("zero length vector\n");
  70       return;
  71    }
  72    v[0] /= d;
  73    v[1] /= d;
  74    v[2] /= d;
  75 }
  76
  77 void drawTriangle(float *v1, float *v2, float *v3)
  78 {
  79    xglBegin(GL_POLYGON);
  80    //xglBegin(GL_POINTS);
  81       xglNormal3fv(v1);
  82       xglVertex3fv(v1);
  83       xglNormal3fv(v2);
  84       xglVertex3fv(v2);
  85       xglNormal3fv(v3);
  86       xglVertex3fv(v3);
  87    xglEnd();
  88 }
  89
  90 void subdivide(float *v1, float *v2, float *v3, long depth)
  91 {
  92    GLfloat v12[3], v23[3], v31[3];
  93    GLint i;
  94
  95    if (!depth)
  96    {
  97      drawTriangle(v1, v2, v3);
  98      return;
  99    }
 100    for (i = 0; i < 3; i++)
 101    {
 102        v12[i] = (v1[i] + v2[i]);
 103        v23[i] = (v2[i] + v3[i]);
 104        v31[i] = (v3[i] + v1[i]);
 105    }
 106    NormalizeVector(v12);
 107    NormalizeVector(v23);
 108    NormalizeVector(v31);
 109    subdivide(v1, v12, v31, depth - 1);
 110    subdivide(v2, v23, v12, depth - 1);
 111    subdivide(v3, v31, v23, depth - 1);
 112    subdivide(v12, v23, v31,depth - 1);
 113
 114 } */
 115 /*
 116 void display(void)
 117 {
 118    int i;
 119    xglClear(GL_COLOR_BUFFER_BIT);
 120    xglPushMatrix();
 121    xglRotatef(spin, 0.0, 0.0, 0.0);
 122    xglColor3f(1.0, 1.0, 0.0);
 123 //   xglBegin(GL_LINE_LOOP);
 124    for (i = 0; i < 20; i++)
 125    {
 126
 127        //xglVertex3fv(&vdata[tindices[i][0]][0]);
 128        //xglVertex3fv(&vdata[tindices[i][1]][0]);
 129        //xglVertex3fv(&vdata[tindices[i][2]][0]);
 130
 131        subdivide(&vdata[tindices[i][0]][0],
 132                  &vdata[tindices[i][1]][0],
 133                  &vdata[tindices[i][2]][0], 3);
 134
 135
 136    }
 137 //   xglEnd();
 138   // xglFlush();
 139   xglPopMatrix();
 140   glutSwapBuffers();
 141 } */
 142
 143 /* --------------------------------------------------------------
 144
 145       This section contains the code that calculates the actual
 146       position of the moon in the night sky.
 147
 148 ----------------------------------------------------------------*/
 149
 150 struct CelestialCoord fgCalculateMoon(struct OrbElements params,
 151                                       struct OrbElements sunParams,
 152                                       struct fgTIME t)
 153 {
 154   struct CelestialCoord
 155     geocCoord, topocCoord;
 156
 157
 158   double
 159     eccAnom, ecl, lonecl, latecl, actTime,
 160     xv, yv, v, r, xh, yh, zh, xg, yg, zg, xe, ye, ze,
 161     Ls, Lm, D, F, mpar, gclat, rho, HA, g;
 162
 163   struct fgAIRCRAFT *a;
 164   struct fgFLIGHT *f;
 165
 166   a = &current_aircraft;
 167   f = &a->flight;
 168
 169 /* calculate the angle between ecliptic and equatorial coordinate system */
 170   actTime = fgCalcActTime(t);
 171   ecl = fgDegToRad(23.4393 - 3.563E-7 * actTime);  // in radians of course
 172
 173   /* calculate the eccentric anomaly */
 174   eccAnom = fgCalcEccAnom(params.M, params.e);
 175
 176   /* calculate the moon's distance (d) and  true anomaly (v) */
 177   xv = params.a * ( cos(eccAnom) - params.e);
 178   yv = params.a * ( sqrt(1.0 - params.e*params.e) * sin(eccAnom));
 179   v =atan2(yv, xv);
 180   r = sqrt(xv*xv + yv*yv);
 181
 182   /* estimate the geocentric rectangular coordinates here */
 183   xh = r * (cos(params.N) * cos(v + params.w) - sin(params.N) * sin(v + params.w) * cos(params.i));
 184   yh = r * (sin(params.N) * cos(v + params.w) + cos(params.N) * sin(v + params.w) * cos(params.i));
 185   zh = r * (sin(v + params.w) * sin(params.i));
 186
 187   /* calculate the ecliptic latitude and longitude here */
 188   lonecl = atan2( yh, xh);
 189   latecl = atan2( zh, sqrt( xh*xh + yh*yh));
 190
 191   /* calculate a number of perturbations */
 192   Ls = sunParams.M + sunParams.w;
 193   Lm =    params.M +    params.w + params.N;
 194   D = Lm - Ls;
 195   F = Lm - params.N;
 196
 197   lonecl += fgDegToRad(
 198                        - 1.274 * sin (params.M - 2*D)                           // the Evection
 199                        + 0.658 * sin (2 * D)                                                    // the Variation
 200                        - 0.186 * sin (sunParams.M)                                      // the yearly variation
 201                        - 0.059 * sin (2*params.M - 2*D)
 202                        - 0.057 * sin (params.M - 2*D + sunParams.M)
 203                        + 0.053 * sin (params.M + 2*D)
 204                        + 0.046 * sin (2*D - sunParams.M)
 205                        + 0.041 * sin (params.M - sunParams.M)
 206                        - 0.035 * sin (D)                             // the Parallactic Equation
 207                        - 0.031 * sin (params.M + sunParams.M)
 208                        - 0.015 * sin (2*F - 2*D)
 209                        + 0.011 * sin (params.M - 4*D)
 210                        ); /* Pheeuuwwww */
 211   latecl += fgDegToRad(
 212                        - 0.173 * sin (F - 2*D)
 213                        - 0.055 * sin (params.M - F - 2*D)
 214                        - 0.046 * sin (params.M + F - 2*D)
 215                        + 0.033 * sin (F + 2*D)
 216                        + 0.017 * sin (2 * params.M + F)
 217                        );  /* Yep */
 218
 219   r += (
 220         - 0.58 * cos(params.M - 2*D)
 221         - 0.46 * cos(2*D)
 222         ); /* Ok! */
 223
 224   xg = r * cos(lonecl) * cos(latecl);
 225   yg = r * sin(lonecl) * cos(latecl);
 226   zg = r *               sin(latecl);
 227
 228   xe  = xg;
 229   ye = yg * cos(ecl) - zg * sin(ecl);
 230   ze = yg * sin(ecl) + zg * cos(ecl);
 231
 232
 233
 234
 235   geocCoord.RightAscension = atan2(ye, xe);
 236   geocCoord.Declination = atan2(ze, sqrt(xe*xe + ye*ye));
 237
 238   /* New since 25 december 1997 */
 239   /* Calculate the moon's topocentric position instead of it's geocentric! */
 240
 241   mpar = asin( 1 / r); /* calculate the moon's parrallax, i.e. the apparent size of the
 242                           (equatorial) radius of the Earth, as seen from the moon */
 243   gclat = FG_Latitude - 0.083358 * sin (2 * fgDegToRad( FG_Latitude));
 244   rho = 0.99883 + 0.00167 * cos(2 * fgDegToRad(FG_Latitude));
 245
 246   if (geocCoord.RightAscension < 0)
 247     geocCoord.RightAscension += (2*FG_PI);
 248
 249   HA = t.lst - (3.8197186 * geocCoord.RightAscension);
 250
 251   g = atan (tan(gclat) / cos( (HA / 3.8197186)));
 252
 253
 254
 255   topocCoord.RightAscension = geocCoord.RightAscension - mpar * rho * cos(gclat) * sin(HA) / cos(geocCoord.Declination);
 256   topocCoord.Declination = geocCoord.Declination - mpar * rho * sin(gclat) * sin(g - geocCoord.Declination) / sin(g);
 257   return topocCoord;
 258 }
 259
 260
 261 void fgMoonInit( void ) {
 262     struct fgLIGHT *l;
 263     static int dl_exists = 0;
 264
 265     fgPrintf( FG_ASTRO, FG_INFO, "Initializing the Moon\n");
 266
 267     l = &cur_light_params;
 268
 269     /* position the moon */
 270     fgSolarSystemUpdate(&(pltOrbElements[1]), cur_time_params);
 271     moonPos = fgCalculateMoon(pltOrbElements[1], pltOrbElements[0],
 272                               cur_time_params);
 273 #ifdef DEBUG
 274     fgPrintf( FG_ASTRO, FG_DEBUG,
 275               "Moon found at %f (ra), %f (dec)\n", moonPos.RightAscension,
 276               moonPos.Declination);
 277 #endif
 278
 279     xMoon = 60000.0 * cos(moonPos.RightAscension) * cos(moonPos.Declination);
 280     yMoon = 60000.0 * sin(moonPos.RightAscension) * cos(moonPos.Declination);
 281     zMoon = 60000.0 * sin(moonPos.Declination);
 282
 283     if ( !dl_exists ) {
 284         dl_exists = 1;
 285
 286         /* printf("First time through, creating moon display list\n"); */
 287
 288         moon = xglGenLists(1);
 289         xglNewList(moon, GL_COMPILE );
 290
 291         /* xglMaterialfv(GL_FRONT, GL_AMBIENT, l->scene_clear);
 292            xglMaterialfv(GL_FRONT, GL_DIFFUSE, moon_color); */
 293
 294
 295         glutSolidSphere(1.0, 10, 10);
 296
 297         xglEndList();
 298     }
 299 }
 300
 301
 302 /* Draw the moon */
 303 void fgMoonRender( void ) {
 304     struct fgLIGHT *l;
 305     GLfloat white[4] = { 1.0, 1.0, 1.0, 1.0 };
 306
 307     /* printf("Rendering moon\n"); */
 308
 309     l = &cur_light_params;
 310
 311     xglMaterialfv(GL_FRONT, GL_AMBIENT, l->sky_color );
 312     xglMaterialfv(GL_FRONT, GL_DIFFUSE, white);
 313
 314     xglPushMatrix();
 315     xglTranslatef(xMoon, yMoon, zMoon);
 316     xglScalef(1400, 1400, 1400);
 317
 318     xglCallList(moon);
 319
 320     xglPopMatrix();
 321 }
 322
 323
 324 /* $Log$
 325 /* Revision 1.4  1998/01/27 00:47:46  curt
 326 /* Incorporated Paul Bleisch's <bleisch@chromatic.com> new debug message
 327 /* system and commandline/config file processing code.
 328 /*
 329  * Revision 1.3  1998/01/19 19:26:57  curt
 330  * Merged in make system changes from Bob Kuehne <rpk@sgi.com>
 331  * This should simplify things tremendously.
 332  *
 333  * Revision 1.2  1998/01/19 18:40:16  curt
 334  * Tons of little changes to clean up the code and to remove fatal errors
 335  * when building with the c++ compiler.
 336  *
 337  * Revision 1.1  1998/01/07 03:16:16  curt
 338  * Moved from .../Src/Scenery/ to .../Src/Astro/
 339  *
 340  * Revision 1.16  1998/01/06 01:20:24  curt
 341  * Tweaks to help building with MSVC++
 342  *
 343  * Revision 1.15  1998/01/05 18:44:35  curt
 344  * Add an option to advance/decrease time from keyboard.
 345  *
 346  * Revision 1.14  1997/12/30 20:47:50  curt
 347  * Integrated new event manager with subsystem initializations.
 348  *
 349  * Revision 1.13  1997/12/30 16:41:00  curt
 350  * Added log at end of file.
 351  *
 352  */