?

Log in

No account? Create an account
Previous Entry Share Flag Next Entry
HP Prime program to find orbital elements from four observations by the method of Gauss
Jenab6
jenab6
ORBIT4 is a program for the HP Prime. It solves, using the method of Gauss, for the Keplerian orbital elements and for the (initially unknown) geocentric and heliocentric distances of an asteroid by using observational data for four known times. The input for each time consists of the time of observation in Julian Date, the position vector of Earth in rectangular heliocentric ecliptic coordinates, and the right ascension and declination of the asteroid in geocentric coordinates.

It works like a charm.

You can get ORBIT4 as a free download from CloudMe, the superior alternative file host to the unreliable, censorship-minded Google.

If you download and install the following two apps, you can run ORBIT4 on your PC. You won't need to own a physical HP Prime calculator.

HP Connectivity Kit (for Windows PC)
HP Prime Virtual Calculator Emulator (for Windows PC)

Test data for ORBIT4 can be obtained from JPL Horizons Web Interface.

I've placed ORBIT4 in the public domain. It is available to anyone, worldwide, for free at no charge. All I want is credit for it and an acknowledgement in derivative products made by others. The SOURCE CODE is in the following section.

Acknowledgement to the Russian astronomer Alexander Dmitriyevich Dubyago (Александр Дмитриевич Дубяго), whose book The Determination of Orbits taught me the method used in the program.


Source Code

EXPORT ORBIT4()
BEGIN

PRINT("ORBIT4 by David Sims");
PRINT("Method of Gauss with four observed positions");
PRINT("to find the Keplerian orbital elements.");
PRINT("User provides input by adjusting inline data.");
PRINT(" ");

LOCAL α,Δ,δ,ξ,η,κ,σ,β;
LOCAL a,b,c,d,e,f,g,h;
LOCAL i,j,k,l,m,p,q,t;
LOCAL Ĉ,Ŝ,Ω,ω,x,y,z,r;
LOCAL s,ε;

1_rad;

// _____ Data format _____
// {Time in Julian date,
// HEC Earth position X in AU,
// HEC Earth position Y in AU,
// HEC Earth position Z in AU,
// Object GCC rt ascension,
// Object GCC declination}

// Data for time 1
L1:=
{2457204.625,
0.155228396,
−1.004732775,
0.00003295786,
HMS→(20°46′57.02″),
HMS→(−27°41′33.9″)};

// Data for time 2
L2:=
{2457214.625,
0.319493277,
−0.965116604,
0.0000311269,
HMS→(20°39′57.10″),
HMS→(−28°47′21.5″)};

// Data for time 3
L3:=
{2457224.625,
0.4747795623,
−0.8983801739,
0.00002841127,
HMS→(20°31′22.81″),
HMS→(−29°49′22.7″)};

// Data for time 4
L4:=
{2457234.625,
0.616702829,
−0.8063620175,
0.00002486325,
HMS→(20°22′06.57″),
HMS→(−30°41′57.3″)};

// Time intervals should be about
// 0.5% to 1% of the period, should
// be near opposition with the sun,
// but should not span an apside.
// The precision in RA: 0.01 seconds.
// The precision in DEC: 0.1 arcsec.

// κ : Mean motion of Earth, rad/day
0.01720209895▶κ;

// Convert α⒤,δ⒤ to radians
(π/12)*L1[5]▶L1[5];
(π/12)*L2[5]▶L2[5];
(π/12)*L3[5]▶L3[5];
(π/12)*L4[5]▶L4[5];
(π/180)*L1[6]▶L1[6];
(π/180)*L2[6]▶L2[6];
(π/180)*L3[6]▶L3[6];
(π/180)*L4[6]▶L4[6];

// ε : Obliquity (Laskar)
(L4[1]+L1[1])/2▶t;
(t-2451545)/3652500▶T;
84381.448▶ε;
ε-4680.93*T▶ε;
ε-1.55*T^2▶ε;
ε+1999.25*T^3▶ε;
ε-51.38*T^4▶ε;
ε-249.67*T^5▶ε;
ε-39.05*T^6▶ε;
ε+7.12*T^7▶ε;
ε+27.87*T^8▶ε;
ε+5.79*T^9▶ε;
ε+2.45*T^10▶ε;
π*ε/648000▶ε;

