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
  27 #include "orbits.h"
  28 #include "moon.h"
  29
  30 #include "../general.h"
  31 #include "../Main/views.h"
  32 #include "../Time/fg_time.h"
  33
  34 struct CelestialCoord
  35     moonPos;
  36
  37 float xMoon, yMoon, zMoon;
  38
  39 /*
  40 static GLfloat vdata[12][3] =
  41 {
  42    {-X, 0.0, Z }, { X, 0.0, Z }, {-X, 0.0, -Z}, {X, 0.0, -Z },
  43    { 0.0, Z, X }, { 0.0, Z, -X}, {0.0, -Z, -X}, {0.0, -Z, -X},
  44    { Z, X, 0.0 }, { -Z, X, 0.0}, {Z, -X, 0.0 }, {-Z, -X, 0.0}
  45 };
  46
  47 static GLuint tindices[20][3] =
  48 {
  49    {0,4,1}, {0,9,4}, {9,5,4}, {4,5,8}, {4,8,1},
  50    {8,10,1}, {8,3,10}, {5,3,8}, {5,2,3}, {2,7,3},
  51    {7,10,3}, {7,6,10}, {7,11,6}, {11,0,6}, {0,1,6},
  52    {6,1,10}, {9,0,11}, {9,11,2}, {9,2,5}, {7,2,11}
  53 };*/
  54
  55 GLint moon;
  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    glBegin(GL_POLYGON);
  78    //glBegin(GL_POINTS);
  79       glNormal3fv(v1);
  80       glVertex3fv(v1);
  81       glNormal3fv(v2);
  82       glVertex3fv(v2);
  83       glNormal3fv(v3);
  84       glVertex3fv(v3);
  85    glEnd();
  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    glClear(GL_COLOR_BUFFER_BIT);
 118    glPushMatrix();
 119    glRotatef(spin, 0.0, 0.0, 0.0);
 120    glColor3f(1.0, 1.0, 0.0);
 121 //   glBegin(GL_LINE_LOOP);
 122    for (i = 0; i < 20; i++)
 123    {
 124
 125        //glVertex3fv(&vdata[tindices[i][0]][0]);
 126        //glVertex3fv(&vdata[tindices[i][1]][0]);
 127        //glVertex3fv(&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 //   glEnd();
 136   // glFlush();
 137   glPopMatrix();
 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           result;
 154
 155     double
 156           eccAnom, ecl, lonecl, latecl, actTime,
 157         xv, yv, v, r, xh, yh, zh, xg, yg, zg, xe, ye, ze,
 158         Ls, Lm, D, F;
 159
 160     /* calculate the angle between ecliptic and equatorial coordinate system */
 161     actTime = fgCalcActTime(t);
 162     ecl = fgDegToRad(23.4393 - 3.563E-7 * actTime);  // in radians of course
 163
 164     /* calculate the eccentric anomaly */
 165     eccAnom = fgCalcEccAnom(params.M, params.e);
 166
 167     /* calculate the moon's distance (d) and  true anomaly (v) */
 168          xv = params.a * ( cos(eccAnom) - params.e);
 169     yv = params.a * ( sqrt(1.0 - params.e*params.e) * sin(eccAnom));
 170     v =atan2(yv, xv);
 171     r = sqrt(xv*xv + yv*yv);
 172
 173     /* estimate the geocentric rectangular coordinates here */
 174     xh = r * (cos(params.N) * cos(v + params.w) - sin(params.N) * sin(v + params.w) * cos(params.i));
 175     yh = r * (sin(params.N) * cos(v + params.w) + cos(params.N) * sin(v + params.w) * cos(params.i));
 176     zh = r * (sin(v + params.w) * sin(params.i));
 177
 178     /* calculate the ecliptic latitude and longitude here */
 179     lonecl = atan2( yh, xh);
 180     latecl = atan2( zh, sqrt( xh*xh + yh*yh));
 181
 182     /* calculate a number of perturbations */
 183     Ls = sunParams.M + sunParams.w;
 184     Lm =    params.M +    params.w + params.N;
 185     D = Lm - Ls;
 186     F = Lm - params.N;
 187
 188     lonecl += fgDegToRad(
 189                             - 1.274 * sin (params.M - 2*D)                      // the Evection
 190                 + 0.658 * sin (2 * D)                                                   // the Variation
 191                 - 0.186 * sin (sunParams.M)                                     // the yearly variation
 192                 - 0.059 * sin (2*params.M - 2*D)
 193                 - 0.057 * sin (params.M - 2*D + sunParams.M)
 194                 + 0.053 * sin (params.M + 2*D)
 195                 + 0.046 * sin (2*D - sunParams.M)
 196                 + 0.041 * sin (params.M - sunParams.M)
 197                 - 0.035 * sin (D)                             // the Parallactic Equation
 198                 - 0.031 * sin (params.M + sunParams.M)
 199                 - 0.015 * sin (2*F - 2*D)
 200                 + 0.011 * sin (params.M - 4*D)
 201               ); /* Pheeuuwwww */
 202     latecl += fgDegToRad(
 203                 - 0.173 * sin (F - 2*D)
 204                 - 0.055 * sin (params.M - F - 2*D)
 205                 - 0.046 * sin (params.M + F - 2*D)
 206                 + 0.033 * sin (F + 2*D)
 207                 + 0.017 * sin (2 * params.M + F)
 208               );  /* Yep */
 209
 210     r += (
 211                      - 0.58 * cos(params.M - 2*D)
 212               - 0.46 * cos(2*D)
 213           ); /* Ok! */
 214
 215     xg = r * cos(lonecl) * cos(latecl);
 216     yg = r * sin(lonecl) * cos(latecl);
 217     zg = r *               sin(latecl);
 218
 219     xe  = xg;
 220     ye = yg * cos(ecl) - zg * sin(ecl);
 221     ze = yg * sin(ecl) + zg * cos(ecl);
 222
 223          result.RightAscension = atan2(ye, xe);
 224     result.Declination = atan2(ze, sqrt(xe*xe + ye*ye));
 225     return result;
 226 }
 227
 228
 229 void fgMoonInit()
 230 {
 231 //   int i;
 232 //   moon = glGenLists(1);
 233 //   glNewList(moon, GL_COMPILE );
 234
 235    fgSolarSystemUpdate(&(pltOrbElements[1]), cur_time_params);
 236    moonPos = fgCalculateMoon(pltOrbElements[1], pltOrbElements[0], cur_time_params);
 237    #ifdef DEBUG
 238    printf("Moon found at %f (ra), %f (dec)\n", moonPos.RightAscension, moonPos.Declination);
 239    #endif
 240    glColor3f(1.0, 1.0, 0.0);
 241
 242    /* xMoon = 90000.0 * cos(moonPos.RightAscension) * cos(moonPos.Declination);
 243    yMoon = 90000.0 * sin(moonPos.RightAscension) * cos(moonPos.Declination);
 244    zMoon = 90000.0 * sin(moonPos.Declination); */
 245
 246    xMoon = 60000.0 * cos(moonPos.RightAscension) * cos(moonPos.Declination);
 247    yMoon = 60000.0 * sin(moonPos.RightAscension) * cos(moonPos.Declination);
 248    zMoon = 60000.0 * sin(moonPos.Declination);
 249
 250 //   glPushMatrix();
 251 //   glBegin(GL_TRIANGLES);
 252    /*
 253    for (i = 0; i < 20; i++)
 254       subdivide(&vdata[tindices[i][0]][0],
 255                 &vdata[tindices[i][1]][0],
 256                 &vdata[tindices[i][2]][0], 3);*/
 257 //     glutSolidSphere(1.0, 25, 25);
 258
 259 //   glEnd();
 260 //   glPopMatrix();
 261 //   glEndList();
 262 }
 263
 264
 265 /* Draw the moon */
 266 void fgMoonRender() {
 267     GLfloat color[4] = { 1.0, 1.0, 1.0, 1.0 };
 268
 269     /* set lighting parameters */
 270     glLightfv(GL_LIGHT0, GL_AMBIENT, color );
 271     glLightfv(GL_LIGHT0, GL_DIFFUSE, color );
 272     glMaterialfv(GL_FRONT, GL_AMBIENT, fgClearColor);
 273     glMaterialfv(GL_FRONT, GL_DIFFUSE, color);
 274
 275     glPushMatrix();
 276     glTranslatef(xMoon, yMoon, zMoon);
 277     glScalef(1400, 1400, 1400);
 278     /* glutSolidSphere(1.0, 25, 25); */
 279     glutSolidSphere(1.0, 15, 15);
 280     glPopMatrix();
 281 }
 282