Astro/planets.cxx

   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 HAVE_CONFIG_H
  27 #  include <config.h>
  28 #endif
  29
  30 #ifdef HAVE_WINDOWS_H
  31 #  include <windows.h>
  32 #endif
  33
  34 #include <GL/glut.h>
  35 #include <XGL/xgl.h>
  36
  37 #include <Debug/fg_debug.h>
  38 #include <Include/fg_constants.h>
  39 #include <Time/fg_time.hxx>
  40 #include <Time/light.hxx>
  41
  42 #include "orbits.hxx"
  43 #include "planets.hxx"
  44 #include "sun.hxx"
  45
  46 GLint planets = 0;
  47
  48 void fgCalculatePlanet(struct OrbElements planet, struct OrbElements theSun,
  49                        fgTIME t, int idx, struct CelestialCoord *result)
  50 {
  51     // struct CelestialCoord result;
  52
  53     fgSUNPOS SolarPosition;
  54
  55     double eccAnom, r, v, ecl, actTime, R, s, ir, Nr, B, FV, ring_magn,
  56         xv, yv, xh, yh, zh, xg, yg, zg, xe, ye, ze;
  57
  58     actTime = fgCalcActTime(t);
  59     /*calculate the angle between ecliptic and equatorial coordinate system */
  60     ecl = DEG_TO_RAD * (23.4393 - 3.563E-7 * actTime);
  61
  62     /* calculate the eccentric anomaly */
  63     eccAnom  = fgCalcEccAnom(planet.M, planet.e);
  64
  65     /* calculate the planets distance (r) and true anomaly (v) */
  66     xv = planet.a * (cos(eccAnom) - planet.e);
  67     yv = planet.a * (sqrt(1.0 - planet.e*planet.e) * sin(eccAnom));
  68     v = atan2(yv, xv);
  69     r = sqrt ( xv*xv + yv*yv);
  70
  71     /* calculate the planets position in 3-dimensional space */
  72     xh = r * (cos (planet.N) * cos (v + planet.w) -
  73               sin (planet.N) * sin (v + planet.w) * cos (planet.i));
  74     yh = r * (sin (planet.N) * cos (v + planet.w) +
  75               cos (planet.N) * sin (v + planet.w) * cos (planet.i));
  76     zh = r * ( sin(v+planet.w) * sin(planet.i));
  77
  78     /* calculate the ecliptic longitude and latitude */
  79
  80     xg = xh + solarPosition.xs;
  81     yg = yh + solarPosition.ys;
  82     zg = zh;
  83
  84     xe = xg;
  85     ye = yg * cos(ecl) - zg * sin(ecl);
  86     ze = yg * sin(ecl) + zg * cos(ecl);
  87
  88     result->RightAscension = atan2(ye,xe);
  89     result->Declination = atan2(ze, sqrt(xe*xe + ye*ye));
  90
  91     /* Let's calculate the brightness of the planet */
  92     R = sqrt ( xg*xg + yg*yg + zg*zg);
  93     s = SolarPosition.dist;
  94     FV = acos(  (r*r + R*R - s*s) / (2*r*R));
  95     FV  *= 57.29578;  /* convert radians to degrees */
  96     switch(idx)
  97       {
  98       case 2: /* mercury */
  99         result->magnitude = -0.36 + 5*log10( r*R ) + 0.027 * FV + 2.2E-13 * pow(FV, 6);
 100         break;
 101       case 3: /*venus */
 102         result->magnitude = -4.34 + 5*log10( r*R ) + 0.013 * FV + 4.2E-07 * pow(FV,3);
 103         break;
 104       case 4: /* mars */
 105         result->magnitude = -1.51 + 5*log10( r*R ) + 0.016 * FV;
 106         break;
 107       case 5: /* Jupiter */
 108         result->magnitude = -9.25 + 5*log10( r*R ) + 0.014 * FV;
 109         break;
 110       case 6: /* Saturn */
 111
 112         ir = 0.4897394;
 113         Nr = 2.9585076 + 6.6672E-7*actTime;
 114
 115         B = asin (sin (result->Declination) * cos (ir) -
 116                   cos (result->Declination) * sin (ir) *
 117                   sin (result->RightAscension - Nr));
 118         ring_magn = -2.6 * sin (fabs(B)) + 1.2 * pow(sin(B),2);
 119         result->magnitude = -9.0 + 5*log10( r*R ) + 0.044 * FV + ring_magn;
 120         break;
 121       case 7: /* Uranus */
 122         result->magnitude = -7.15 + 5*log10( r*R) + 0.001 * FV;
 123         break;
 124       case 8:  /* Neptune */
 125         result->magnitude = -6.90 + 5*log10 (r*R) + 0.001 *FV;
 126         break;
 127       default:
 128         fgPrintf( FG_ASTRO, FG_ALERT, "index %d out of range !!!!\n", idx);
 129       }
 130     fgPrintf( FG_ASTRO, FG_DEBUG,
 131               "    Planet found at %f (ra), %f (dec)\n",
 132               result->RightAscension, result->Declination);
 133     fgPrintf( FG_ASTRO, FG_DEBUG,
 134               "      Geocentric dist     %f\n"
 135               "      Heliocentric dist   %f\n"
 136               "      Distance to the sun %f\n"
 137               "      Phase angle         %f\n"
 138               "      Brightness          %f\n", R, r, s, FV, result->magnitude);
 139
 140     // return result;
 141 }
 142
 143
 144 void fgPlanetsInit( void )
 145 {
 146   fgLIGHT *l;
 147   int i;
 148   struct CelestialCoord pltPos;
 149   double magnitude;
 150
 151   l = &cur_light_params;
 152   /* if the display list was already built during a previous init,
 153      then recycle it */
 154
 155   if (planets) {
 156       xglDeleteLists(planets, 1);
 157   }
 158
 159   planets = xglGenLists(1);
 160   xglNewList( planets, GL_COMPILE );
 161   xglBegin( GL_POINTS );
 162
 163
 164   /* Add the planets to all four display lists */
 165   for ( i = 2; i < 9; i++ ) {
 166     fgSolarSystemUpdate(&(pltOrbElements[i]), cur_time_params);
 167     fgCalculatePlanet(pltOrbElements[i], pltOrbElements[0],
 168                       cur_time_params, i, &pltPos);
 169
 170     magnitude = (0.0 - pltPos.magnitude) / 5.0 + 1.0;
 171
 172     /* scale magnitudes again so they look ok */
 173     /* magnitude =
 174        magnitude * 0.7 + (((FG_STAR_LEVELS - 1) - i) * 0.1); */
 175
 176     /* the following statement could be made a little more sopisticated
 177        for the moment: Stick to this one */
 178     magnitude = magnitude * 0.7 + (3 * 0.1);
 179     if ( magnitude > 1.0 ) { magnitude = 1.0; }
 180     if ( magnitude < 0.0 ) { magnitude = 0.0; }
 181
 182     /* Add planets to the display list, based on sun_angle and current
 183        magnitude It's pretty experimental... */
 184
 185     if ((double) (l->sun_angle - FG_PI_2) >
 186         ((magnitude - 1.0) * -20 * DEG_TO_RAD))
 187         {
 188             xglColor3f (magnitude, magnitude, magnitude);
 189             printf ("Sun Angle to Horizon (in Rads) = %f\n",
 190                     (double) (l->sun_angle - FG_PI_2));
 191             printf ("Transformed Magnitude is :%f  %f\n",
 192                     magnitude, (magnitude - 1.0) * -20 * DEG_TO_RAD);
 193
 194             xglVertex3f (50000.0 * cos (pltPos.RightAscension) *
 195                          cos (pltPos.Declination),
 196                          50000.0 * sin (pltPos.RightAscension) *
 197                          cos (pltPos.Declination),
 198                          50000.0 * sin (pltPos.Declination));
 199         }
 200   }
 201   xglEnd();
 202   xglEndList();
 203
 204 }
 205
 206
 207 void fgPlanetsRender( void ) {
 208         xglCallList(planets);
 209 }
 210
 211
 212 /* $Log$
 213 /* Revision 1.6  1998/06/27 16:51:54  curt
 214 /* In fgCalculatePlanet() pass a pointer to a structure to be modified, rather
 215 /* than returning the entire structure.
 216 /*
 217  * Revision 1.5  1998/04/28 01:19:01  curt
 218  * Type-ified fgTIME and fgVIEW
 219  *
 220  * Revision 1.4  1998/04/26 05:10:01  curt
 221  * "struct fgLIGHT" -> "fgLIGHT" because fgLIGHT is typedef'd.
 222  *
 223  * Revision 1.3  1998/04/25 22:06:25  curt
 224  * Edited cvs log messages in source files ... bad bad bad!
 225  *
 226  * Revision 1.2  1998/04/24 00:45:02  curt
 227  * Wrapped "#include <config.h>" in "#ifdef HAVE_CONFIG_H"
 228  * Fixed a bug when generating sky colors.
 229  *
 230  * Revision 1.1  1998/04/22 13:21:31  curt
 231  * C++ - ifing the code a bit.
 232  *
 233  * Revision 1.9  1998/04/18 04:13:57  curt
 234  * Moved fg_debug.c to it's own library.
 235  *
 236  * Revision 1.8  1998/04/03 21:52:50  curt
 237  * Converting to Gnu autoconf system.
 238  *
 239  * Revision 1.7  1998/02/23 19:07:55  curt
 240  * Incorporated Durk's Astro/ tweaks.  Includes unifying the sun position
 241  * calculation code between sun display, and other FG sections that use this
 242  * for things like lighting.
 243  *
 244  * Revision 1.6  1998/02/12 21:59:36  curt
 245  * Incorporated code changes contributed by Charlie Hotchkiss
 246  * <chotchkiss@namg.us.anritsu.com>
 247  *
 248  * Revision 1.5  1998/02/03 23:20:12  curt
 249  * Lots of little tweaks to fix various consistency problems discovered by
 250  * Solaris' CC.  Fixed a bug in fg_debug.c with how the fgPrintf() wrapper
 251  * passed arguments along to the real printf().  Also incorporated HUD changes
 252  * by Michele America.
 253  *
 254  * Revision 1.4  1998/02/02 20:53:23  curt
 255  * To version 0.29
 256  *
 257  * Revision 1.3  1998/01/27 00:47:47  curt
 258  * Incorporated Paul Bleisch's <pbleisch@acm.org> new debug message
 259  * system and commandline/config file processing code.
 260  *
 261  * Revision 1.2  1998/01/19 19:26:59  curt
 262  * Merged in make system changes from Bob Kuehne <rpk@sgi.com>
 263  * This should simplify things tremendously.
 264  *
 265  * Revision 1.1  1998/01/07 03:16:18  curt
 266  * Moved from .../Src/Scenery/ to .../Src/Astro/
 267  *
 268  * Revision 1.4  1997/12/30 20:47:52  curt
 269  * Integrated new event manager with subsystem initializations.
 270  *
 271  * Revision 1.3  1997/12/30 16:36:52  curt
 272  * Merged in Durk's changes ...
 273  *
 274  * Revision 1.2  1997/12/12 21:41:29  curt
 275  * More light/material property tweaking ... still a ways off.
 276  *
 277  * Revision 1.1  1997/10/25 03:16:10  curt
 278  * Initial revision of code contributed by Durk Talsma.
 279  *
 280  */
 281
 282