// GCC position of the sun
// X⒤:12, Y⒤:13, Z⒤:14
−L1[2]▶L1[12];
−L1[3]*COS(ε)+L1[4]*SIN(ε)▶L1[13];
−L1[3]*SIN(ε)-L1[4]*COS(ε)▶L1[14];
−L2[2]▶L2[12];
−L2[3]*COS(ε)+L2[4]*SIN(ε)▶L2[13];
−L2[3]*SIN(ε)-L2[4]*COS(ε)▶L2[14];
−L3[2]▶L3[12];
−L3[3]*COS(ε)+L3[4]*SIN(ε)▶L3[13];
−L3[3]*SIN(ε)-L3[4]*COS(ε)▶L3[14];
−L4[2]▶L4[12];
−L4[3]*COS(ε)+L4[4]*SIN(ε)▶L4[13];
−L4[3]*SIN(ε)-L4[4]*COS(ε)▶L4[14];

// PRINT("X⒤ Y⒤ Z⒤");
// PRINT(STRING(L1[12])+" "+STRING(L1[13])+" "+STRING(L1[14]));
// PRINT(STRING(L2[12])+" "+STRING(L2[13])+" "+STRING(L2[14]));
// PRINT(STRING(L3[12])+" "+STRING(L3[13])+" "+STRING(L3[14]));
// PRINT(STRING(L4[12])+" "+STRING(L4[13])+" "+STRING(L4[14]));
// PRINT(" ");

// a⒤:7, b⒤:8, c⒤:9
COS(L1[5])*COS(L1[6])▶L1[7];
SIN(L1[5])*COS(L1[6])▶L1[8];
SIN(L1[6])▶L1[9];
COS(L2[5])*COS(L2[6])▶L2[7];
SIN(L2[5])*COS(L2[6])▶L2[8];
SIN(L2[6])▶L2[9];
COS(L3[5])*COS(L3[6])▶L3[7];
SIN(L3[5])*COS(L3[6])▶L3[8];
SIN(L3[6])▶L3[9];
COS(L4[5])*COS(L4[6])▶L4[7];
SIN(L4[5])*COS(L4[6])▶L4[8];
SIN(L4[6])▶L4[9];

// PRINT("a⒤ b⒤ c⒤");
// PRINT(STRING(L1[7])+" "+STRING(L1[8])+" "+STRING(L1[9]));
// PRINT(STRING(L2[7])+" "+STRING(L2[8])+" "+STRING(L2[9]));
// PRINT(STRING(L3[7])+" "+STRING(L3[8])+" "+STRING(L3[9]));
// PRINT(STRING(L4[7])+" "+STRING(L4[8])+" "+STRING(L4[9]));
// PRINT(" ");

// 2R cos θ : 10
−2*(L1[7]*L1[12]+L1[8]*L1[13]+L1[9]*L1[14])▶L1[10];
−2*(L4[7]*L4[12]+L4[8]*L4[13]+L4[9]*L4[14])▶L4[10];

// PRINT("2R cos θ, for t₁ & t₄");
// PRINT(STRING(L1[10])+" "+STRING(L4[10]));
// PRINT(" ");

// R² : 11
L1[12]^2+L1[13]^2+L1[14]^2▶L1[11];
L4[12]^2+L4[13]^2+L4[14]^2▶L4[11];

// PRINT("R², for t₁ & t₄");
// PRINT(STRING(L1[11])+" "+STRING(L4[11]));
// PRINT(" ");

// L5: τ₁ τ₄ τ₂ τ₁' τ₄'
κ*(L4[1]-L2[1])▶L5[1];
κ*(L2[1]-L1[1])▶L5[2];
κ*(L4[1]-L1[1])▶L5[3];
κ*(L4[1]-L3[1])▶L5[4];
κ*(L3[1]-L1[1])▶L5[5];

