Vanguard News Network
VNN Media
VNN Digital Library
VNN Reader Mail
VNN Broadcasts

Old September 10th, 2018 #1
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default A new HP Prime program to determine asteroid orbits

Hello, fellow white nationalists! Those of you who are interested in astronomy, or in the math pertaining to the motions of planets and asteroids, might find my new calculator program interesting, useful, and maybe both.

I use orbit determination by the method of Gauss in my program ORBIT4. There is no need to have a priori range information. It is available for free download from

https://my.cloudme.com/jenab6/ORBIT4

The orbits that the program is designed to treat are elliptical (two-body) orbits around the sun. (It wasn't designed for objects in orbit around Earth.)

ORBIT4 is written in Prime Programming Language (PPL) for the HP Prime Calculator or its emulators. The input, in this case for the asteroid 1 Ceres, is in the program (inline code), and looks like this:

// 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″)};

The first number (after the open curly bracket) is the time of observation in Julian Date. It should be accurate to about 0.0001 days. The observations should be separated in time by somewhere between 0.5% to 1.0% of the object's orbital period. The observation interval should be reasonably near the opposition of the asteroid with the Sun, but it should not span an apside of the asteroid's orbit.

The 2nd number is the X component of the Earth's position in heliocentric ecliptic coordinates. The 3rd number is the Y component. The 4th number is the Z component. These distances are in astronomical units, should be accurate to eight significant figures, and can be obtained from JPL Horizons.

The 5th number, inside the HMS operand parentheses, is the asteroid's geocentric right ascension in HH°MM'SS.SS" format. Right ascension should be accurate to 0.01 seconds of time.

The 6th number, also inside HMS operand parantheses, is the asteroid's geocentric declination in degrees, arcminutes, arcseconds format. Declination should be accurate to 0.1 arcsec.

If observational data aren't available, then test data for the RA & DEC of known asteroids or planets can be obtained from JPL Horizons.

The output for the data shown in the PPL code above follows:

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 [for Ceres]
a 2.76694735 AU
e 0.076026341
i 10.5918141°
Ω 80.3183813°
ω 72.6265868°
T 2456552.87 JD
T+P 2458234 JD

For comparison, here's the official orbit from NASA/JPL for Ceres.

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

You can examine the source code for ORBIT4 on my LiveJournal website. And here's the download link for the executable code for ORBIT4.

To run ORBIT4, you can buy an HP Prime calculator and install the program on it. For example, I bought two HP Prime calculators from this eBay seller, so I know that he does deliver them by mail. You can choose a different seller if you have another preference.

But you don't need to own the physical calculator, if you download two applications for the Windows PC from the Hewlett Packard company.

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

The Connectivity Kit does the file transfers...
...from a physical HP Prime calculator to the computer desktop, via a special USB cable
...from the computer desktop to a physical HP Prime calculator, via a special USB cable
...from an HP Prime emulator running on your computer to the computer desktop
...from your computer desktop to an HP Prime emulator running on your computer

The Virtual Calculator Emulator runs a simulated calculator on your computer. It's like an extreme upgrade to the normal Windows calculator. One that you can write long PPL programs on and run them. (PPL = Prime Programming Language.) You can also install someone else's programs in your emulator, once you have the executable code.

Last edited by Jerry Abbott; September 10th, 2018 at 10:27 PM.
 
Old September 23rd, 2018 #2
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default The math from the program (in plain algebra)

Determining an Orbit from Four Observations in Right Ascension and Declination

September 19th, 23:31

.
The initial Data consist of

• The time of observation in Julian date
Preferably accurate to at least 0.00001 days (or to the nearest second)

• The X component of Earth's position in heliocentric ecliptic coordinates
In astronomical units, accurate to at least six significant figures

• The Y component of Earth's position in heliocentric ecliptic coordinates
In astronomical units, accurate to at least six significant figures

• The Z component of Earth's position in heliocentric ecliptic coordinates
In astronomical units, accurate to at least six significant figures

• The asteroid's observed right ascension
Accurate to 0.01 seconds

• The asteroid's observed declination
Accurate to 0.1 arcsec


for four observations of the asteroid. Each observation should be separated from the next in time by 0.5% to 1% of the asteroid's orbital period. The observations should be made near an opposition of the asteroid with the sun, relative to Earth, but they should not span an apside of the asteroid's orbit.


Initial data

For an example case involving the asteroid Ceres.

Data for time 1
t₁ = 2457204.62500
X₁' = 0.155228396
Y₁' = −1.004732775
Z₁' = 0.00003295786
α₁ = 20:46:57.02s
δ₁ = −27°41′33.9″

Data for time 2
t₂ = 2457214.62500
X₂' = 0.319493277
Y₂' = −0.965116604
Z₂' = 0.0000311269
α₂ = 20:39:57.10
δ₂ = −28°47′21.5″

Data for time 3
t₃ = 2457224.62500
X₃' = 0.4747795623
Y₃' = −0.8983801739
Z₃' = 0.00002841127
α₃ = 20:31:22.81
δ₃ = −29°49′22.7″

Data for time 4
t₄ = 2457234.62500
X₄' = 0.616702829
Y₄' = −0.8063620175
Z₄' = 0.00002486325
α₄ = 20:22:06.57
δ₄ = −30°41′57.3″



Handy constants

k : Mean motion of Earth
k = 0.01720209895 rad/day

γ : Reciprocal speed of light
γ = 0.00577551833 days/AU

GM : Solar gravitational parameter
GM = 1.3271244ᴇ20 m³ kg⁻¹ sec⁻²

U : convert distance from AU to meters
U = 1.4959787ᴇ11 m/AU

β : converts speed from AU/day to m/s
β = 1731456.837 AU sec m⁻¹ day⁻¹



Convert the celestial angles to radians

In the example case,

α₁ = 5.44084723 radians δ₁ = −0.483329666 radians
α₂ = 5.41030979 radians δ₂ = −0.502468171 radians
α₃ = 5.37290956 radians δ₃ = −0.520509058 radians
α₄ = 5.33245865 radians δ₄ = −0.53580299 radians



Find the geocentric celestial unit vectors toward the asteroid

a₁ = cos α₁ cos δ₁ ⇨ example ⇨ +0.589463371
b₁ = sin α₁ cos δ₁ ⇨ example ⇨ −0.660726081
c₁ = sin δ₁ ⇨ example ⇨ −0.464730008

a₂ = cos α₂ cos δ₂ ⇨ example ⇨ +0.563195268
b₂ = sin α₂ cos δ₂ ⇨ example ⇨ −0.671477524
c₂ = sin δ₂ ⇨ example ⇨ −0.4815901

a₃ = cos α₃ cos δ₃ ⇨ example ⇨ +0.532276132
b₃ = sin α₃ cos δ₃ ⇨ example ⇨ −0.685093499
c₃ = sin δ₃ ⇨ example ⇨ −0.497321844

a₄ = cos α₄ cos δ₄ ⇨ example ⇨ +0.49965704
b₄ = sin α₄ cos δ₄ ⇨ example ⇨ −0.69978587
c₄ = sin δ₄ ⇨ example ⇨ −0.510531662



Find the obliquity of the ecliptic, ε, in radians

t₀ = (t₁+t₄)/2 ⇨ example ⇨ 2457219.625 JD

τ = (t₀−2451545)/3652500 ⇨ example ⇨ 0.001553627652

ε = (π/648000){84381.448 − 4680.93 τ − 1.55 τ² + 1999.25 τ³ − 51.38 τ⁴
− 249.67 τ⁵ − 39.05 τ⁶ + 7.12 τ⁷ + 27.87 τ⁸ + 5.79 τ⁹ + 2.45 τ¹⁰}

In the example, ε = 0.409057547 radians.



Find the position of the sun in geocentric celestial coordinates

X₁ = −X₁' ⇨ example ⇨ −0.155228396
Y₁ = −Y₁' cos ε + Z₁' sin ε ⇨ example ⇨ +0.921851498
Z₁ = −Y₁' sin ε − Z₁' cos ε ⇨ example ⇨ +0.399597005

X₂ = −X₂' ⇨ example ⇨ −0.319493277
Y₂ = −Y₂' cos ε + Z₂' sin ε ⇨ example ⇨ +0.885503087
Z₂ = −Y₂' sin ε − Z₂' cos ε ⇨ example ⇨ +0.383841559

X₃ = −X₃' ⇨ example ⇨ −0.474779562
Y₃ = −Y₃' cos ε + Z₃' sin ε ⇨ example ⇨ +0.824271594
Z₃ = −Y₃' sin ε − Z₃' cos ε ⇨ example ⇨ +0.357299982

X₄ = −X₄' ⇨ example ⇨ −0.616702829
Y₄ = −Y₄' cos ε + Z₄' sin ε ⇨ example ⇨ +0.739843884
Z₄ = −Y₄' sin ε − Z₄' cos ε ⇨ example ⇨ +0.320703493



Find some intermediate quantities

2R₁ cos ζ₁ = −2(a₁X₁ + b₁Y₁ + c₁Z₁) ⇨ example ⇨ 1.772595
R₁² = X₁² + Y₁² + Z₁² ⇨ example ⇨ 1.97920299

2R₄ cos ζ₄ = −2(a₄X₄ + b₄Y₄ + c₄Z₄) ⇨ example ⇨ 1.03358381
R₄² = X₄² + Y₄² + Z₄² ⇨ example ⇨ 1.03054208



Time differences

τ₁ = k(t₄−t₂) ⇨ example ⇨ 0.344041979
τ₄ = k(t₂−t₁) ⇨ example ⇨ 0.17202099
τ₂ = k(t₄−t₁) ⇨ example ⇨ 0.516062969
τ₁'= k(t₄−t₃) ⇨ example ⇨ 0.17202099
τ₄'= k(t₃−t₁) ⇨ example ⇨ 0.344041979



Determine some determinants

Δ = a₂b₄−b₂a₄ ⇨ example ⇨ −0.058607619
Δ'= a₃b₄−b₃a₄ ⇨ example ⇨ −0.030167527



Find some constants

A = (a₁b₂−b₁a₂)/Δ A' = (a₁b₃−b₁a₃)/Δ'
B = (a₂Y₁−b₂X₁)/Δ B' = (a₃Y₁−b₃X₁)/Δ'
C = (b₂X₂−a₂Y₂)/Δ C' = (b₃X₃−a₃Y₃)/Δ'
D = (a₂Y₄−b₂X₄)/Δ D' = (a₃Y₄−b₃X₄)/Δ'
E = τ₁/τ₄ E' = τ₁'/τ₄'
F = (4/3) τ₁τ₂ F' = (4/3) τ₁'τ₂
G = AE G' = A'E'
H = F(A−G) H' = F'(A'−G')
I = 4Aτ₁² I' = 4A'(τ₁')²
K = E(B+C)+C+D K' = E'(B'+C')+C'+D'
L = F(B−C+D−K) L' = F'(B'−C'+D'−K')
M = 4[Bτ₁² + τ₁τ₄C] M' = 4[B'(τ₁')² + τ₁'τ₄'C']

In our example case,

A = 0.404275125 A' = 1.72864027
B = −7.08013785 B' = −12.7399767
C = 4.84883362 C' = 3.76138574
D = −0.043927505 D' = 0.951283093
E = 2 E' = 0.5
F = 0.236729767 F' = 0.118364883
G = 0.80855025 G' = 0.864320133
H = −0.095703956 H' = 0.102305152
I = 0.191407912 I' = 0.204610303
K = 0.342297675 K' = 0.223373365
L = −2.91537363 L' = −1.86702289
M = −2.20429551 M' = −0.617533885



Make first guess at the asteroid's heliocentric distances at time 1 and at time 4

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



Converge r₁ , r₄ , ρ₁ , ρ₄ by successive approximations

OldSum = 9e+99
NewSum = r₁ + r₄

WHILE |NewSum−OldSum|/NewSum > 1e-9 DO

ξ = 1 / (r₁ + r₄)³
η = (r₄ − r₁) / (r₁ + r₄)

P = G + Hξ + Iξη
Q = K + Lξ + Mξη
P' = G' + H'ξ + I'ξη
Q' = K' + L'ξ + M'ξη

ρ₁ = (Q'−Q)/(P−P')
ρ₄ = Pρ₁ + Q

r₁ = √(R₁² + 2ρ₁R₁ cos ζ₁ + ρ₁²)
r₄ = √(R₄² + 2ρ₄R₄ cos ζ₄ + ρ₄²)

OldSum = NewSum
NewSum = r₁ + r₄

WHILE-END


It is possible that r₁ and r₄ will not converge, but will bounce around erratically in value. This usually means that your observations weren't sufficiently close to the opposition of the asteroid with the sun, with respect to Earth. On the other hand, the reason might be that you didn't note the asteroid's right ascension and declination with enough accuracy.

However, in our test case:

r₁ = 2.75 (initial guess) r₄ = 2.75 (initial guess)
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

Converged values:

r₁ = 2.93349421 r₄ = 2.94612568
ρ₁ = 2.00460681 ρ₄ = 1.94781669



The heliocentric position vectors at time 1 and time 4

x₁ = a₁ρ₁ − X₁ ⇨ example ⇨ +1.33687069
y₁ = b₁ρ₁ − Y₁ ⇨ example ⇨ −2.2463475
z₁ = c₁ρ₁ − Z₁ ⇨ example ⇨ −1.33119794

x₄ = a₄ρ₄ − X₄ ⇨ example ⇨ +1.58994315
y₄ = b₄ρ₄ − Y₄ ⇨ example ⇨ −2.10289848
z₄ = c₄ρ₄ − Z₄ ⇨ example ⇨ −1.31512559



Correct times of observation for planetary aberration

t₁° = t₁ − γρ₁
t₄° = t₄ − γρ₄

In our example case,

−γρ₁ = −0.011577643 days
−γρ₄ = −0.011249651 days

t₁° = 2457204.61342
t₄° = 2457234.61375



Epoch of state vector

t₀ = (t₁° + t₄°)/2 ⇨ example ⇨ 2457219.61358 JD



Position part of state vector, heliocentric celestial coordinates, AU

r₀ = (r₁ + r₄)/2
x₀' = (x₁ + x₄)/2
y₀' = (y₁ + y₄)/2
z₀' = (z₁ + z₄)/2
r₀' = √{(x₀')² + (y₀')² + (z₀')²}
x₀ = x₀'(r₀/r₀') ⇨ example ⇨ +1.465203439 AU
y₀ = y₀'(r₀/r₀') ⇨ example ⇨ −2.177292621 AU
z₀ = z₀'(r₀/r₀') ⇨ example ⇨ −1.324786117 AU



