1 /**************************************************************************
4 * Written 1997 by Durk Talsma, started October, 1997. For the flight gear
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.
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.
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.
22 * (Log is kept at end of this file)
23 **************************************************************************/
26 #include "../Time/fg_time.h"
32 struct CelestialCoord fgCalculatePlanet(struct OrbElements planet,
33 struct OrbElements sun,
43 eccAnom, r, v, ecl, actTime,
44 xv, yv, xh, yh, zh, xg, yg, zg, xe, ye, ze;
46 actTime = fgCalcActTime(t);
47 /* calculate the angle between ecliptic and equatorial coordinate system */
48 ecl = fgDegToRad(23.4393 - 3.563E-7 * actTime);
51 /* calculate the eccentric anomaly */
52 eccAnom = fgCalcEccAnom(planet.M, planet.e);
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));
58 r = sqrt ( xv*xv + yv*yv);
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));
65 /* calculate the ecleptic longitude and latitude */
68 lonecl = atan2(yh, xh);
69 latecl = atan2(zh, sqrt ( xh*xh + yh*yh));
71 /* calculate the solar position */
73 SolarPosition = fgCalcSunPos(sun);
74 xg = xh + SolarPosition.xs;
75 yg = yh + SolarPosition.ys;
79 ye = yg * cos(ecl) - zg * sin(ecl);
80 ze = yg * sin(ecl) + zg * cos(ecl);
81 result.RightAscension = atan2(ye,xe);
82 result.Declination = atan2(ze, sqrt(xe*xe + ye*ye));
88 /* Revision 1.1 1997/10/25 03:16:10 curt
89 /* Initial revision of code contributed by Durk Talsma.