Scenery/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 <string.h>
  27
  28 #include "orbits.h"
  29
  30 #include "../general.h"
  31 #include "../Time/fg_time.h"
  32
  33
  34 struct OrbElements pltOrbElements[9];
  35
  36
  37 double fgCalcActTime(struct fgTIME t)
  38 {
  39    double
  40          actTime, UT;
  41    int year;
  42
  43    /* a hack. This one introduces the 2000 problem into the program */
  44    year = t.gmt->tm_year + 1900;
  45
  46    /* calculate the actual time */
  47    actTime = 367 * year - 7 *
  48                   (year + (t.gmt->tm_mon + 9) / 12) / 4 + 275 *
  49                    t.gmt->tm_mon / 9 + t.gmt->tm_mday - 730530;
  50
  51     UT = t.gmt->tm_hour + ((double) t.gmt->tm_min / 60);
  52     /*printf("UT = %f\n", UT); */
  53     actTime += (UT / 24.0);
  54     printf("current day = %f\t", actTime);
  55     printf("GMT = %d, %d, %d, %d, %d, %d\n",           year, t.gmt->tm_mon, t.gmt->tm_mday,
  56                                              t.gmt->tm_hour, t.gmt->tm_min, t.gmt->tm_sec);
  57     return actTime;
  58 }
  59
  60 /* convert degrees to radians */
  61 double fgDegToRad(double angle)
  62 {
  63         return (angle * PIOVER180);
  64 }
  65
  66 double fgCalcEccAnom(double M, double e)
  67 {
  68     double
  69         eccAnom, E0, E1, diff;
  70
  71     eccAnom = M + e * sin(M) * (1.0 + e * cos(M));
  72     /* iterate to achieve a greater precision for larger eccentricities */
  73     if (e > 0.05)
  74     {
  75         E0 = eccAnom;
  76         do
  77         {
  78                 E1 = E0 - (E0 - e * sin(E0) - M) / (1 - e * cos(E0));
  79             diff = abs(E0 - E1);
  80             E0 = E1;
  81                 }
  82         while (diff > fgDegToRad(0.001));
  83         return E0;
  84         }
  85     return eccAnom;
  86 }
  87
  88
  89
  90 void fgReadOrbElements(struct OrbElements *dest, FILE *src)
  91 {
  92         char line[256];
  93     int i,j;
  94     j = 0;
  95     do
  96     {
  97         fgets(line, 256,src);
  98         for (i = 0; i < 256; i++)
  99         {
 100                 if (line[i] == '#')
 101                 line[i] = 0;
 102         }
 103         /*printf("Reading line %d\n", j++); */
 104     }
 105     while (!(strlen(line)));
 106     sscanf(line, "%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf\n",
 107                 &dest->NFirst, &dest->NSec,
 108             &dest->iFirst, &dest->iSec,
 109             &dest->wFirst, &dest->wSec,
 110             &dest->aFirst, &dest->aSec,
 111             &dest->eFirst, &dest->eSec,
 112             &dest->MFirst, &dest->MSec);
 113 }
 114
 115
 116 void fgSolarSystemInit(struct fgTIME t)
 117 {
 118    struct GENERAL *g;
 119    char path[80];
 120    int i;
 121    FILE *data;
 122
 123    /* build the full path name to the orbital elements database file */
 124    g = &general;
 125    path[0] = '\0';
 126    strcat(path, g->root_dir);
 127    strcat(path, "/Scenery/");
 128    strcat(path, "planets.dat");
 129
 130    if ( (data = fopen(path, "r")) == NULL )
 131    {
 132             printf("Cannot open data file: '%s'\n", path);
 133             return;
 134    }
 135    printf("reading datafile %s", path);
 136
 137    /* for all the objects... */
 138    for (i = 0; i < 9; i ++)
 139    {
 140       /* ...read from the data file ... */
 141       fgReadOrbElements(&pltOrbElements[i], data);
 142       /* ...and calculate the actual values */
 143       fgSolarSystemUpdate(&pltOrbElements[i], t);
 144    }
 145
 146 }
 147
 148
 149 void fgSolarSystemUpdate(struct OrbElements *planet, struct fgTIME t)
 150 {
 151    double
 152          actTime;
 153
 154    actTime = fgCalcActTime(t);
 155
 156    /* calculate the actual orbital elements */
 157     planet->M = fgDegToRad(planet->MFirst + (planet->MSec * actTime));  // angle in radians
 158     planet->w = fgDegToRad(planet->wFirst + (planet->wSec * actTime));  // angle in radians
 159     planet->N = fgDegToRad(planet->NFirst + (planet->NSec * actTime));  // angle in radians
 160     planet->i = fgDegToRad(planet->iFirst + (planet->iSec * actTime));  // angle in radians
 161     planet->e = planet->eFirst + (planet->eSec * actTime);
 162     planet->a = planet->aFirst + (planet->aSec * actTime);
 163 }
 164
 165
 166 /* $Log$
 167 /* Revision 1.1  1997/10/25 03:16:10  curt
 168 /* Initial revision of code contributed by Durk Talsma.
 169 /*
 170  */