Astro/planets.c

   1 /**************************************************************************
   2  * planets.c
   3  *
   4  * Written 1997 by Durk Talsma, started October, 1997.  For the flight gear
   5  * project.
   6  *
   7  * This program is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License as
   9  * published by the Free Software Foundation; either version 2 of the
  10  * License, or (at your option) any later version.
  11  *
  12  * This program is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15  * General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with this program; if not, write to the Free Software
  19  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20  *
  21  * $Id$
  22  * (Log is kept at end of this file)
  23  **************************************************************************/
  24
  25
  26 #ifdef WIN32
  27 #  include <windows.h>
  28 #endif
  29
  30 #include <GL/glut.h>
  31 #include <XGL/xgl.h>
  32
  33 #include <Time/fg_time.h>
  34 #include <Astro/orbits.h>
  35 #include <Astro/planets.h>
  36 #include <Astro/sun.h>
  37 #include <Include/fg_constants.h>
  38 #include <Main/fg_debug.h>
  39
  40 GLint planets = 0;
  41
  42 struct CelestialCoord fgCalculatePlanet(struct OrbElements planet,
  43                                         struct OrbElements theSun,
  44                                         struct fgTIME t, int idx)
  45 {
  46    struct CelestialCoord
  47                 result;
  48
  49     struct SunPos
  50         SolarPosition;
  51
  52         double
  53         eccAnom, r, v, ecl, actTime, R, s, ir, Nr, B, FV, ring_magn,
  54         xv, yv, xh, yh, zh, xg, yg, zg, xe, ye, ze;
  55
  56       actTime = fgCalcActTime(t);
  57       /* calculate the angle between ecliptic and equatorial coordinate system */
  58       /* ecl = DEG_TO_RAD * (23.4393 - 3.563E-7 * actTime); */
  59       ecl = 0.409093 - 6.2186E-9 * actTime;
  60
  61     /* calculate the eccentric anomaly */
  62         eccAnom  = fgCalcEccAnom(planet.M, planet.e);
  63
  64     /* calculate the planets distance (r) and true anomaly (v) */
  65     xv = planet.a * (cos(eccAnom) - planet.e);
  66     yv = planet.a * (sqrt(1.0 - planet.e*planet.e) * sin(eccAnom));
  67     v = atan2(yv, xv);
  68     r = sqrt ( xv*xv + yv*yv);
  69
  70     /* calculate the planets position in 3-dimensional space */
  71     xh = r * ( cos(planet.N) * cos(v+planet.w) - sin(planet.N) * sin(v+planet.w) * cos(planet.i));
  72     yh = r * ( sin(planet.N) * cos(v+planet.w) + cos(planet.N) * sin(v+planet.w) * cos(planet.i));
  73     zh = r * ( sin(v+planet.w) * sin(planet.i));
  74
  75     /* calculate the ecliptic longitude and latitude */
  76
  77     /*
  78     lonecl = atan2(yh, xh);
  79     latecl = atan2(zh, sqrt ( xh*xh + yh*yh));
  80     */
  81     /* calculate the solar position */
  82
  83     SolarPosition = fgCalcSunPos(theSun);
  84     xg = xh + SolarPosition.xs;
  85     yg = yh + SolarPosition.ys;
  86     zg = zh;
  87
  88     xe = xg;
  89     ye = yg * cos(ecl) - zg * sin(ecl);
  90     ze = yg * sin(ecl) + zg * cos(ecl);
  91
  92     result.RightAscension = atan2(ye,xe);
  93     result.Declination = atan2(ze, sqrt(xe*xe + ye*ye));
  94
  95
  96     /* Let's calculate the brightness of the planet */
  97     R = sqrt ( xg*xg + yg*yg + zg*zg);
  98     s = SolarPosition.dist;
  99     FV = acos(  (r*r + R*R - s*s) / (2*r*R));
 100     FV  *= 57.29578;  /* convert radians to degrees */
 101     switch(idx)
 102       {
 103       case 2: /* mercury */
 104         result.magnitude = -0.36 + 5*log10( r*R ) + 0.027 * FV + 2.2E-13 * pow(FV, 6);
 105         break;
 106       case 3: /*venus */
 107         result.magnitude = -4.34 + 5*log10( r*R ) + 0.013 * FV + 4.2E-07 * pow(FV,3);
 108         break;
 109       case 4: /* mars */
 110         result.magnitude = -1.51 + 5*log10( r*R ) + 0.016 * FV;
 111         break;
 112       case 5: /* Jupiter */
 113         result.magnitude = -9.25 + 5*log10( r*R ) + 0.014 * FV;
 114         break;
 115       case 6: /* Saturn */
 116
 117         ir = 0.4897394;
 118         Nr = 2.9585076 + 6.6672E-7*actTime;
 119
 120         B = asin ( sin (result.Declination) * cos(ir) - cos(result.Declination) * sin (ir) * sin (result.RightAscension - Nr));
 121         ring_magn = -2.6 * sin (fabs(B)) + 1.2 * pow(sin(B),2);
 122         result.magnitude = -9.0 + 5*log10( r*R ) + 0.044 * FV + ring_magn;
 123         break;
 124       case 7: /* Uranus */
 125         result.magnitude = -7.15 + 5*log10( r*R) + 0.001 * FV;
 126         break;
 127       case 8:  /* Neptune */
 128         result.magnitude = -6.90 + 5*log10 (r*R) + 0.001 *FV;
 129         break;
 130       default:
 131         fgPrintf( FG_ASTRO, FG_ALERT, "index %d out of range !!!!\n", idx);
 132       }
 133     fgPrintf( FG_ASTRO, FG_DEBUG,
 134               "    Planet found at %f (ra), %f (dec)\n",
 135               result.RightAscension, result.Declination);
 136     fgPrintf( FG_ASTRO, FG_DEBUG,
 137               "      Geocentric dist     %f\n"
 138               "      Heliocentric dist   %f\n"
 139               "      Distance to the sun %f\n"
 140               "      Phase angle         %f\n"
 141               "      Brightness          %f\n", R, r, s, FV, result.magnitude);
 142     return result;
 143 }
 144
 145
 146 void fgPlanetsInit( void )
 147 {
 148   int i;
 149   struct CelestialCoord pltPos;
 150   double magnitude;
 151
 152   /* if the display list was already built during a previous init,
 153      recycle it */
 154
 155   if (planets)
 156       xglDeleteLists(planets, 1);
 157
 158   planets = xglGenLists(1);
 159   xglNewList( planets, GL_COMPILE );
 160   xglBegin( GL_POINTS );
 161
 162
 163   /* Add the planets to all four display lists */
 164   for ( i = 2; i < 9; i++ ) {
 165     fgSolarSystemUpdate(&(pltOrbElements[i]), cur_time_params);
 166     pltPos = fgCalculatePlanet(pltOrbElements[i],
 167                                pltOrbElements[0], cur_time_params, i);
 168
 169     /* give the planets a temporary color, for testing purposes */
 170     /* xglColor3f( 1.0, 0.0, 0.0); */
 171     /* scale magnitudes to (0.0 - 1.0) */
 172
 173     magnitude = (0.0 - pltPos.magnitude) / 5.0 + 1.0;
 174
 175     /* scale magnitudes again so they look ok */
 176     /* magnitude =
 177        magnitude * 0.7 + (((FG_STAR_LEVELS - 1) - i) * 0.1); */
 178
 179     /* the following statement could be made a little more sopisticated
 180        for the moment: Stick to this one */
 181     magnitude = magnitude * 0.7 + (3 * 0.1);
 182     if ( magnitude > 1.0 ) { magnitude = 1.0; }
 183     if ( magnitude < 0.0 ) { magnitude = 0.0; }
 184
 185
 186     xglColor3f(magnitude, magnitude, magnitude);
 187     /* xglColor3f(1.0, 1.0,1.0); */
 188
 189     xglVertex3f( 50000.0 * cos(pltPos.RightAscension) *
 190                            cos(pltPos.Declination),
 191                  50000.0 * sin(pltPos.RightAscension) *
 192                            cos(pltPos.Declination),
 193                  50000.0 * sin(pltPos.Declination) );
 194   }
 195   xglEnd();
 196   xglEndList();
 197
 198 }
 199
 200
 201 void fgPlanetsRender( void ) {
 202         xglCallList(planets);
 203 }
 204
 205
 206 /* $Log$
 207 /* Revision 1.6  1998/02/12 21:59:36  curt
 208 /* Incorporated code changes contributed by Charlie Hotchkiss
 209 /* <chotchkiss@namg.us.anritsu.com>
 210 /*
 211  * Revision 1.5  1998/02/03 23:20:12  curt
 212  * Lots of little tweaks to fix various consistency problems discovered by
 213  * Solaris' CC.  Fixed a bug in fg_debug.c with how the fgPrintf() wrapper
 214  * passed arguments along to the real printf().  Also incorporated HUD changes
 215  * by Michele America.
 216  *
 217  * Revision 1.4  1998/02/02 20:53:23  curt
 218  * To version 0.29
 219  *
 220  * Revision 1.3  1998/01/27 00:47:47  curt
 221  * Incorporated Paul Bleisch's <bleisch@chromatic.com> new debug message
 222  * system and commandline/config file processing code.
 223  *
 224  * Revision 1.2  1998/01/19 19:26:59  curt
 225  * Merged in make system changes from Bob Kuehne <rpk@sgi.com>
 226  * This should simplify things tremendously.
 227  *
 228  * Revision 1.1  1998/01/07 03:16:18  curt
 229  * Moved from .../Src/Scenery/ to .../Src/Astro/
 230  *
 231  * Revision 1.4  1997/12/30 20:47:52  curt
 232  * Integrated new event manager with subsystem initializations.
 233  *
 234  * Revision 1.3  1997/12/30 16:36:52  curt
 235  * Merged in Durk's changes ...
 236  *
 237  * Revision 1.2  1997/12/12 21:41:29  curt
 238  * More light/material property tweaking ... still a ways off.
 239  *
 240  * Revision 1.1  1997/10/25 03:16:10  curt
 241  * Initial revision of code contributed by Durk Talsma.
 242  *
 243  */
 244
 245