elliptic_to_rectangular.c | elliptic_to_rectangular.c | |||
---|---|---|---|---|

/************************************************************************ | /************************************************************************ | |||

The code in this file is heavily inspired by the TASS17 and GUST86 theories | The code in this file is heavily inspired by the TASS17 and GUST86 theories | |||

found on | found on | |||

ftp://ftp.imcce.fr/pub/ephem/satel | ftp://ftp.imcce.fr/pub/ephem/satel | |||

I (Johannes Gajdosik) have just taken the Fortran code and data | I (Johannes Gajdosik) have just taken the Fortran code and data | |||

obtained from above and rearranged it into this piece of software. | obtained from above and rearranged it into this piece of software. | |||

I can neigther allow nor forbid the above theories. | I can neither allow nor forbid the above theories. | |||

The copyright notice below covers just my work, | The copyright notice below covers just my work, | |||

that is the compilation of the data obtained from above | that is the compilation of the data obtained from above | |||

into the software supplied in this file. | into the software supplied in this file. | |||

Copyright (c) 2005 Johannes Gajdosik | Copyright (c) 2005 Johannes Gajdosik | |||

Permission is hereby granted, free of charge, to any person obtaining a | Permission is hereby granted, free of charge, to any person obtaining a | |||

copy of this software and associated documentation files (the "Software"), | copy of this software and associated documentation files (the "Software"), | |||

to deal in the Software without restriction, including without limitation | to deal in the Software without restriction, including without limitation | |||

the rights to use, copy, modify, merge, publish, distribute, sublicense, | the rights to use, copy, modify, merge, publish, distribute, sublicense, | |||

skipping to change at line 45 | skipping to change at line 45 | |||

****************************************************************/ | ****************************************************************/ | |||

#include "elliptic_to_rectangular.h" | #include "elliptic_to_rectangular.h" | |||

#include <math.h> | #include <math.h> | |||

#ifndef M_PI | #ifndef M_PI | |||

#define M_PI 3.14159265358979323846 | #define M_PI 3.14159265358979323846 | |||

#endif | #endif | |||

/* | ||||

Given the orbital elements at some time t0 calculate the | ||||

rectangular coordinates at time (t0+dt). | ||||

mu = G*(m1+m2) .. gravitational constant of the two body problem | ||||

a .. semi major axis | ||||

n = mean motion = 2*M_PI/(orbit period) | ||||

elem[0] .. irrelevant (either a (called by EllipticToRectangularA()) or | ||||

n (called by EllipticToRectangularN()) | ||||

elem[1] .. L | ||||

elem[2] .. K=e*cos(Omega+omega) | ||||

elem[3] .. H=e*sin(Omega+omega) | ||||

elem[4] .. Q=sin(i/2)*cos(Omega) | ||||

elem[5] .. P=sin(i/2)*sin(Omega) | ||||

Omega = longitude of ascending node | ||||

omega = argument of pericenter | ||||

L = mean longitude = Omega + omega + M | ||||

M = mean anomaly | ||||

i = inclination | ||||

e = excentricity | ||||

Units are suspected to be: Julian days, AU, rad | ||||

*/ | ||||

static void | static void | |||