// PRINT("τ₁ "+STRING(L5[1])+" τ₄ "+STRING(L5[2]));
// PRINT("τ₂ "+STRING(L5[3]));
// PRINT("τ₁' "+STRING(L5[4])+" τ₄' "+STRING(L5[5]));
// PRINT(" ");

L2[7]*L4[8]-L2[8]*L4[7]▶Δ;
L3[7]*L4[8]-L3[8]*L4[7]▶δ;

// PRINT("Δ "+STRING(Δ)+" Δ' "+STRING(δ));
// PRINT(" ");

(L1[7]*L2[8]-L1[8]*L2[7])/Δ▶A;
(L2[7]*L1[13]-L2[8]*L1[12])/Δ▶B;
(L2[8]*L2[12]-L2[7]*L2[13])/Δ▶C;
(L2[7]*L4[13]-L2[8]*L4[12])/Δ▶D;

(L1[7]*L3[8]-L1[8]*L3[7])/δ▶a;
(L3[7]*L1[13]-L3[8]*L1[12])/δ▶b;
(L3[8]*L3[12]-L3[7]*L3[13])/δ▶c;
(L3[7]*L4[13]-L3[8]*L4[12])/δ▶d;

// PRINT("A "+STRING(A)+" A' "+STRING(a));
// PRINT("B "+STRING(B)+" B' "+STRING(b));
// PRINT("C "+STRING(C)+" C' "+STRING(c));
// PRINT("D "+STRING(D)+" D' "+STRING(d));

L5[1]/L5[2]▶E;
L5[4]/L5[5]▶e;
(4/3)*L5[1]*L5[3]▶F;
(4/3)*L5[4]*L5[3]▶f;
A*E▶G;
a*e▶g;
F*(A-G)▶H;
f*(a-g)▶h;
4*A*L5[1]^2▶I;
4*a*L5[4]^2▶i;
E*(B+C)+C+D▶K;
e*(b+c)+c+d▶k;
F*(B-C+D-K)▶L;
f*(b-c+d-k)▶l;
4*(B*L5[1]^2+L5[1]*L5[2]*C)▶M;
4*(b*L5[4]^2+L5[4]*L5[5]*c)▶m;

// PRINT("E "+STRING(E)+" E' "+STRING(e));
// PRINT("F "+STRING(F)+" F' "+STRING(f));
// PRINT("G "+STRING(G)+" G' "+STRING(g));
// PRINT("H "+STRING(H)+" H' "+STRING(h));
// PRINT("I "+STRING(I)+" I' "+STRING(i));
// PRINT("K "+STRING(K)+" K' "+STRING(k));
// PRINT("L "+STRING(L)+" L' "+STRING(l));
// PRINT("M "+STRING(M)+" M' "+STRING(m));
// PRINT(" ");

// Initial guesses for r₁ & r₄
2.75▶L1[15];
2.75▶L4[15];

PRINT("r₁ "+STRING(L1[15])+" (initial guess)");
PRINT("r₄ "+STRING(L4[15])+" (initial guess)");
PRINT(" ");
PRINT("Successive approximations");

9ᴇ99▶O;
L1[15]+L4[15]▶N;

WHILE ABS(N-O)/N>1ᴇ−9 DO

(L1[15]+L4[15])^(−3)▶ξ;
(L4[15]-L1[15])/(L1[15]+L4[15])▶η;

// PRINT("ξ "+STRING(ξ)+" η "+STRING(η));

G+H*ξ+I*ξ*η▶P;
K+L*ξ+M*ξ*η▶Q;
g+h*ξ+i*ξ*η▶p;
k+l*ξ+m*ξ*η▶q;

// PRINT("P "+STRING(P)+" Q "+STRING(Q));
// PRINT("P' "+STRING(p)+" Q' "+STRING(q));

// Find ρ₁ & ρ₄
(q-Q)/(P-p)▶L1[16];
P*L1[16]+Q▶L4[16];

// Next approximation for ρ₁ & ρ₄
// PRINT("ρ₁ "+STRING(L1[16])+" ρ₄ "+STRING(L4[16]));

√(L1[11]+L1[10]*L1[16]+L1[16]^2)▶L1[15];
√(L4[11]+L4[10]*L4[16]+L4[16]^2)▶L4[15];

