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