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>
27 #include <Astro/orbits.h>
28 #include <Astro/planets.h>
29 #include <Astro/sun.h>
32 struct CelestialCoord fgCalculatePlanet(struct OrbElements planet,
33 struct OrbElements theSun,
34 struct fgTIME t, int idx)
43 eccAnom, r, v, ecl, actTime, R, s, ir, Nr, B, FV, ring_magn,
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(theSun);
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);
82 result.RightAscension = atan2(ye,xe);
83 result.Declination = atan2(ze, sqrt(xe*xe + ye*ye));
86 /* Let's calculate the brightness of the planet */
87 R = sqrt ( xg*xg + yg*yg + zg*zg);
88 s = SolarPosition.dist;
89 FV = acos( (r*r + R*R - s*s) / (2*r*R));
90 FV *= 57.29578; /* convert radians to degrees */
94 result.magnitude = -0.36 + 5*log10( r*R ) + 0.027 * FV + 2.2E-13 * pow(FV, 6);
97 result.magnitude = -4.34 + 5*log10( r*R ) + 0.013 * FV + 4.2E-07 * pow(FV,3);
100 result.magnitude = -1.51 + 5*log10( r*R ) + 0.016 * FV;
102 case 5: /* Jupiter */
103 result.magnitude = -9.25 + 5*log10( r*R ) + 0.014 * FV;
108 Nr = 2.9585076 + 6.6672E-7*actTime;
110 B = asin ( sin (result.Declination) * cos(ir) - cos(result.Declination) * sin (ir) * sin (result.RightAscension - Nr));
111 ring_magn = -2.6 * sin (abs(B)) + 1.2 * pow(sin(B),2);
112 result.magnitude = -9.0 + 5*log10( r*R ) + 0.044 * FV + ring_magn;
115 result.magnitude = -7.15 + 5*log10( r*R) + 0.001 * FV;
117 case 8: /* Neptune */
118 result.magnitude = -6.90 + 5*log10 (r*R) + 0.001 *FV;
121 printf("index %d out of range !!!!\n", idx);
123 printf(" Planet found at %f (ra), %f (dec)\n",
124 result.RightAscension, result.Declination);
125 printf(" Geocentric dist %f\n"
126 " Heliocentric dist %f\n"
127 " Distance to the sun %f\n"
129 " Brightness %f\n", R, r, s, FV, result.magnitude);
136 /* Revision 1.2 1998/01/19 19:26:59 curt
137 /* Merged in make system changes from Bob Kuehne <rpk@sgi.com>
138 /* This should simplify things tremendously.
140 * Revision 1.1 1998/01/07 03:16:18 curt
141 * Moved from .../Src/Scenery/ to .../Src/Astro/
143 * Revision 1.4 1997/12/30 20:47:52 curt
144 * Integrated new event manager with subsystem initializations.
146 * Revision 1.3 1997/12/30 16:36:52 curt
147 * Merged in Durk's changes ...
149 * Revision 1.2 1997/12/12 21:41:29 curt
150 * More light/material property tweaking ... still a ways off.
152 * Revision 1.1 1997/10/25 03:16:10 curt
153 * Initial revision of code contributed by Durk Talsma.