Astro/orbits.c

   1 /**************************************************************************
   2  * orbits.c - calculates the orbital elements of the sun, moon and planets.
   3  *            For inclusion in flight gear
   4  *
   5  * Written 1997 by Durk Talsma, started October 19, 1997.
   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 <math.h>
  27 #include <string.h>
  28
  29 #include <Astro/orbits.h>
  30
  31 #include <Include/fg_constants.h>
  32 #include <Include/general.h>
  33 #include <Time/fg_time.h>
  34 #include <Main/fg_debug.h>
  35
  36 struct OrbElements pltOrbElements[9];
  37
  38
  39 //double fgCalcActTime(struct fgTIME t)
  40 //{
  41 //   double
  42 //         actTime, UT;
  43 //   int year;
  44 //
  45 //   /* a hack. This one introduces the 2000 problem into the program */
  46 //   year = t.gmt->tm_year + 1900;
  47 //
  48 //   /* calculate the actual time, remember to add 1 to tm_mon! */
  49 //   actTime = 367 * year - 7 *
  50 //                (year + ((t.gmt->tm_mon+1) + 9) / 12) / 4 + 275 *
  51 //                 (t.gmt->tm_mon+1) / 9 + t.gmt->tm_mday - 730530;
  52 //
  53 //    UT = t.gmt->tm_hour + ((double) t.gmt->tm_min / 60);
  54 //    /*printf("UT = %f\n", UT); */
  55 //    actTime += (UT / 24.0);
  56 //    #define DEBUG 1
  57 //    #ifdef DEBUG
  58 //    /* printf("  Actual Time:\n"); */
  59 //    /* printf("  current day = %f\t", actTime); */
  60 //    /* printf("  GMT = %d, %d, %d, %d, %d, %d\n",
  61 //         year, t.gmt->tm_mon, t.gmt->tm_mday,
  62 //         t.gmt->tm_hour, t.gmt->tm_min, t.gmt->tm_sec); */
  63 //    #endif
  64 //    return actTime;
  65 //}
  66
  67
  68 double fgCalcActTime(struct fgTIME t)
  69 {
  70   return (t.mjd - 36523.5);
  71 }
  72
  73
  74 /* convert degrees to radians */
  75 /*
  76 double fgDegToRad(double angle)
  77 {
  78         return (angle * PIOVER180);
  79 }
  80 */
  81 double fgCalcEccAnom(double M, double e)
  82 {
  83     double
  84         eccAnom, E0, E1, diff;
  85
  86     eccAnom = M + e * sin(M) * (1.0 + e * cos(M));
  87     /* iterate to achieve a greater precision for larger eccentricities */
  88     if (e > 0.05)
  89     {
  90         E0 = eccAnom;
  91         do
  92         {
  93                 E1 = E0 - (E0 - e * sin(E0) - M) / (1 - e * cos(E0));
  94             diff = fabs(E0 - E1);
  95             E0 = E1;
  96                 }
  97         while (diff > (DEG_TO_RAD * 0.001));
  98         return E0;
  99         }
 100     return eccAnom;
 101 }
 102
 103
 104
 105 void fgReadOrbElements(struct OrbElements *dest, FILE *src)
 106 {
 107         char line[256];
 108     int i,j;
 109     j = 0;
 110     do
 111     {
 112         fgets(line, 256,src);
 113         for (i = 0; i < 256; i++)
 114         {
 115                 if (line[i] == '#')
 116                 line[i] = 0;
 117         }
 118         /*printf("Reading line %d\n", j++); */
 119     }
 120     while (!(strlen(line)));
 121     sscanf(line, "%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf\n",
 122                 &dest->NFirst, &dest->NSec,
 123             &dest->iFirst, &dest->iSec,
 124             &dest->wFirst, &dest->wSec,
 125             &dest->aFirst, &dest->aSec,
 126             &dest->eFirst, &dest->eSec,
 127             &dest->MFirst, &dest->MSec);
 128 }
 129
 130
 131 void fgSolarSystemInit(struct fgTIME t)
 132 {
 133    struct fgGENERAL *g;
 134    char path[80];
 135    int i;
 136    FILE *data;
 137
 138    fgPrintf( FG_ASTRO, FG_INFO, "Initializing solar system\n");
 139
 140    /* build the full path name to the orbital elements database file */
 141    g = &general;
 142    path[0] = '\0';
 143    strcat(path, g->root_dir);
 144    strcat(path, "/Scenery/");
 145    strcat(path, "Planets.dat");
 146
 147    if ( (data = fopen(path, "r")) == NULL )
 148    {
 149             fgPrintf( FG_ASTRO, FG_ALERT,
 150                       "Cannot open data file: '%s'\n", path);
 151             return;
 152    }
 153    /* printf("  reading datafile %s\n", path); */
 154    fgPrintf( FG_ASTRO, FG_INFO, "  reading datafile %s\n", path);
 155
 156    /* for all the objects... */
 157    for (i = 0; i < 9; i ++)
 158    {
 159       /* ...read from the data file ... */
 160       fgReadOrbElements(&pltOrbElements[i], data);
 161       /* ...and calculate the actual values */
 162       fgSolarSystemUpdate(&pltOrbElements[i], t);
 163    }
 164
 165 }
 166
 167
 168 void fgSolarSystemUpdate(struct OrbElements *planet, struct fgTIME t)
 169 {
 170    double
 171          actTime;
 172
 173    actTime = fgCalcActTime(t);
 174
 175    /* calculate the actual orbital elements */
 176     planet->M = DEG_TO_RAD * (planet->MFirst + (planet->MSec * actTime));       // angle in radians
 177     planet->w = DEG_TO_RAD * (planet->wFirst + (planet->wSec * actTime));       // angle in radians
 178     planet->N = DEG_TO_RAD * (planet->NFirst + (planet->NSec * actTime));       // angle in radians
 179     planet->i = DEG_TO_RAD * (planet->iFirst + (planet->iSec * actTime));  // angle in radians
 180     planet->e = planet->eFirst + (planet->eSec * actTime);
 181     planet->a = planet->aFirst + (planet->aSec * actTime);
 182 }
 183
 184
 185 /* $Log$
 186 /* Revision 1.6  1998/02/03 23:20:11  curt
 187 /* Lots of little tweaks to fix various consistency problems discovered by
 188 /* Solaris' CC.  Fixed a bug in fg_debug.c with how the fgPrintf() wrapper
 189 /* passed arguments along to the real printf().  Also incorporated HUD changes
 190 /* by Michele America.
 191 /*
 192  * Revision 1.5  1998/02/02 20:53:22  curt
 193  * To version 0.29
 194  *
 195  * Revision 1.4  1998/01/27 00:47:47  curt
 196  * Incorporated Paul Bleisch's <bleisch@chromatic.com> new debug message
 197  * system and commandline/config file processing code.
 198  *
 199  * Revision 1.3  1998/01/22 02:59:27  curt
 200  * Changed #ifdef FILE_H to #ifdef _FILE_H
 201  *
 202  * Revision 1.2  1998/01/19 19:26:58  curt
 203  * Merged in make system changes from Bob Kuehne <rpk@sgi.com>
 204  * This should simplify things tremendously.
 205  *
 206  * Revision 1.1  1998/01/07 03:16:17  curt
 207  * Moved from .../Src/Scenery/ to .../Src/Astro/
 208  *
 209  * Revision 1.6  1997/12/30 20:47:52  curt
 210  * Integrated new event manager with subsystem initializations.
 211  *
 212  * Revision 1.5  1997/12/15 23:55:02  curt
 213  * Add xgl wrappers for debugging.
 214  * Generate terrain normals on the fly.
 215  *
 216  * Revision 1.4  1997/12/10 22:37:51  curt
 217  * Prepended "fg" on the name of all global structures that didn't have it yet.
 218  * i.e. "struct WEATHER {}" became "struct fgWEATHER {}"
 219  *
 220  * Revision 1.3  1997/11/25 23:20:44  curt
 221  * Changed planets.dat Planets.dat
 222  *
 223  * Revision 1.2  1997/11/25 19:25:36  curt
 224  * Changes to integrate Durk's moon/sun code updates + clean up.
 225  *
 226  * Revision 1.1  1997/10/25 03:16:10  curt
 227  * Initial revision of code contributed by Durk Talsma.
 228  *
 229  */