1
/************************************************************************
3
The code in this file is heavily inspired by the TASS17 and GUST86 theories
5
ftp://ftp.imcce.fr/pub/ephem/satel
7
I (Johannes Gajdosik) have just taken the Fortran code and data
8
obtained from above and rearranged it into this piece of software.
10
I can neigther allow nor forbid the above theories.
11
The copyright notice below covers just my work,
12
that is the compilation of the data obtained from above
13
into the software supplied in this file.
16
Copyright (c) 2005 Johannes Gajdosik
18
Permission is hereby granted, free of charge, to any person obtaining a
19
copy of this software and associated documentation files (the "Software"),
20
to deal in the Software without restriction, including without limitation
21
the rights to use, copy, modify, merge, publish, distribute, sublicense,
22
and/or sell copies of the Software, and to permit persons to whom the
23
Software is furnished to do so, subject to the following conditions:
25
The above copyright notice and this permission notice shall be included
26
in all copies or substantial portions of the Software.
28
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
29
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
30
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
31
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
32
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
33
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
36
****************************************************************/
38
#include "elliptic_to_rectangular.h"
43
EllipticToRectangular(const double mu,const double a,const double n,
44
const double elem[6],const double dt,double xyz[]) {
45
const double L = fmod(elem[1]+n*dt,2.0*M_PI);
46
/* solve Keplers equation
47
x = L - elem[2]*sin(x) + elem[3]*cos(x)
48
not by trivially iterating
50
x_{j+1} = L - elem[2]*sin(x_j) + elem[3]*cos(x_j)
51
but instead by Newton approximation:
52
0 = f(x) = x - L - elem[2]*sin(x) + elem[3]*cos(x)
53
f'(x) = 1 - elem[2]*cos(x) - elem[3]*sin(x)
54
x_0 = L or whatever, perhaps first step of trivial iteration
55
x_{j+1} = x_j - f(x_j)/f'(x_j)
57
double Le = L - elem[2]*sin(L) + elem[3]*cos(L);
59
const double cLe = cos(Le);
60
const double sLe = sin(Le);
61
/* for excenticity < 1 we have denominator > 0 */
62
const double dLe = (L - Le + elem[2]*sLe - elem[3]*cLe)
63
/ (1.0 - elem[2]*cLe - elem[3]*sLe);
65
if (fabs(dLe) <= 1e-14) break; /* L1: <1e-12 */
68
const double cLe = cos(Le);
69
const double sLe = sin(Le);
71
const double dlf = -elem[2]*sLe + elem[3]*cLe;
72
const double phi = sqrt(1.0 - elem[2]*elem[2] - elem[3]*elem[3]);
73
const double psi = 1.0 / (1.0 + phi);
75
const double x1 = a * (cLe - elem[2] - psi*dlf*elem[3]);
76
const double y1 = a * (sLe - elem[3] + psi*dlf*elem[2]);
78
const double elem_4q = elem[4] * elem[4];
79
const double elem_5q = elem[5] * elem[5];
80
const double dwho = 2.0 * sqrt(1.0 - elem_4q - elem_5q);
81
const double rtp = 1.0 - elem_5q - elem_5q;
82
const double rtq = 1.0 - elem_4q - elem_4q;
83
const double rdg = 2.0 * elem[5] * elem[4];
85
xyz[0] = x1 * rtp + y1 * rdg;
86
xyz[1] = x1 * rdg + y1 * rtq;
87
xyz[2] = (-x1 * elem[5] + y1 * elem[4]) * dwho;
90
const double rsam1 = -elem[2]*cLe - elem[3]*sLe;
91
const double h = a*n / (1.0 + rsam1);
92
const double vx1 = h * (-sLe - psi*rsam1*elem[3]);
93
const double vy1 = h * ( cLe + psi*rsam1*elem[2]);
95
xyz[3] = vx1 * rtp + vy1 * rdg;
96
xyz[4] = vx1 * rdg + vy1 * rtq;
97
xyz[5] = (-vx1 * elem[5] + vy1 * elem[4]) * dwho;
101
void EllipticToRectangularN(double mu,const double elem[6],double dt,
103
const double n = elem[0];
104
const double a = cbrt(mu/(n*n));
105
EllipticToRectangular(mu,a,n,elem,dt,xyz);
108
void EllipticToRectangularA(double mu,const double elem[6],double dt,
110
const double a = elem[0];
111
const double n = sqrt(mu/(a*a*a));
112
EllipticToRectangular(mu,a,n,elem,dt,xyz);
added
removed
src/stellplanet/elliptic_to_rectangular.c
