Scenery/moon.c

   1 /**************************************************************************
   2  * moon.c
   3  *
   4  * This program is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU General Public License as
   6  * published by the Free Software Foundation; either version 2 of the
   7  * License, or (at your option) any later version.
   8  *
   9  * This program is distributed in the hope that it will be useful, but
  10  * WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  * General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU General Public License
  15  * along with this program; if not, write to the Free Software
  16  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  17  *
  18  * $Id$
  19  * (Log is kept at end of this file)
  20  **************************************************************************/
  21
  22
  23 #include <math.h>
  24 #include <GL/glut.h>
  25
  26 #include "orbits.h"
  27 #include "moon.h"
  28
  29 #include "../Main/views.h"
  30 #include "../Time/fg_time.h"
  31 /* #include "../Aircraft/aircraft.h"*/
  32 #include "../general.h"
  33
  34
  35 static GLint moon;
  36
  37 struct CelestialCoord fgCalculateMoon(struct OrbElements params,
  38                                       struct OrbElements sunParams,
  39                                       struct fgTIME t)
  40 {
  41         struct CelestialCoord
  42           result;
  43
  44     double
  45         eccAnom, ecl, lonecl, latecl, actTime,
  46         xv, yv, v, r, xh, yh, zh, xg, yg, zg, xe, ye, ze,
  47         Ls, Lm, D, F;
  48
  49     /* calculate the angle between ecliptic and equatorial coordinate system */
  50     actTime = fgCalcActTime(t);
  51     ecl = fgDegToRad(23.4393 - 3.563E-7 * actTime);  // in radians of course
  52
  53     /* calculate the eccentric anomaly */
  54     eccAnom = fgCalcEccAnom(params.M, params.e);
  55
  56     /* calculate the moon's distance (d) and  true anomaly (v) */
  57         xv = params.a * ( cos(eccAnom) - params.e);
  58     yv = params.a * ( sqrt(1.0 - params.e*params.e) * sin(eccAnom));
  59     v =atan2(yv, xv);
  60     r = sqrt(xv*xv + yv*yv);
  61
  62     /* estimate the geocentric rectangular coordinates here */
  63     xh = r * (cos(params.N) * cos(v + params.w) - sin(params.N) * sin(v + params.w) * cos(params.i));
  64     yh = r * (sin(params.N) * cos(v + params.w) + cos(params.N) * sin(v + params.w) * cos(params.i));
  65     zh = r * (sin(v + params.w) * sin(params.i));
  66
  67     /* calculate the ecliptic latitude and longitude here */
  68     lonecl = atan2( yh, xh);
  69     latecl = atan2( zh, sqrt( xh*xh + yh*yh));
  70
  71     /* calculate a number of perturbations */
  72     Ls = sunParams.M + sunParams.w;
  73     Lm =    params.M +    params.w + params.N;
  74     D = Lm - Ls;
  75     F = Lm - params.N;
  76
  77     lonecl += fgDegToRad(
  78                         - 1.274 * sin (params.M - 2*D)                          // the Evection
  79                 + 0.658 * sin (2 * D)                                                   // the Variation
  80                 - 0.186 * sin (sunParams.M)                                     // the yearly variation
  81                 - 0.059 * sin (2*params.M - 2*D)
  82                 - 0.057 * sin (params.M - 2*D + sunParams.M)
  83                 + 0.053 * sin (params.M + 2*D)
  84                 + 0.046 * sin (2*D - sunParams.M)
  85                 + 0.041 * sin (params.M - sunParams.M)
  86                 - 0.035 * sin (D)                             // the Parallactic Equation
  87                 - 0.031 * sin (params.M + sunParams.M)
  88                 - 0.015 * sin (2*F - 2*D)
  89                 + 0.011 * sin (params.M - 4*D)
  90               ); /* Pheeuuwwww */
  91     latecl += fgDegToRad(
  92                 - 0.173 * sin (F - 2*D)
  93                 - 0.055 * sin (params.M - F - 2*D)
  94                 - 0.046 * sin (params.M + F - 2*D)
  95                 + 0.033 * sin (F + 2*D)
  96                 + 0.017 * sin (2 * params.M + F)
  97               );  /* Yep */
  98
  99     r += (
 100                   - 0.58 * cos(params.M - 2*D)
 101               - 0.46 * cos(2*D)
 102          );
 103     xg = r * cos(lonecl) * cos(latecl);
 104     yg = r * sin(lonecl) * cos(latecl);
 105     zg = r *               sin(latecl);
 106
 107     xe  = xg;
 108     ye = yg * cos(ecl) - zg * sin(ecl);
 109     ze = yg * sin(ecl) + zg * cos(ecl);
 110
 111         result.RightAscension = atan2(ye, xe);
 112     result.Declination = atan2(ze, sqrt(xe*xe + ye*ye));
 113
 114     return result;
 115 }
 116
 117
 118 void fgMoonInit()
 119 {
 120    struct CelestialCoord
 121        moonPos;
 122
 123    moon = glGenLists(1);
 124    glNewList(moon, GL_COMPILE );
 125         glBegin( GL_POINTS );
 126    moonPos = fgCalculateMoon(pltOrbElements[1], pltOrbElements[0], cur_time_params);
 127    printf("Moon found at %f (ra), %f (dec)\n", moonPos.RightAscension, moonPos.Declination);
 128    /* give the moon a temporary color, for testing purposes */
 129    glColor3f( 0.0, 1.0, 0.0);
 130    glVertex3f( 190000.0 * cos(moonPos.RightAscension) * cos(moonPos.Declination),
 131                190000.0 * sin(moonPos.RightAscension) * cos(moonPos.Declination),
 132                  190000.0 * sin(moonPos.Declination) );
 133    glEnd();
 134    glEndList();
 135 }
 136
 137 void fgMoonRender()
 138 {
 139     double angle;
 140     static double warp = 0;
 141     struct VIEW *v;
 142     struct fgTIME *t;
 143
 144     t = &cur_time_params;
 145     v = &current_view;
 146
 147
 148    glDisable( GL_FOG );
 149         glDisable( GL_LIGHTING );
 150         glPushMatrix();
 151         glTranslatef( v->view_pos.x, v->view_pos.y, v->view_pos.z );
 152         angle = t->gst * 15.0;  /* 15 degrees per hour rotation */
 153         /* warp += 1.0; */
 154         /* warp = 15.0; */
 155         warp = 0.0;
 156         glRotatef( (angle+warp), 0.0, 0.0, -1.0 );
 157         printf("Rotating moon by %.2f degrees + %.2f degrees\n",angle,warp);
 158
 159         glCallList(moon);
 160
 161         glPopMatrix();
 162         glEnable( GL_LIGHTING );
 163         glEnable( GL_FOG );
 164 }
 165
 166
 167 /* $Log$
 168 /* Revision 1.2  1997/10/28 21:00:21  curt
 169 /* Changing to new terrain format.
 170 /*
 171  * Revision 1.1  1997/10/25 03:16:08  curt
 172  * Initial revision of code contributed by Durk Talsma.
 173  *
 174  */