// Next approximation for r₁ & r₄
PRINT("r₁ "+STRING(L1[15])+" r₄ "+STRING(L4[15]));

N▶O;
L1[15]+L4[15]▶N;

END;

PRINT(" ");
PRINT("Heliocentric distances in AU at t₁ & t₄");
PRINT("r₁ "+STRING(L1[15]));
PRINT("r₄ "+STRING(L4[15]));
PRINT(" ");
PRINT("Geocentric distances in AU at t₁ & t₄");
PRINT("ρ₁ "+STRING(L1[16]));
PRINT("ρ₄ "+STRING(L4[16]));
PRINT(" ");

// L1&L4[17,18,19] : x,y,z
L1[7]*L1[16]-L1[12]▶L1[17];
L1[8]*L1[16]-L1[13]▶L1[18];
L1[9]*L1[16]-L1[14]▶L1[19];
L4[7]*L4[16]-L4[12]▶L4[17];
L4[8]*L4[16]-L4[13]▶L4[18];
L4[9]*L4[16]-L4[14]▶L4[19];

PRINT("HEC positions in AU at t₁ & t₄");
PRINT("x₁ "+STRING(L1[17]));
PRINT("y₁ "+STRING(L1[18]));
PRINT("z₁ "+STRING(L1[19]));
PRINT("x₄ "+STRING(L4[17]));
PRINT("y₄ "+STRING(L4[18]));
PRINT("z₄ "+STRING(L4[19]));
PRINT(" ");

// Reciprocal speed of light, days/AU
0.00577551833▶α;

// G : Solar gravitational parameter
1.3271244ᴇ20▶G;

// U : convert AU to meters
1.4959787ᴇ11▶U;

// β : converts AU/day to m/s
1731456.837▶β;

// 20 : Aberration correction: t₁°,t₄°
L1[1]-α*L1[16]▶L1[20];
L4[1]-α*L4[16]▶L4[20];

PRINT("Aberration corrections to time");
PRINT("Aρ₁ "+STRING(α*L1[16])+" days");
PRINT("Aρ₄ "+STRING(α*L4[16])+" days");
PRINT(" ");

// Nominal time of state vector
(L4[20]+L1[20])/2▶t;

PRINT("Epoch of state vector & obliquity");
PRINT("t₀ "+STRING(t)+" JD");
PRINT("ε₀ "+STRING(ε)+" radians");
PRINT(" ");

// State vector, celestial
(L1[15]+L4[15])/2▶r;
(L4[17]+L1[17])/2▶X;
(L4[18]+L1[18])/2▶Y;
(L4[19]+L1[19])/2▶Z;
√(X^2+Y^2+Z^2)▶R;
(r/R)*U*X▶X;
(r/R)*U*Y▶Y;
(r/R)*U*Z▶Z;
√((L4[17]-X/U)^2+(L4[18]-Y/U)^2+(L4[19]-Z/U)^2)▶S;
√((X/U-L1[17])^2+(Y/U-L1[18])^2+(Z/U-L1[19])^2)▶s;
s+S▶s;
√((L4[17]-L1[17])^2+(L4[18]-L1[18])^2+(L4[19]-L1[19])^2)▶S;
β*(s/S)*(L4[17]-L1[17])/(L4[20]-L1[20])▶I;
β*(s/S)*(L4[18]-L1[18])/(L4[20]-L1[20])▶J;
β*(s/S)*(L4[19]-L1[19])/(L4[20]-L1[20])▶K;

// State vector, ecliptic
X▶x;
Y*COS(ε)+Z*SIN(ε)▶y;
−Y*SIN(ε)+Z*COS(ε)▶z;
I▶i;
J*COS(ε)+K*SIN(ε)▶j;
−J*SIN(ε)+K*COS(ε)▶k;

PRINT("HEC state vector");
PRINT("x₀ "+STRING(x/U)+" AU");
PRINT("y₀ "+STRING(y/U)+" AU");
PRINT("z₀ "+STRING(z/U)+" AU");
PRINT("Vx₀ "+STRING(i)+" m/s");
PRINT("Vy₀ "+STRING(j)+" m/s");
PRINT("Vz₀ "+STRING(k)+" m/s");
PRINT(" ");

