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     l = &cur_light_params;
 264
 265     /* position the moon */
 266     fgSolarSystemUpdate(&(pltOrbElements[1]), cur_time_params);
 267     moonPos = fgCalculateMoon(pltOrbElements[1], pltOrbElements[0],
 268                               cur_time_params);
 269 #ifdef DEBUG
 270     printf("Moon found at %f (ra), %f (dec)\n", moonPos.RightAscension,
 271            moonPos.Declination);
 272 #endif
 273
 274     if ( !dl_exists ) {
 275         dl_exists = 1;
 276
 277         /* printf("First time through, creating moon display list\n"); */
 278
 279         moon = xglGenLists(1);
 280         xglNewList(moon, GL_COMPILE );
 281
 282         /* xglMaterialfv(GL_FRONT, GL_AMBIENT, l->scene_clear);
 283            xglMaterialfv(GL_FRONT, GL_DIFFUSE, moon_color); */
 284
 285
 286         xMoon = 60000.0 * cos(moonPos.RightAscension) *
 287             cos(moonPos.Declination);
 288         yMoon = 60000.0 * sin(moonPos.RightAscension) *
 289             cos(moonPos.Declination);
 290         zMoon = 60000.0 * sin(moonPos.Declination);
 291
 292         glutSolidSphere(1.0, 10, 10);
 293
 294         xglEndList();
 295     }
 296 }
 297
 298
 299 /* Draw the moon */
 300 void fgMoonRender() {
 301     struct fgLIGHT *l;
 302     GLfloat white[4] = { 1.0, 1.0, 1.0, 1.0 };
 303
 304     l = &cur_light_params;
 305
 306     xglMaterialfv(GL_FRONT, GL_AMBIENT, l->sky_color );
 307     xglMaterialfv(GL_FRONT, GL_DIFFUSE, white);
 308
 309     xglPushMatrix();
 310     xglTranslatef(xMoon, yMoon, zMoon);
 311     xglScalef(1400, 1400, 1400);
 312
 313     xglCallList(moon);
 314
 315     xglPopMatrix();
 316 }
 317
 318
 319 /* $Log$
 320 /* Revision 1.13  1997/12/30 16:41:00  curt
 321 /* Added log at end of file.
 322 /*
 323  */