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 #include <Time/fg_time.h>
  27 #include <Astro/orbits.h>
  28 #include <Astro/planets.h>
  29 #include <Astro/sun.h>
  30 #include <Main/fg_debug.h>
  31
  32 struct CelestialCoord fgCalculatePlanet(struct OrbElements planet,
  33                                         struct OrbElements theSun,
  34                                         struct fgTIME t, int idx)
  35 {
  36    struct CelestialCoord
  37                 result;
  38
  39     struct SunPos
  40         SolarPosition;
  41
  42         double
  43         eccAnom, r, v, ecl, actTime, R, s, ir, Nr, B, FV, ring_magn,
  44         xv, yv, xh, yh, zh, xg, yg, zg, xe, ye, ze;
  45
  46       actTime = fgCalcActTime(t);
  47       /* calculate the angle between ecliptic and equatorial coordinate system */
  48       ecl = fgDegToRad(23.4393 - 3.563E-7 * actTime);
  49
  50
  51     /* calculate the eccentric anomaly */
  52         eccAnom  = fgCalcEccAnom(planet.M, planet.e);
  53
  54     /* calculate the planets distance (r) and true anomaly (v) */
  55     xv = planet.a * (cos(eccAnom) - planet.e);
  56     yv = planet.a * (sqrt(1.0 - planet.e*planet.e) * sin(eccAnom));
  57     v = atan2(yv, xv);
  58     r = sqrt ( xv*xv + yv*yv);
  59
  60     /* calculate the planets position in 3-dimensional space */
  61     xh = r * ( cos(planet.N) * cos(v+planet.w) - sin(planet.N) * sin(v+planet.w) * cos(planet.i));
  62     yh = r * ( sin(planet.N) * cos(v+planet.w) + cos(planet.N) * sin(v+planet.w) * cos(planet.i));
  63     zh = r * ( sin(v+planet.w) * sin(planet.i));
  64
  65     /* calculate the ecleptic longitude and latitude */
  66
  67     /*
  68     lonecl = atan2(yh, xh);
  69     latecl = atan2(zh, sqrt ( xh*xh + yh*yh));
  70     */
  71     /* calculate the solar position */
  72
  73     SolarPosition = fgCalcSunPos(theSun);
  74     xg = xh + SolarPosition.xs;
  75     yg = yh + SolarPosition.ys;
  76     zg = zh;
  77
  78     xe = xg;
  79     ye = yg * cos(ecl) - zg * sin(ecl);
  80     ze = yg * sin(ecl) + zg * cos(ecl);
  81
  82     result.RightAscension = atan2(ye,xe);
  83     result.Declination = atan2(ze, sqrt(xe*xe + ye*ye));
  84
  85
  86     /* Let's calculate the brightness of the planet */
  87     R = sqrt ( xg*xg + yg*yg + zg*zg);
  88     s = SolarPosition.dist;
  89     FV = acos(  (r*r + R*R - s*s) / (2*r*R));
  90     FV  *= 57.29578;  /* convert radians to degrees */
  91     switch(idx)
  92       {
  93       case 2: /* mercury */
  94         result.magnitude = -0.36 + 5*log10( r*R ) + 0.027 * FV + 2.2E-13 * pow(FV, 6);
  95         break;
  96       case 3: /*venus */
  97         result.magnitude = -4.34 + 5*log10( r*R ) + 0.013 * FV + 4.2E-07 * pow(FV,3);
  98         break;
  99       case 4: /* mars */
 100         result.magnitude = -1.51 + 5*log10( r*R ) + 0.016 * FV;
 101         break;
 102       case 5: /* Jupiter */
 103         result.magnitude = -9.25 + 5*log10( r*R ) + 0.014 * FV;
 104         break;
 105       case 6: /* Saturn */
 106
 107         ir = 0.4897394;
 108         Nr = 2.9585076 + 6.6672E-7*actTime;
 109
 110         B = asin ( sin (result.Declination) * cos(ir) - cos(result.Declination) * sin (ir) * sin (result.RightAscension - Nr));
 111         ring_magn = -2.6 * sin (abs(B)) + 1.2 * pow(sin(B),2);
 112         result.magnitude = -9.0 + 5*log10( r*R ) + 0.044 * FV + ring_magn;
 113         break;
 114       case 7: /* Uranus */
 115         result.magnitude = -7.15 + 5*log10( r*R) + 0.001 * FV;
 116         break;
 117       case 8:  /* Neptune */
 118         result.magnitude = -6.90 + 5*log10 (r*R) + 0.001 *FV;
 119         break;
 120       default:
 121         fgPrintf( FG_ASTRO, FG_ALERT, "index %d out of range !!!!\n", idx);
 122       }
 123     fgPrintf( FG_ASTRO, FG_DEBUG,
 124               "    Planet found at %f (ra), %f (dec)\n",
 125               result.RightAscension, result.Declination);
 126     fgPrintf( FG_ASTRO, FG_DEBUG,
 127               "      Geocentric dist     %f\n"
 128               "      Heliocentric dist   %f\n"
 129               "      Distance to the sun %f\n"
 130               "      Phase angle         %f\n"
 131               "      Brightness          %f\n", R, r, s, FV, result.magnitude);
 132     return result;
 133 }
 134
 135
 136
 137 /* $Log$
 138 /* Revision 1.3  1998/01/27 00:47:47  curt
 139 /* Incorporated Paul Bleisch's <bleisch@chromatic.com> new debug message
 140 /* system and commandline/config file processing code.
 141 /*
 142  * Revision 1.2  1998/01/19 19:26:59  curt
 143  * Merged in make system changes from Bob Kuehne <rpk@sgi.com>
 144  * This should simplify things tremendously.
 145  *
 146  * Revision 1.1  1998/01/07 03:16:18  curt
 147  * Moved from .../Src/Scenery/ to .../Src/Astro/
 148  *
 149  * Revision 1.4  1997/12/30 20:47:52  curt
 150  * Integrated new event manager with subsystem initializations.
 151  *
 152  * Revision 1.3  1997/12/30 16:36:52  curt
 153  * Merged in Durk's changes ...
 154  *
 155  * Revision 1.2  1997/12/12 21:41:29  curt
 156  * More light/material property tweaking ... still a ways off.
 157  *
 158  * Revision 1.1  1997/10/25 03:16:10  curt
 159  * Initial revision of code contributed by Durk Talsma.
 160  *
 161  */
 162
 163