√(x^2+y^2+z^2)▶r;
√(i^2+j^2+k^2)▶V;

PRINT("Heliocentric distance");
PRINT("r "+STRING(r/U)+" AU");
PRINT(" ");
PRINT("Sun−relative speed");
PRINT("v "+STRING(V)+" m/s");
PRINT(" ");

(2/r-V^2/G)⁻¹/U▶a;

// L6[1,2,3] : r※V
y*k-z*j▶L6[1];
z*i-x*k▶L6[2];
x*j-y*i▶L6[3];

// H : Angular momentum per unit mass
√(L6[1]^2+L6[2]^2+L6[3]^2)▶H;

// e : Eccentricity
√(1-H*H/(a*U*G))▶e;

// N : Inclination
(180/π)*ACOS(L6[3]/H)▶N;

// Ω : Longitude of ascending node
(180/π)*ATAN(−L6[1]/L6[2])▶Ω;
IF −L6[2]<0 THEN
Ω+180▶Ω
END;
IF −L6[2]>0 AND L6[1]<0 THEN
Ω+360▶Ω
END;

// θ : True anomaly
H^2/(r*G)-1▶Ĉ;
H*(x*i+y*j+z*k)/(r*G)▶Ŝ;
(180/π)*ATAN(Ŝ/Ĉ)▶θ;
IF Ĉ<0 THEN
θ+180▶θ
END;
IF Ĉ>0 AND Ŝ<0 THEN
θ+360▶θ
END;

PRINT("True anomaly "+STRING(θ)+"°");

// M : Sum θ+ω
(x*COS(Ω*π/180)+y*SIN(Ω*π/180))/r▶Ĉ;
z/(r*SIN(N*π/180))▶Ŝ;
(180/π)*ATAN(Ŝ/Ĉ)▶M;
IF Ĉ<0 THEN
M+180▶M
END;
IF Ĉ>0 AND Ŝ<0 THEN
M+360▶M
END;

// ω : Argument of perihelion
M-θ▶ω;
IF ω<0 THEN
ω+360▶ω
END;

// Eccentric anomaly
1-r/(a*U)▶Ĉ;
(x*i+y*j+z*k)/√(a*U*G)▶Ŝ;
(180/π)*ATAN(Ŝ/Ĉ)▶Q;
IF Ĉ<0 THEN
Q+180▶Q
END;
IF Ĉ>0 AND Ŝ<0 THEN
Q+360▶Q
END;

PRINT("Ecc. anomaly "+STRING(Q)+"°");

// Mean anomaly
Q*π/180▶Q;
Q-e*SIN(Q)▶Q;
Q*180/π▶Q;

PRINT("Mean anomaly "+STRING(Q)+"°");
PRINT(" ");

// Period of orbit
365.256898326*a^1.5▶P;

PRINT("Period of orbit "+STRING(P)+" days");
PRINT(" ");

// T : Time of perihelion passage
t-P*Q/360▶T;

PRINT("Orbital elements");
PRINT("a "+STRING(a)+" AU");
PRINT("e "+STRING(e));
PRINT("i "+STRING(N)+"°");
PRINT("Ω "+STRING(Ω)+"°");
PRINT("ω "+STRING(ω)+"°");
PRINT("T "+STRING(T)+" JD");
PRINT("T+P "+STRING(T+P)+" JD");

END;


Complete Output from a Test Run with Data for Ceres

ORBIT4 by David Sims
Method of Gauss with four observed positions
to find the Keplerian orbital elements.
User provides input by adjusting inline data.

r₁ 2.75 (initial guess)
r₄ 2.75 (initial guess)

Successive approximations
r₁ 2.90652064 r₄ 2.92071388
r₁ 2.93008666 r₄ 2.94292188
r₁ 2.93307059 r₄ 2.94572742
r₁ 2.93344165 r₄ 2.94607627
r₁ 2.93348769 r₄ 2.94611956
r₁ 2.9334934 r₄ 2.94612492
r₁ 2.93349411 r₄ 2.94612559
r₁ 2.9334942 r₄ 2.94612567
r₁ 2.93349421 r₄ 2.94612568
r₁ 2.93349421 r₄ 2.94612568