Velocity part of the state vector, heliocentric celestial coordinates, m/s

S₁ = √{(x₁−x₀)² + (y₁−y₀)² + (z₁−z₀)²}
S₄ = √{(x₄−x₀)² + (y₄−y₀)² + (z₄−z₀)²}
S = S₁ + S₄
s = √{(x₄−x₁)² + (y₄−y₁)² + (z₄−z₁)²}
Vx₀ = β(S/s)(x₄−x₁)/(t₄°−t₁°) ⇨ example ⇨ +14610.4367 m/s
Vy₀ = β(S/s)(y₄−y₁)/(t₄°−t₁°) ⇨ example ⇨ +8281.6311 m/s
Vz₀ = β(S/s)(z₄−z₁)/(t₄°−t₁°) ⇨ example ⇨ +927.8930 m/s



State vector in heliocentric ecliptic coordinates, MKS

x = Ux₀ ⇨ example ⇨ +2.191913145e+11 meters
y = U(y₀ cos ε + z₀ sin ε) ⇨ example ⇨ −3.776724292e+11 meters
z = U(z₀ cos ε − y₀ sin ε) ⇨ example ⇨ −5.228135061e+10 meters
Vx = Vx₀ ⇨ example ⇨ +14610.4367 m/s
Vy = Vy₀ cos ε + Vz₀ sin ε ⇨ example ⇨ +7967.4288 m/s
Vz = Vz₀ cos ε − Vy₀ sin ε ⇨ example ⇨ −2442.6375 m/s



