home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Gambler 19
/
GAMBLERCD19.BIN
/
UTILS
/
RIVA128
/
Uzytki
/
ssystem-1.2
/
SOURCE
/
positions.c
< prev
next >
Wrap
C/C++ Source or Header
|
1998-05-02
|
8KB
|
252 lines
/*
Calculate planets positions in a given date.
The main algorithm uses elements dated Julian 2450680.5 (8/20/1997) so
it's quite accurate for dates near this.
To be improved in the future ...
Galilean and Saturn satellites computed using code extracted with
permission from XEphem, (c) 1997 Elwood Charles Downey.
Keith Burnett's planet positioning algorithm used with persmission.
*/
#include "planets.h"
static void GalileanSat(void),SaturnSat(void),NeptuneSat(void),PlutoSat(void);
static double SolveKepler(double, double);
/* Iterative solution of Kepler equation */
static double SolveKepler(double M, double ec)
{
static double E,d,delta;
E=M;
d=E-ec*sin(E)-M;
while (fabs(d)>1.0E-10) {
delta=d/(1.0-ec*cos(E));
E=E-delta;
d=E-ec*sin(E)-M;
}
return E;
}
/* This function is nearly a clone of the jupiter_data function in xephem 3.0
uses elements dated 12/31/1899
*/
static void GalileanSat()
{
#define dsin(x) sin(DEG2RAD(x))
#define dcos(x) cos(DEG2RAD(x))
static double A, B, Del, J, K, M, N, R, V;
static double cor_u1, cor_u2, cor_u3, cor_u4;
static double tmp, GG, H, psi, r, r1, r2, r3, r4;
static double u1, u2, u3, u4;
static double lam, Ds;
static double z1, z2, z3, z4;
static double De, dsinDe,d;
d=days+35295.0;
V = 134.63 + 0.00111587 * d;
M = (358.47583 + 0.98560003*d);
N = (225.32833 + 0.0830853*d) + 0.33 * dsin (V);
J = 221.647 + 0.9025179*d - 0.33 * dsin(V);;
A = 1.916*dsin(M) + 0.02*dsin(2*M);
B = 5.552*dsin(N) + 0.167*dsin(2*N);
K = (J+A-B);
R = 1.00014 - 0.01672 * dcos(M) - 0.00014 * dcos(2*M);
r = 5.20867 - 0.25192 * dcos(N) - 0.00610 * dcos(2*N);
Del = sqrt (R*R + r*r - 2*R*r*dcos(K));
psi = RAD2DEG (asin(R/Del*dsin(K)));
tmp = psi - B;
u1 = 84.5506 + 203.4058630 * d + tmp;
u2 = 41.5015 + 101.2916323 * d + tmp;
u3 = 109.9770 + 50.2345169 * d + tmp;
u4 = 176.3586 + 21.4879802 * d + tmp;
GG = 187.3 + 50.310674 * d;
H = 311.1 + 21.569229 * d;
cor_u1 = 0.472 * dsin (2*(u1-u2));
cor_u2 = 1.073 * dsin (2*(u2-u3));
cor_u3 = 0.174 * dsin (GG);
cor_u4 = 0.845 * dsin (H);
r1 = RADIUSSCALE(planets[JUPITER].Radius)*5.9061 - 0.0244 * dcos (2*(u1-u2));
r2 = RADIUSSCALE(planets[JUPITER].Radius)*9.3972 - 0.0889 * dcos (2*(u2-u3));
r3 = RADIUSSCALE(planets[JUPITER].Radius)*14.9894 - 0.0227 * dcos (GG);
r4 = RADIUSSCALE(planets[JUPITER].Radius)*26.3649 - 0.1944 * dcos (H);
planets[IO].posx = r1 * dsin (u1+cor_u1);
planets[EUROPA].posx = r2 * dsin (u2+cor_u2);
planets[GANYMEDE].posx = r3 * dsin (u3+cor_u3);
planets[CALLISTO].posx = r4 * dsin (u4+cor_u4);
lam = 238.05 + 0.083091*d + 0.33*dsin(V) + B;
Ds = 3.07*dsin(lam + 44.5);
De = Ds - 2.15*dsin(psi)*dcos(lam+24.)
- 1.31*(r-Del)/Del*dsin(lam-99.4);
dsinDe = dsin(De);
z1 = r1 * dcos(u1+cor_u1);
z2 = r2 * dcos(u2+cor_u2);
z3 = r3 * dcos(u3+cor_u3);
z4 = r4 * dcos(u4+cor_u4);
planets[IO].posy = z1*dsinDe;
planets[EUROPA].posy = z2*dsinDe;
planets[GANYMEDE].posy = z3*dsinDe;
planets[CALLISTO].posy = z4*dsinDe;
planets[IO].posz = z1;
planets[EUROPA].posz = z2;
planets[GANYMEDE].posz = z3;
planets[CALLISTO].posz = z4;
#undef dsin
#undef dcos
}
/* Another xephem based algorithm
Elements dated Jan 0.0 1980 = December 31,1979 0:0:0 UT
*/
static void SaturnSat()
{
static double SMA[11], U[11], U0[11], PD[11];
static double NN;
static int i;
/* Semimajor Axis of the Satellites' Orbit in planet radius */
SMA[1] = 3.08; SMA[2] = 3.95; SMA[3] = 4.89; SMA[4] = 6.26;
SMA[5] = 8.74; SMA[6] = 20.27; SMA[7] = 24.58; SMA[8] = 59.09;
/* Eccentricity of Satellites' Orbit [Program uses 0] */
/* Synodic Orbital Period of Moons in Days */
PD[1] = .9425049; /* Mimas */
PD[2] = 1.3703731; /* Enceladus */
PD[3] = 1.8880926; /* Tethys */
PD[4] = 2.7375218; /* Dione */
PD[5] = 4.5191631; /* Rhea */
PD[6] = 15.9669028; /* Titan */
PD[7] = 21.3174647; /* Hyperion */
PD[8] = 79.9190206; /* Iapetus */
U0[1] = DEG2RAD(18.2919);
U0[2] = DEG2RAD(174.2135);
U0[3] = DEG2RAD(172.8546);
U0[4] = DEG2RAD(76.8438);
U0[5] = DEG2RAD(37.2555);
U0[6] = DEG2RAD(57.7005);
U0[7] = DEG2RAD(266.6977);
U0[8] = DEG2RAD(195.3513);
NN=days+7305.0;
for (i=TETHYS;i<=TITAN;i++) {
U[i-TETHYS+1] = U0[i-TETHYS+1] + (NN * 2 * PI / PD[i-TETHYS+1]);
U[i-TETHYS+1] = ((U[i-TETHYS+1] / (2 * PI)) - (int)(U[i-TETHYS+1] / (2 * PI))) * 2 * PI;
planets[i].posx = SMA[i-TETHYS+1] * sin(U[i-TETHYS+1]) * RADIUSSCALE(planets[SATURN].Radius);
planets[i].posz = SMA[i-TETHYS+1] * cos(U[i-TETHYS+1]) * RADIUSSCALE(planets[SATURN].Radius);
planets[i].posy = -SMA[i-TETHYS+1] * cos(U[i-TETHYS+1]) * sin(DEG2RAD(1.0)) * RADIUSSCALE(planets[SATURN].Radius);
}
}
static void NeptuneSat()
{
static double SMA[11], U[11], U0[11], PD[11];
static double NN;
/* Semimajor Axis of the Satellites' Orbit in planet radius */
SMA[1] = 14.33;
/* Eccentricity of Satellites' Orbit [Program uses 0] */
/* Synodic Orbital Period of Moons in Days */
PD[1] = 5.876854; /* Triton */
U0[1] = 0.0; /* ? */
NN=days;
U[1] = U0[1] + (NN * 2 * PI / PD[1]);
U[1] = ((U[1] / (2 * PI)) - (int)(U[1] / (2 * PI))) * 2 * PI;
planets[TRITON].posx = SMA[1] * sin(U[1]) * RADIUSSCALE(planets[NEPTUNE].Radius);
planets[TRITON].posz = SMA[1] * cos(U[1]) * RADIUSSCALE(planets[NEPTUNE].Radius);
planets[TRITON].posy = -SMA[1] * cos(U[1]) * sin(DEG2RAD(1.0)) * RADIUSSCALE(planets[NEPTUNE].Radius);
}
static void PlutoSat()
{
static double SMA[11], U[11], U0[11], PD[11];
static double NN;
/* Semimajor Axis of the Satellites' Orbit in planet radius */
SMA[1] = 17.0;
/* Eccentricity of Satellites' Orbit [Program uses 0] */
/* Synodic Orbital Period of Moons in Days */
PD[1] = 6.38725; /* Charon */
U0[1] = 0.0; /* ? */
NN=days;
U[1] = U0[1] + (NN * 2 * PI / PD[1]);
U[1] = ((U[1] / (2 * PI)) - (int)(U[1] / (2 * PI))) * 2 * PI;
planets[CHARON].posx = SMA[1] * RADIUSSCALE(planets[PLUTO].Radius)
* sin(U[1]);
planets[CHARON].posz = SMA[1] * RADIUSSCALE(planets[PLUTO].Radius)
* cos(U[1]);
planets[CHARON].posy = -SMA[1] * RADIUSSCALE(planets[PLUTO].Radius)
* cos(U[1]) * sin(DEG2RAD(1.0));
}
/* Based on Keith Burnett's QBASIC code found here:
http://www.xylem.demon.co.uk/kepler/
*/
void UpdatePositions(void)
{
static int j;
static double e,M,E,r,v,o,p,i,x,y,z;
for (j=1;j<=MOON;j++) {
planets[j].DeltaRotation+=planets[j].Rotation*timefactor;
e=planets[j].Eccentricity;
M=planets[j].DailyMotion*days+planets[j].MeanLongitude-planets[j].Perihelion;
E=SolveKepler(M,e);
r=planets[j].MeanDistance*(1.0-e*cos(E));
v=2.0*atan(sqrt((1+e)/(1-e))*tan(E/2.0));
o=planets[j].AscendingNode;
p=planets[j].Perihelion;
i=planets[j].Inclination;
planets[j].posx=r*(cos(o)*cos(v+p-o)-sin(o)*sin(v+p-o)*cos(i));
planets[j].posz=-r*(sin(o)*cos(v+p-o)+cos(o)*sin(v+p-o)*cos(i));
planets[j].posy=r*(sin(v+p-o)*sin(i));
}
GalileanSat();
SaturnSat();
NeptuneSat();
PlutoSat();
for (j=MOON;j<=CHARON;j++) {
x=planets[j].posx;
y=planets[j].posy;
z=planets[j].posz;
o=cos(DEG2RAD(planets[planets[j].Sat].Degrees)+planets[j].Inclination);
p=sin(DEG2RAD(planets[planets[j].Sat].Degrees)+planets[j].Inclination);
planets[j].posx = planets[planets[j].Sat].posx+x;
planets[j].posy = planets[planets[j].Sat].posy+y*o+z*p;
planets[j].posz = planets[planets[j].Sat].posz-y*p+z*o;
}
i=atan2(planets[EARTH].posz,planets[EARTH].posx)*180.0/PI;
planets[EARTH].DeltaRotation=(days-floor(days))*planets[EARTH].Rotation-i+90.0;
}