Heliocentric distances in AU at t₁ & t₄
r₁ 2.93349421
r₄ 2.94612568

Geocentric distances in AU at t₁ & t₄
ρ₁ 2.00460681
ρ₄ 1.94781669

HEC positions in AU at t₁ & t₄
x₁ 1.33687069
y₁ −2.2463475
z₁ −1.33119794
x₄ 1.58994315
y₄ −2.10289848
z₄ −1.31512559

Aberration corrections to time
Aρ₁ 0.011577643 days
Aρ₄ 0.011249651 days

Epoch of state vector & obliquity
t₀ 2457219.61 JD
ε₀ 0.409057547 radians

HEC state vector
x₀ 1.46520344 AU
y₀ −2.52458426 AU
z₀ −0.349479243 AU
Vx₀ 14610.4367 m/s
Vy₀ 7967.42879 m/s
Vz₀ −2442.63758 m/s

Heliocentric distance
r 2.93980995 AU

Sun−relative speed
v 16819.9661 m/s

True anomaly 147.669798°
Ecc. anomaly 145.259666°
Mean anomaly 142.77737°

Period of orbit 1681.12408 days

Orbital elements
a 2.76694735 AU
e 0.076026341
i 10.5918141°
Ω 80.3183813°
ω 72.6265868°
T 2456552.87 JD
T+P 2458234 JD

For comparison, the orbit of asteroid 1 Ceres for 2457219.61 JD from JPL Horizons.

a 2.768008676 AU
e 0.075773357
i 10.59221734°
Ω 80.32683297°
ω 72.66267214°
T 2456552.644 JD
P 1682.091 days


Input and Output from a Test Run for Vesta

// Data for time 1
L1:=
{2458319.625,
0.46327126,
-0.90443214,
0.00003691,
HMS→(17°27'35.43"),
HMS→(-21°57'42.0")};

// Data for time 2
L2:=
{2458329.625,
0.60636977,
-0.81427714,
0.0000325,
HMS→(17°24'52.39"),
HMS→(-22°38'47.9")};

// Data for time 3
L3:=
{2458340.625,
0.74387086,
-0.6886877,
0.00002659,
HMS→(17°25'53.52"),
HMS→(-23°21'49.6")};

// Data for time 4
L4:=
{2458350.625,
0.84685421,
-0.55390808,
0.00002041,
HMS→(17°30'18.90"),
HMS→(-23°58'18.1")};

The orbital elements calculated are:

a 2.359195134 AU
e 0.0885562327
i 7.146912208°
Ω 103.797078°
ω 149.3285751°
T 2458243.236 JD

For comparison, JPL's orbit is:

a 2.3616071252 AU
e 0.088874634
i 7.141583463°
Ω 103.8114586°
ω 150.719275°
T 2458247.557 JD


Input and Output from a Test Run for Mars

// Data for time 1
L1:=
{2458309.7,
0.308344662,
-0.968753701,
0.0000402562,
HMS→(20°47'23.26"),
HMS→(-23°45'36.7")};

// Data for time 2
L2:=
{2458314.7,
0.387736428,
-0.93960836,
0.0000387037,
HMS→(20°43'46.98"),
HMS→(-24°18'18.9")};

// Data for time 3
L3:=
{2458319.7,
0.464398349,
-0.90384802,
0.000036878,
HMS→(20°39'10.17"),
HMS→(-24°51'14.8")};

// Data for time 4
L4:=
{2458324.7,
0.537787852,
-0.861718947,
0.0000347918,
HMS→(20°33'49.66"),
HMS→(-25°22'18.2")};

The orbital elements calculated are:

a 1.519845257 AU
e 0.0919892461
i 1.83846542°
Ω 49.55960152°
ω 285.6210229°
T 2458376.367 JD

For comparison, JPL's orbit is:

a 1.5236555436 AU
e 0.0933250338
i 1.8481512723°
Ω 49.506875204°
ω 286.65369589°
T 2458377.999 JD