--- /dev/null
+/**************************************************************************
+ * 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 <math.h>
+#include <GL/glut.h>
+
+#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.
+/*
+ */
--- /dev/null
+/**************************************************************************
+ * 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.
+/*
+ */
--- /dev/null
+/**************************************************************************
+ * 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 <string.h>
+
+#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.
+/*
+ */
--- /dev/null
+/**************************************************************************
+ * 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 <stdio.h>
+#include <math.h>
+
+#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.
+/*
+ */
--- /dev/null
+/**************************************************************************
+ * 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.
+/*
+ */
--- /dev/null
+/**************************************************************************
+ * 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.
+/*
+ */
--- /dev/null
+/**************************************************************************
+ * 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.
+/*
+ */
--- /dev/null
+/**************************************************************************
+ * 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.
+/*
+ */
projects / simgear.git / commitdiff