Heliocentric distance and sun-relative speed

r = √(x² + y² + z²) ⇨ example ⇨ 4.39789309e+11 meters
V = √(Vx² + Vy² + Vz²) ⇨ example ⇨ 16819.9661 m/s



The vector of angular momentum per unit mass

hx = y Vz − z Vy ⇨ example ⇨ +1.33906480091ᴇ15 m²/s
hy = z Vx − x Vz ⇨ example ⇨ −2.28448420525ᴇ14 m²/s
hz = x Vy − y Vx ⇨ example ⇨ +7.26435027984ᴇ15 m²/s

h = √( hx² + hy² + hz² ) ⇨ example ⇨ 7.39026848024ᴇ15 m²/s



The semimajor axis of the asteroid's orbit

a = [ U (2/r − V²/GM) ]⁻¹

Example: a = 2.766947347 AU



The eccentricity of the asteroid's orbit

e = √{ 1 − h² / (U a GM) }

Example: e = 0.07602634107



The inclination of the asteroid's orbit

i = (180/π) arccos( hz / h )

Example: i = 10.5918141°



The longitude the ascending node

Ĉ = −hy
Ŝ = hx
Ω' = (180/π) arctan(Ŝ/Ĉ)
If Ĉ>0 and Ŝ>0 then Ω = Ω'
If Ĉ<0 then Ω = Ω' + 180°
If Ĉ>0 and Ŝ<0 then Ω = Ω' + 360°

