
September 10th, 2018  #1 
Senior Member
Join Date: Nov 2007
Posts: 942

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 (twobody) 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. 
September 23rd, 2018  #2 
Senior Member
Join Date: Nov 2007
Posts: 942

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 > 1e9 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₄ WHILEEND 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 sunrelative 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 10degree polynomial to calculate the obliquity of the ecliptic. 
September 24th, 2018  #3 
Banned
Join Date: Aug 2018
Location: Living on Cestus III
Posts: 98

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