EllipticToRectangular(const double a,const double n, | EllipticToRectangular(const double a,const double n, | |||

const double elem[6],const double dt,double xyz[]) { | const double elem[6],const double dt,double xyz[]) { | |||

const double L = fmod(elem[1]+n*dt,2.0*M_PI); | const double L = fmod(elem[1]+n*dt,2.0*M_PI); | |||

/* solve Keplers equation | /* solve Keplers equation | |||

x = L - elem[2]*sin(x) + elem[3]*cos(x) | x = L - elem[2]*sin(x) + elem[3]*cos(x) | |||

not by trivially iterating | not by trivially iterating | |||

x_0 = L | x_0 = L | |||

x_{j+1} = L - elem[2]*sin(x_j) + elem[3]*cos(x_j) | x_{j+1} = L - elem[2]*sin(x_j) + elem[3]*cos(x_j) | |||

but instead by Newton approximation: | but instead by Newton approximation: | |||

0 = f(x) = x - L - elem[2]*sin(x) + elem[3]*cos(x) | 0 = f(x) = x - L - elem[2]*sin(x) + elem[3]*cos(x) | |||

f'(x) = 1 - elem[2]*cos(x) - elem[3]*sin(x) | f'(x) = 1 - elem[2]*cos(x) - elem[3]*sin(x) | |||

x_0 = L or whatever, perhaps first step of trivial iteration | x_0 = L or whatever, perhaps first step of trivial iteration | |||

x_{j+1} = x_j - f(x_j)/f'(x_j) | x_{j+1} = x_j - f(x_j)/f'(x_j) | |||

*/ | */ | |||

double Le = L - elem[2]*sin(L) + elem[3]*cos(L); | double Le = L - elem[2]*sin(L) + elem[3]*cos(L); | |||

for (;;) { | for (;;) { | |||

const double cLe = cos(Le); | const double cLe = cos(Le); | |||

const double sLe = sin(Le); | const double sLe = sin(Le); | |||

/* for excenticity < 1 we have denominator > 0 */ | /* for excentricity < 1 we have denominator > 0 */ | |||

const double dLe = (L - Le + elem[2]*sLe - elem[3]*cLe) | const double dLe = (L - Le + elem[2]*sLe - elem[3]*cLe) | |||

/ (1.0 - elem[2]*cLe - elem[3]*sLe); | / (1.0 - elem[2]*cLe - elem[3]*sLe); | |||

Le += dLe; | Le += dLe; | |||

if (fabs(dLe) <= 1e-14) break; /* L1: <1e-12 */ | if (fabs(dLe) <= 1e-14) break; /* L1: <1e-12 */ | |||

} | } | |||

{ | { | |||

const double cLe = cos(Le); | const double cLe = cos(Le); | |||

const double sLe = sin(Le); | const double sLe = sin(Le); | |||

const double dlf = -elem[2]*sLe + elem[3]*cLe; | const double dlf = -elem[2]*sLe + elem[3]*cLe; | |||

const double phi = sqrt(1.0 - elem[2]*elem[2] - elem[3]*elem[3]); | const double phi = sqrt(1.0 - elem[2]*elem[2] - elem[3]*elem[3]); | |||

const double psi = 1.0 / (1.0 + phi); | const double psi = 1.0 / (1.0 + phi); | |||

const double x1 = a * (cLe - elem[2] - psi*dlf*elem[3]); | const double x1 = a * (cLe - elem[2] - psi*dlf*elem[3]); | |||

const double y1 = a * (sLe - elem[3] + psi*dlf*elem[2]); | const double y1 = a * (sLe - elem[3] + psi*dlf*elem[2]); | |||

const double elem_4q = elem[4] * elem[4]; | const double elem_4q = elem[4] * elem[4]; // Q² | |||

const double elem_5q = elem[5] * elem[5]; | const double elem_5q = elem[5] * elem[5]; // P² | |||

const double dwho = 2.0 * sqrt(1.0 - elem_4q - elem_5q); | const double dwho = 2.0 * sqrt(1.0 - elem_4q - elem_5q); | |||

const double rtp = 1.0 - elem_5q - elem_5q; | const double rtp = 1.0 - elem_5q - elem_5q; | |||

const double rtq = 1.0 - elem_4q - elem_4q; | const double rtq = 1.0 - elem_4q - elem_4q; | |||

const double rdg = 2.0 * elem[5] * elem[4]; | const double rdg = 2.0 * elem[5] * elem[4]; | |||

xyz[0] = x1 * rtp + y1 * rdg; | xyz[0] = x1 * rtp + y1 * rdg; | |||

xyz[1] = x1 * rdg + y1 * rtq; | xyz[1] = x1 * rdg + y1 * rtq; | |||

xyz[2] = (-x1 * elem[5] + y1 * elem[4]) * dwho; | xyz[2] = (-x1 * elem[5] + y1 * elem[4]) * dwho; | |||

/* | /* | |||

skipping to change at line 121 | skipping to change at line 145 | |||

const double a = cbrt(mu/(n*n)); | const double a = cbrt(mu/(n*n)); | |||

#else | #else | |||

const double a = exp(log(mu/(n*n))/3.0); | const double a = exp(log(mu/(n*n))/3.0); | |||

#endif | #endif | |||

EllipticToRectangular(a,n,elem,dt,xyz); | EllipticToRectangular(a,n,elem,dt,xyz); | |||

} | } | |||

void EllipticToRectangularA(double mu,const double elem[6],double dt, | void EllipticToRectangularA(double mu,const double elem[6],double dt, | |||

double xyz[]) { | double xyz[]) { | |||

const double a = elem[0]; | const double a = elem[0]; | |||

const double n = sqrt(mu/(a*a*a)); | const double n = sqrt(mu/(a*a*a)); // mean motion | |||

EllipticToRectangular(a,n,elem,dt,xyz); | EllipticToRectangular(a,n,elem,dt,xyz); | |||

} | } | |||

End of changes. 5 change blocks. | ||||

5 lines changed or deleted | | 30 lines changed or added | ||

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/ |