Example: Ω = 80.3183813°



The true anomaly of the asteroid in its orbit

Ĉ = h²/(r GM) − 1
Ŝ = h ( x Vx + y Vy + z Vz ) / (r GM)
θ' = (180/π) arctan(Ŝ/Ĉ)
If Ĉ>0 and Ŝ>0 then θ = θ'
If Ĉ<0 then θ = θ' + 180°
If Ĉ>0 and Ŝ<0 then θ = θ' + 360°

Example: θ = 147.6697976°



The sum of the true anomaly and the argument of perihelion

Ĉ = (x cos Ω + y sin Ω) / r
Ŝ = z / (r sin i)
φ' = (180/π) arctan(Ŝ/Ĉ)
If Ĉ>0 and Ŝ>0 then φ = φ'
If Ĉ<0 then φ = φ' + 180°
If Ĉ>0 and Ŝ<0 then φ = φ' + 360°

Example: φ = 220.2963844°

Be sure to take account of the angular units of Ω and i, in this block.



The argument of the perihelion

ω' = φ − θ
If ω'≥0 then ω = ω'
If ω'<0 then ω = ω'+360°

Example: ω = 72.6265868°



The eccentric anomaly

Ĉ = 1 − r/(U a)
Ŝ = (x Vx + y Vy + z Vz) / √(U a GM)
u' = (180/π) arctan(Ŝ/Ĉ)
If Ĉ>0 and Ŝ>0 then u = u'
If Ĉ<0 then u = u' + 180°
If Ĉ>0 and Ŝ<0 then u = u' + 360°

