From 0c270d61d85e818fefdb5be155c79e0b743c347c Mon Sep 17 00:00:00 2001 From: curt Date: Sat, 25 Oct 1997 03:16:08 +0000 Subject: [PATCH] Initial revision of code contributed by Durk Talsma. --- Scenery/moon.c | 171 ++++++++++++++++++++++++++++++++++++++++++++++ Scenery/moon.h | 54 +++++++++++++++ Scenery/orbits.c | 170 +++++++++++++++++++++++++++++++++++++++++++++ Scenery/orbits.h | 85 +++++++++++++++++++++++ Scenery/planets.c | 91 ++++++++++++++++++++++++ Scenery/planets.h | 41 +++++++++++ Scenery/sun.c | 61 +++++++++++++++++ Scenery/sun.h | 40 +++++++++++ 8 files changed, 713 insertions(+) create mode 100644 Scenery/moon.c create mode 100644 Scenery/moon.h create mode 100644 Scenery/orbits.c create mode 100644 Scenery/orbits.h create mode 100644 Scenery/planets.c create mode 100644 Scenery/planets.h create mode 100644 Scenery/sun.c create mode 100644 Scenery/sun.h diff --git a/Scenery/moon.c b/Scenery/moon.c new file mode 100644 index 00000000..9b10b3c1 --- /dev/null +++ b/Scenery/moon.c @@ -0,0 +1,171 @@ +/************************************************************************** + * moon.c + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * + * $Id$ + * (Log is kept at end of this file) + **************************************************************************/ + + +#include +#include + +#include "orbits.h" +#include "moon.h" + +#include "../Time/fg_time.h" +#include "../GLUT/views.h" +/* #include "../Aircraft/aircraft.h"*/ +#include "../general.h" + + +static GLint moon; + +struct CelestialCoord fgCalculateMoon(struct OrbElements params, + struct OrbElements sunParams, + struct fgTIME t) +{ + struct CelestialCoord + result; + + double + eccAnom, ecl, lonecl, latecl, actTime, + xv, yv, v, r, xh, yh, zh, xg, yg, zg, xe, ye, ze, + Ls, Lm, D, F; + + /* calculate the angle between ecliptic and equatorial coordinate system */ + actTime = fgCalcActTime(t); + ecl = fgDegToRad(23.4393 - 3.563E-7 * actTime); // in radians of course + + /* calculate the eccentric anomaly */ + eccAnom = fgCalcEccAnom(params.M, params.e); + + /* calculate the moon's distance (d) and true anomaly (v) */ + xv = params.a * ( cos(eccAnom) - params.e); + yv = params.a * ( sqrt(1.0 - params.e*params.e) * sin(eccAnom)); + v =atan2(yv, xv); + r = sqrt(xv*xv + yv*yv); + + /* estimate the geocentric rectangular coordinates here */ + xh = r * (cos(params.N) * cos(v + params.w) - sin(params.N) * sin(v + params.w) * cos(params.i)); + yh = r * (sin(params.N) * cos(v + params.w) + cos(params.N) * sin(v + params.w) * cos(params.i)); + zh = r * (sin(v + params.w) * sin(params.i)); + + /* calculate the ecliptic latitude and longitude here */ + lonecl = atan2( yh, xh); + latecl = atan2( zh, sqrt( xh*xh + yh*yh)); + + /* calculate a number of perturbations */ + Ls = sunParams.M + sunParams.w; + Lm = params.M + params.w + params.N; + D = Lm - Ls; + F = Lm - params.N; + + lonecl += fgDegToRad( + - 1.274 * sin (params.M - 2*D) // the Evection + + 0.658 * sin (2 * D) // the Variation + - 0.186 * sin (sunParams.M) // the yearly variation + - 0.059 * sin (2*params.M - 2*D) + - 0.057 * sin (params.M - 2*D + sunParams.M) + + 0.053 * sin (params.M + 2*D) + + 0.046 * sin (2*D - sunParams.M) + + 0.041 * sin (params.M - sunParams.M) + - 0.035 * sin (D) // the Parallactic Equation + - 0.031 * sin (params.M + sunParams.M) + - 0.015 * sin (2*F - 2*D) + + 0.011 * sin (params.M - 4*D) + ); /* Pheeuuwwww */ + latecl += fgDegToRad( + - 0.173 * sin (F - 2*D) + - 0.055 * sin (params.M - F - 2*D) + - 0.046 * sin (params.M + F - 2*D) + + 0.033 * sin (F + 2*D) + + 0.017 * sin (2 * params.M + F) + ); /* Yep */ + + r += ( + - 0.58 * cos(params.M - 2*D) + - 0.46 * cos(2*D) + ); + xg = r * cos(lonecl) * cos(latecl); + yg = r * sin(lonecl) * cos(latecl); + zg = r * sin(latecl); + + xe = xg; + ye = yg * cos(ecl) - zg * sin(ecl); + ze = yg * sin(ecl) + zg * cos(ecl); + + result.RightAscension = atan2(ye, xe); + result.Declination = atan2(ze, sqrt(xe*xe + ye*ye)); + + return result; +} + + +void fgMoonInit() +{ + struct CelestialCoord + moonPos; + + moon = glGenLists(1); + glNewList(moon, GL_COMPILE ); + glBegin( GL_POINTS ); + moonPos = fgCalculateMoon(pltOrbElements[1], pltOrbElements[0], cur_time_params); + printf("Moon found at %f (ra), %f (dec)\n", moonPos.RightAscension, moonPos.Declination); + /* give the moon a temporary color, for testing purposes */ + glColor3f( 0.0, 1.0, 0.0); + glVertex3f( 190000.0 * cos(moonPos.RightAscension) * cos(moonPos.Declination), + 190000.0 * sin(moonPos.RightAscension) * cos(moonPos.Declination), + 190000.0 * sin(moonPos.Declination) ); + glEnd(); + glEndList(); +} + +void fgMoonRender() +{ + double angle; + static double warp = 0; + struct VIEW *v; + struct fgTIME *t; + + t = &cur_time_params; + v = ¤t_view; + + + glDisable( GL_FOG ); + glDisable( GL_LIGHTING ); + glPushMatrix(); + glTranslatef( v->view_pos.x, v->view_pos.y, v->view_pos.z ); + angle = t->gst * 15.0; /* 15 degrees per hour rotation */ + /* warp += 1.0; */ + /* warp = 15.0; */ + warp = 0.0; + glRotatef( (angle+warp), 0.0, 0.0, -1.0 ); + printf("Rotating moon by %.2f degrees + %.2f degrees\n",angle,warp); + + glCallList(moon); + + glPopMatrix(); + glEnable( GL_LIGHTING ); + glEnable( GL_FOG ); +} + + +/* $Log$ +/* Revision 1.1 1997/10/25 03:16:08 curt +/* Initial revision of code contributed by Durk Talsma. +/* + */ diff --git a/Scenery/moon.h b/Scenery/moon.h new file mode 100644 index 00000000..c43f6dca --- /dev/null +++ b/Scenery/moon.h @@ -0,0 +1,54 @@ +/************************************************************************** + * moon.h + * + * Written 1997 by Durk Talsma, started October, 1997. + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * + * $Id$ + * (Log is kept at end of this file) + **************************************************************************/ + + +#ifndef _MOON_H_ +#define _MOON_H_ + + +#include "orbits.h" + +#include "../Time/fg_time.h" + + + /* Initialize the Moon Display management Subsystem */ +void fgMoonInit(); + +/* Draw the Stars */ + +void fgMoonRender(); +struct CelestialCoord fgCalculateMoon(struct OrbElements Params, + struct OrbElements sunParams, + struct fgTIME t); + +extern struct OrbElements pltOrbElements[9]; + + +#endif /* _MOON_H_ */ + + +/* $Log$ +/* Revision 1.1 1997/10/25 03:16:09 curt +/* Initial revision of code contributed by Durk Talsma. +/* + */ diff --git a/Scenery/orbits.c b/Scenery/orbits.c new file mode 100644 index 00000000..9440b756 --- /dev/null +++ b/Scenery/orbits.c @@ -0,0 +1,170 @@ +/************************************************************************** + * orbits.c - calculates the orbital elements of the sun, moon and planets. + * For inclusion in flight gear + * + * Written 1997 by Durk Talsma, started October 19, 1997. + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * + * $Id$ + * (Log is kept at end of this file) + **************************************************************************/ + + +#include + +#include "orbits.h" + +#include "../general.h" +#include "../Time/fg_time.h" + + +struct OrbElements pltOrbElements[9]; + + +double fgCalcActTime(struct fgTIME t) +{ + double + actTime, UT; + int year; + + /* a hack. This one introduces the 2000 problem into the program */ + year = t.gmt->tm_year + 1900; + + /* calculate the actual time */ + actTime = 367 * year - 7 * + (year + (t.gmt->tm_mon + 9) / 12) / 4 + 275 * + t.gmt->tm_mon / 9 + t.gmt->tm_mday - 730530; + + UT = t.gmt->tm_hour + ((double) t.gmt->tm_min / 60); + /*printf("UT = %f\n", UT); */ + actTime += (UT / 24.0); + printf("current day = %f\t", actTime); + printf("GMT = %d, %d, %d, %d, %d, %d\n", year, t.gmt->tm_mon, t.gmt->tm_mday, + t.gmt->tm_hour, t.gmt->tm_min, t.gmt->tm_sec); + return actTime; +} + +/* convert degrees to radians */ +double fgDegToRad(double angle) +{ + return (angle * PIOVER180); +} + +double fgCalcEccAnom(double M, double e) +{ + double + eccAnom, E0, E1, diff; + + eccAnom = M + e * sin(M) * (1.0 + e * cos(M)); + /* iterate to achieve a greater precision for larger eccentricities */ + if (e > 0.05) + { + E0 = eccAnom; + do + { + E1 = E0 - (E0 - e * sin(E0) - M) / (1 - e * cos(E0)); + diff = abs(E0 - E1); + E0 = E1; + } + while (diff > fgDegToRad(0.001)); + return E0; + } + return eccAnom; +} + + + +void fgReadOrbElements(struct OrbElements *dest, FILE *src) +{ + char line[256]; + int i,j; + j = 0; + do + { + fgets(line, 256,src); + for (i = 0; i < 256; i++) + { + if (line[i] == '#') + line[i] = 0; + } + /*printf("Reading line %d\n", j++); */ + } + while (!(strlen(line))); + sscanf(line, "%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf\n", + &dest->NFirst, &dest->NSec, + &dest->iFirst, &dest->iSec, + &dest->wFirst, &dest->wSec, + &dest->aFirst, &dest->aSec, + &dest->eFirst, &dest->eSec, + &dest->MFirst, &dest->MSec); +} + + +void fgSolarSystemInit(struct fgTIME t) +{ + struct GENERAL *g; + char path[80]; + int i; + FILE *data; + + /* build the full path name to the orbital elements database file */ + g = &general; + path[0] = '\0'; + strcat(path, g->root_dir); + strcat(path, "/Scenery/"); + strcat(path, "planets.dat"); + + if ( (data = fopen(path, "r")) == NULL ) + { + printf("Cannot open data file: '%s'\n", path); + return; + } + printf("reading datafile %s", path); + + /* for all the objects... */ + for (i = 0; i < 9; i ++) + { + /* ...read from the data file ... */ + fgReadOrbElements(&pltOrbElements[i], data); + /* ...and calculate the actual values */ + fgSolarSystemUpdate(&pltOrbElements[i], t); + } + +} + + +void fgSolarSystemUpdate(struct OrbElements *planet, struct fgTIME t) +{ + double + actTime; + + actTime = fgCalcActTime(t); + + /* calculate the actual orbital elements */ + planet->M = fgDegToRad(planet->MFirst + (planet->MSec * actTime)); // angle in radians + planet->w = fgDegToRad(planet->wFirst + (planet->wSec * actTime)); // angle in radians + planet->N = fgDegToRad(planet->NFirst + (planet->NSec * actTime)); // angle in radians + planet->i = fgDegToRad(planet->iFirst + (planet->iSec * actTime)); // angle in radians + planet->e = planet->eFirst + (planet->eSec * actTime); + planet->a = planet->aFirst + (planet->aSec * actTime); +} + + +/* $Log$ +/* Revision 1.1 1997/10/25 03:16:10 curt +/* Initial revision of code contributed by Durk Talsma. +/* + */ diff --git a/Scenery/orbits.h b/Scenery/orbits.h new file mode 100644 index 00000000..2ce03b12 --- /dev/null +++ b/Scenery/orbits.h @@ -0,0 +1,85 @@ +/************************************************************************** + * orbits.h - calculates the orbital elements of the sun, moon and planets. + * For inclusion in flight gear + * + * Written 1997 by Durk Talsma, started October 19, 1997. + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * + * $Id$ + * (Log is kept at end of this file) + **************************************************************************/ + + +#ifndef ORBITS_H +#define ORBITS_H + + +#include +#include + +#include "../Time/fg_time.h" + + + +#define STANDARDEPOCH 2000 +#define PIOVER180 1.74532925199433E-002 + +struct SunPos { + double xs; + double ys; +}; + +struct OrbElements { + double NFirst; /* longitude of the ascending node first part */ + double NSec; /* longitude of the ascending node second part */ + double iFirst; /* inclination to the ecliptic first part */ + double iSec; /* inclination to the ecliptic second part */ + double wFirst; /* first part of argument of perihelion */ + double wSec; /* second part of argument of perihelion */ + double aFirst; /* semimayor axis first part*/ + double aSec; /* semimayor axis second part */ + double eFirst; /* eccentricity first part */ + double eSec; /* eccentricity second part */ + double MFirst; /* Mean anomaly first part */ + double MSec; /* Mean anomaly second part */ + + double N, i, w, a, e, M; /* the resultant orbital elements, obtained from the former */ +}; + +struct CelestialCoord { + double RightAscension; + double Declination; + double distance; +}; + + +double fgDegToRad(double angle); +double fgCalcEccAnom(double M, double e); +double fgCalcActTime(struct fgTIME t); + +void fgReadOrbElements(struct OrbElements *dest, FILE *src); +void fgSolarSystemInit(struct fgTIME t); +void fgSolarSystemUpdate(struct OrbElements *planets, struct fgTIME t); + + +#endif /* ORBITS_H */ + + +/* $Log$ +/* Revision 1.1 1997/10/25 03:16:10 curt +/* Initial revision of code contributed by Durk Talsma. +/* + */ diff --git a/Scenery/planets.c b/Scenery/planets.c new file mode 100644 index 00000000..9f7911a0 --- /dev/null +++ b/Scenery/planets.c @@ -0,0 +1,91 @@ +/************************************************************************** + * planets.c + * + * Written 1997 by Durk Talsma, started October, 1997. For the flight gear + * project. + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * + * $Id$ + * (Log is kept at end of this file) + **************************************************************************/ + + +#include "../Time/fg_time.h" +#include "orbits.h" +#include "planets.h" +#include "sun.h" + + +struct CelestialCoord fgCalculatePlanet(struct OrbElements planet, + struct OrbElements sun, + struct fgTIME t) +{ + struct CelestialCoord + result; + + struct SunPos + SolarPosition; + + double + eccAnom, r, v, ecl, actTime, + xv, yv, xh, yh, zh, xg, yg, zg, xe, ye, ze; + + actTime = fgCalcActTime(t); + /* calculate the angle between ecliptic and equatorial coordinate system */ + ecl = fgDegToRad(23.4393 - 3.563E-7 * actTime); + + + /* calculate the eccentric anomaly */ + eccAnom = fgCalcEccAnom(planet.M, planet.e); + + /* calculate the planets distance (r) and true anomaly (v) */ + xv = planet.a * (cos(eccAnom) - planet.e); + yv = planet.a * (sqrt(1.0 - planet.e*planet.e) * sin(eccAnom)); + v = atan2(yv, xv); + r = sqrt ( xv*xv + yv*yv); + + /* calculate the planets position in 3-dimensional space */ + xh = r * ( cos(planet.N) * cos(v+planet.w) - sin(planet.N) * sin(v+planet.w) * cos(planet.i)); + yh = r * ( sin(planet.N) * cos(v+planet.w) + cos(planet.N) * sin(v+planet.w) * cos(planet.i)); + zh = r * ( sin(v+planet.w) * sin(planet.i)); + + /* calculate the ecleptic longitude and latitude */ + + /* + lonecl = atan2(yh, xh); + latecl = atan2(zh, sqrt ( xh*xh + yh*yh)); + */ + /* calculate the solar position */ + + SolarPosition = fgCalcSunPos(sun); + xg = xh + SolarPosition.xs; + yg = yh + SolarPosition.ys; + zg = zh; + + xe = xg; + ye = yg * cos(ecl) - zg * sin(ecl); + ze = yg * sin(ecl) + zg * cos(ecl); + result.RightAscension = atan2(ye,xe); + result.Declination = atan2(ze, sqrt(xe*xe + ye*ye)); + return result; +} + + +/* $Log$ +/* Revision 1.1 1997/10/25 03:16:10 curt +/* Initial revision of code contributed by Durk Talsma. +/* + */ diff --git a/Scenery/planets.h b/Scenery/planets.h new file mode 100644 index 00000000..e0ce0541 --- /dev/null +++ b/Scenery/planets.h @@ -0,0 +1,41 @@ +/************************************************************************** + * planets.h + * + * Written 1997 by Durk Talsma, started October, 1997. For the flight gear + * project. + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * + * $Id$ + * (Log is kept at end of this file) + **************************************************************************/ + + +#ifndef PLANETS_H +#define PLANETS_H + + +struct CelestialCoord fgCalculatePlanet(struct OrbElements planet, + struct OrbElements sun, + struct fgTIME t); + +#endif /* PLANETS_H */ + + +/* $Log$ +/* Revision 1.1 1997/10/25 03:16:11 curt +/* Initial revision of code contributed by Durk Talsma. +/* + */ diff --git a/Scenery/sun.c b/Scenery/sun.c new file mode 100644 index 00000000..6f1a213f --- /dev/null +++ b/Scenery/sun.c @@ -0,0 +1,61 @@ +/************************************************************************** + * sun.c + * + * Written 1997 by Durk Talsma, started October, 1997. For the flight gear + * project. + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * + * $Id$ + * (Log is kept at end of this file) + **************************************************************************/ + + +#include "../Time/fg_time.h" +#include "orbits.h" + +struct SunPos fgCalcSunPos(struct OrbElements params) +{ + double + EccAnom, lonSun, + xv, yv, v, r; + struct SunPos + solarPosition; + + /* calculate the eccentric anomaly */ + EccAnom = fgCalcEccAnom(params.M, params.e); + + /* calculate the Suns distance (r) and its true anomaly (v) */ + xv = cos(EccAnom) - params.e; + yv = sqrt(1.0 - params.e*params.e) * sin(EccAnom); + v = atan2(yv, xv); + r = sqrt(xv*xv + yv*yv); + + /* calculate the the Suns true longitude (lonsun) */ + lonSun = v + params.w; + + /* convert true longitude and distance to ecliptic rectangular geocentric + coordinates (xs, ys) */ + solarPosition.xs = r * cos(lonSun); + solarPosition.ys = r * sin(lonSun); + return solarPosition; +} + + +/* $Log$ +/* Revision 1.1 1997/10/25 03:16:11 curt +/* Initial revision of code contributed by Durk Talsma. +/* + */ diff --git a/Scenery/sun.h b/Scenery/sun.h new file mode 100644 index 00000000..5082a67c --- /dev/null +++ b/Scenery/sun.h @@ -0,0 +1,40 @@ +/************************************************************************** + * sun.h + * + * Written 1997 by Durk Talsma, started October, 1997. For the flight gear + * project. + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * + * $Id$ + * (Log is kept at end of this file) + **************************************************************************/ + + +#ifndef SUN_H +#define SUN_H + + +struct SunPos fgCalcSunPos(struct OrbElements sunParams); + + +#endif /* SUN_H */ + + +/* $Log$ +/* Revision 1.1 1997/10/25 03:16:12 curt +/* Initial revision of code contributed by Durk Talsma. +/* + */