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