Example: u = 145.2596659°



The mean anomaly

M = u − e sin u

In our example, M = 142.7773704°

In this block, u must be in radians. The resulting M will also be in radians,
which you can convert to degrees by multiplying by (180/π).



The period of the asteroid's orbit

P = (365.256898326 days) a¹·⁵

Example: P = 1681.12408 days



The time of perihelion passage

T = t₀ − P M / 360°

Example: T = 2456552.87337 JD

In this block, the mean anomaly M must be in the same angular units as the
denominator, in this case degrees. If M is in radians, then the denominator
must be changed from 360° to 2π radians.



Compare this calculated orbit for Ceres with that from JPL Horizons

Orbital elements of Ceres for JD 2457219.61358

Example JPL Horizons
a 2.766947347 AU 2.768008676 AU
e 0.076026341 0.075773357
i 10.59181413° 10.59221734°
Ω 80.31838125° 80.32683297°
ω 72.62658682° 72.66267214°
T 2456552.873 JD 2456552.644 JD
P 1681.124 days 1682.091 days



Acknowledgements

To the Russian astronomer Alexander Dmitriyevich Dubyago, from whose book The Determination of Orbits about 75% of the method used on this page was taken.

To the French astronomer Jacques Laskar, who supplied the 10-degree polynomial to calculate the obliquity of the ecliptic.
 
Old September 24th, 2018 #3
James T Kirk
Banned
 
Join Date: Aug 2018
Location: Living on Cestus III
Posts: 98
Default

I've fed all your data into shipboard computers. We're now taking care of the asteroid problem...

 
Reply

Share


Thread
Display Modes


All times are GMT -5. The time now is 05:23 PM.
Page generated in 0.09344 seconds.