
September 13th, 2020  #1 
Senior Member
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 996

Determining a Preliminary Orbit from Four Observations
The Determination of a Preliminary Orbit from Four Observations
Presented hereafter is a method for determining a preliminary heliocentric orbit from four geocentric directions of a sunorbiting object at four distinct times of observation. I take this method after that presented in chapter six of "The Determination of Orbits," by A.D. Dubyago. Summary of the method The initial data t₁, X⊕₁, Y⊕₁, Z⊕₁, α₁, δ₁ t₂, X⊕₂, Y⊕₂, Z⊕₂, α₂, δ₂ t₃, X⊕₃, Y⊕₃, Z⊕₃, α₃, δ₃ t₄, X⊕₄, Y⊕₄, Z⊕₄, α₄, δ₄ The times of observation, tᵢ, are given in Julian date format. The vectors [X⊕ᵢ,Y⊕ᵢ,Z⊕ᵢ] are the positions of the Earth in heliocentric ecliptic coordinates, with the components being in astronomical units. The αᵢ are the geocentric right ascensions for Ceres. The δᵢ are the geocentric declinations for Ceres. The time intervals should be about 0.5% to 1% of the object's period (estimated as 8 to 16 days for a main belt asteroid), should be near opposition with the sun, but should NOT span an apside of the object's orbit. The precision in right ascension should be 0.01 seconds or better, and the precision in declination should be 0.1 arcseconds or better. The Earth's orbital mean motion, κ = 0.01720209895 radians/day We will use a single value for the obliquity of the ecliptic to transform all four of the observation angles from celestial coordinates to ecliptic coordinates. First, we find the middle of the observation time window in units of 10000 years since 1 January 2000, and then we use that number to find the obliquity using the 10degree polynomial fit of J. Laskar. t = (t₁ + t₄)/2 T = (t − 2451545)/3652500 The obliquity in seconds of arc is ε" = 84381.448 − 4680.93 T − 1.55 T² + 1999.25 T³ − 51.38 T⁴ − 249.67 T⁵ − 39.05 T⁶ + 7.12 T⁷ + 27.87 T⁸ + 5.79 T⁹ + 2.45 T¹⁰ The obliquity in radians is ε = (π / 648000) ε" Although Earth's axial tilt (i.e. the obliquity of the ecliptic) does change over time, it changes so slowly that the difference will almost always be negligible across the t₁ to t₄ time window. However, in the rare case when this isn't true, separate evaluations of the obliquity will have to be made for each time of observation. The geocentric positions of the sun in celestial coordinates are (for i = 1 to 4) Xᵢ = −X⊕ᵢ Yᵢ = −Y⊕ᵢ cos ε + Z⊕ᵢ sin ε Zᵢ = −Y⊕ᵢ sin ε − Z⊕ᵢ cos ε The geocentric unit vectors in the direction of the target object, in celestial coordinates, are (for i = 1 to 4) aᵢ = cos αᵢ cos δᵢ bᵢ = sin αᵢ cos δᵢ cᵢ = sin δᵢ The squares of the SunEarth distances at times t₁ and t₄ are R₁² = X₁² + Y₁² + Z₁² R₄² = X₄² + Y₄² + Z₄² The values of 2Rᵢ cos θᵢ at times t₁ and t₄ are 2R₁ cos θ₁ = −2 ( a₁X₁ + b₁Y₁ + c₁Z₁ ) 2R₄ cos θ₄ = −2 ( a₄X₄ + b₄Y₄ + c₄Z₄ ) where θ is the supplementary angle to the sunEarthobject angle. We find some time differences: τ₁ = κ (t₄−t₂) τ₂ = κ (t₂−t₁) τ₃ = κ (t₄−t₁) τ₄ = κ (t₄−t₃) τ₅ = κ (t₃−t₁) We find this pair of determinants: Φ = a₂b₄−b₂a₄ φ = a₃b₄−b₃a₄ We calculate some intermediate quantities: A = ( a₁b₂ − b₁a₂ ) / Φ B = ( a₂Y₁ − b₂X₁ ) / Φ C = ( b₂X₂ − a₂Y₂ ) / Φ D = ( a₂Y₄ − b₂X₄ ) / Φ A' = ( a₁b₃ − b₁a₃ ) / φ B' = ( a₃Y₁ − b₃X₁ ) / φ C' = ( b₃X₃ − a₃Y₃ ) / φ D' = ( a₃Y₄ − b₃X₄ ) / φ And then more intermediate quantities: E = τ₁/τ₂ F = (4/3) τ₁τ₃ G = AE H = F(A−G) I = 4Aτ₁² K = E (B + C) + C + D L = F (B − C + D − K) M = 4 (Bτ₁² + τ₁τ₂C) E' = τ₄/τ₅ F' = (4/3) τ₄τ₃ G' = A'E' H' = F'(A'−G') I' = 4A'τ₄² K' = E' (B' + C') + C' + D' L' = F' (B' − C' + D' − K') M' = 4 (B'τ₄² + τ₄τ₅C') Make initial guesses for the sunobject distance r₁ at time t₁ and for the sunobject distance r₄ at time t₄. For main belt asteroids, a reasonable initial guess for both times is 2.75 AU. Then use the loop below to converge, by successive approximations, to the true sunobject distances, r, and for the Earthobject distances, ρ, distances at times t₁° and t₄° (i.e. the times for the first and fourth observations, corrected for the speed of light travel time). O = 9.999e+99 N = r₁ + r₄ while N−O/N > 1ᴇ11 do ξ = (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₁² + (2R₁cos θ₁)ρ₁ + ρ₁²] r₄ = √[R₄² + (2R₄cos θ₄)ρ₄ + ρ₄²] O = N N = r₁ + r₄ endwhile The positions of the object at times t₁ and t₄ in geocentric celestial coordinates are x₁ = a₁ρ₁ − X₁ y₁ = b₁ρ₁ − Y₁ z₁ = c₁ρ₁ − Z₁ x₄ = a₄ρ₄ − X₄ y₄ = b₄ρ₄ − Y₄ z₄ = c₄ρ₄ − Z₄ The reciprocal of the speed of light ç = 0.00577551833 days/AU The sun's gravitational parameter μ = 1.32712440018ᴇ20 m³ sec⁻² The conversion factor from AU to meters is U = 1.495978707ᴇ11 The conversion factor from AU/day to m/sec is β = 1731456.8368 Correcting times of observation for planetary abberation. t₁° = t₁ − çρ₁ t₄° = t₄ − çρ₄ The nominal time associated with the forthcoming state vector is t₀ = ½ (t₁° + t₄°) The nominal heliocentric distance of the object at time t₀ is r₀ = ½ (r₁ + r₄) Find the heliocentric position vector [x',y',z'] for the object at time t₀ in celestial coordinates. x" = ½ (x₁ + x₄) y" = ½ (y₁ + y₄) z" = ½ (z₁ + z₄) r" = √[(x")²+(y")²+(z")²] x' = (r₀/r") U x" y' = (r₀/r") U y" z' = (r₀/r") U z" Find the sunrelative velocity vector for the object at time t₀ in celestial coordinates. S = √[ (x₄ − x'/U)² + (y₄ − y'/U)² + (z₄ − z'/U)² ] s = √[ (x'/U − x₁)² + (y'/U − y₁)² + (z'/U − z₁)² ] Ψ = S + s ψ = √[(x₄−x₁)² + (y₄−y₁)² + (z₄−z₁)²] Vx' = β (Ψ/ψ) (x₄−x₁) / (t₄°−t₁°) Vy' = β (Ψ/ψ) (y₄−y₁) / (t₄°−t₁°) Vz' = β (Ψ/ψ) (z₄−z₁) / (t₄°−t₁°) The object's position in heliocentric ecliptic coordinates x₀ = x' y₀ = y' cos ε + z' sin ε z₀ = −y' sin ε + z' cos ε The object's sunrelative velocity in ecliptic coordinates Vx₀ = Vx' Vy₀ = Vy' cos ε + Vz' sin ε Vz₀ = −Vy' sin ε + Vz' cos ε The object's speed relative to the sun V₀ = √[(Vx₀)² + (Vy₀)² + (Vz₀)²] The semimajor axis of the object's orbit, in AU a = (2/r₀ − V₀²/μ)⁻¹ / U The angular momentum per unit mass in the object's orbit hx = y₀ Vz₀ − z₀ Vy₀ hy = z₀ Vx₀ − x₀ Vz₀ hz = x₀ Vy₀ − y₀ Vx₀ h = √[(hx)² + (hy)² + (hz)²] The eccentricity of the object's orbit e = √[1 − h²/(aμU)] The inclination of the object's orbit i = arccos(hz/h) The longitude of the ascending node of the object's orbit Ω' = arctan(−hx/hy) if hy>0 then Ω = Ω' + π If hy<0 and hx<0 then Ω = Ω' + 2π The true anomaly at time t₀ sin θ₀ = h ( x₀ Vx₀ + y₀ Vy₀ + z₀ Vz₀ ) / (r₀μ) cos θ₀ = h²/(r₀μ) − 1 θ₀' = arctan( sin θ₀ / cos θ₀ ) If cos θ₀ < 0 then θ₀ = θ₀' + π If cos θ₀ > 0 and sin θ₀ < 0 then θ₀ = θ₀' + 2π The sum of the true anomaly at time t₀ and the argument of the perihelion of the object's orbit sin(θ₀+ω) = z₀ / (r₀ sin i) cos(θ₀+ω) = ( x₀ cos Ω + y₀ sin Ω ) / r₀ (θ₀+ω)' = arctan[ sin(θ₀+ω) / cos(θ₀+ω) ] If cos(θ₀+ω) < 0 then θ₀ = θ₀' + π If cos(θ₀+ω) > 0 and sin(θ₀+ω) < 0 then (θ₀+ω) = (θ₀+ω)' + 2π The argument of the perihelion of the object's orbit ω' = (θ₀+ω) − θ₀ if ω'<0 then ω=ω'+2π else ω=ω' The eccentric anomaly of the object at time t₀ cos u₀ = 1 − r₀/(aU) sin u₀ = (x₀ Vx₀ + y₀ Vy₀ + z₀ Vz₀) / √(aμU) u₀' = arctan( sin u₀ / cos u₀ ) If cos u₀ < 0 then u₀ = u₀' + π If cos u₀ > 0 and sin u₀ < 0 then u₀ = u₀' + 2π The mean anomaly of the object at time t₀ M₀ = u₀ − e sin u₀ The period of the object's orbit in days P = 365.256898326 a¹·⁵ The object's time of perihelion passage T = t₀ − PM₀/(2π) Example problem: Find the orbit of Ceres Observation #1 t₁ = JD 2457204.625 X⊕₁ = +0.155228396 AU Y⊕₁ = −1.004732775 AU Z⊕₁ = +0.00003295786 AU α₁ = 20h 46m 57.02s δ₁ = −27°41′33.9″ Observation #2 t₂ = JD 2457214.625 X⊕₂ = +0.319493277 AU Y⊕₂ = −0.965116604 AU Z⊕₂ = +0.0000311269 AU α₂ = 20h 39m 57.10s δ₂ = −28°47′21.5″ Observation #3 t₃ = JD 2457224.625 X⊕₃ = +0.4747795623 AU Y⊕₃ = −0.8983801739 AU Z⊕₃ = +0.00002841127 AU α₃ = 20h 31m 22.81s δ₃ = −29°49′22.7″ Observation #4 t₄ = JD 2457234.625 X⊕₄ = +0.616702829 AU Y⊕₄ = −0.8063620175 AU Z⊕₄ = +0.00002486325 AU α₄ = 20h 22m 06.57s δ₄ = −30°41′57.3″ Converting the observation angles to radians: α₁ = 5.44084723 δ₁ = −0.483329666 α₂ = 5.41030979 δ₂ = −0.502468171 α₃ = 5.37290956 δ₃ = −0.520509058 α₄ = 5.33245865 δ₄ = −0.53580299 Continuing through the procedure, t = 2457219.63 T = 0.001553628 ε = 0.409057547 radians X₁ = −0.155228396 Y₁ = +0.921851498 Z₁ = +0.399597005 X₂ = −0.319493277 Y₂ = +0.885503087 Z₂ = +0.383841559 X₃ = −0.474779562 Y₃ = +0.824271594 Z₃ = +0.357299982 X₄ = −0.616702829 Y₄ = +0.739843884 Z₄ = +0.320703493 a₁ = +0.589463371 b₁ = −0.660726081 c₁ = −0.464730008 a₂ = +0.563195268 b₂ = −0.671477524 c₂ = −0.4815901 a₃ = +0.532276132 b₃ = −0.685093499 c₃ = −0.497321844 a₄ = +0.49965704 b₄ = −0.69978587 c₄ = −0.510531662 R₁² = 1.03358381 R₄² = 1.03054208 2R₁ cos θ₁ = 1.772595 2R₄ cos θ₄ = 1.97920299 τ₁ = 0.344041979 τ₂ = 0.17202099 τ₃ = 0.516062969 τ₄ = 0.17202099 τ₅ = 0.344041979 Φ = −0.058607619 φ = −0.030167527 A = +0.404275125 B = −7.08013785 C = +4.84883362 D = −0.043927505 A' = +1.72864027 B' = −12.7399767 C' = +3.76138574 D' = +0.951283093 E = +2 F = +0.236729767 G = +0.80855025 H = −0.095703956 I = +0.191407912 K = +0.342297675 L = −2.91537363 M = −2.20429551 E' = +0.5 F' = +0.118364883 G' = +0.864320133 H' = +0.102305152 I' = +0.204610303 K' = +0.223373365 L' = −1.86702289 M' = −0.617533885 r₁ = 2.75 (initial guess) r₄ = 2.75 (initial guess) 1st approximation ρ₁ = 1.97723208 ρ₄ = 1.9223289 r₁ = 2.90652064 r₄ = 2.92071388 2nd approximation ρ₁ = 2.001149 ρ₄ = 1.94460337 r₁ = 2.93008666 r₄ = 2.94292188 3rd approximation ρ₁ = 2.00417695 ρ₄ = 1.94741724 r₁ = 2.93307059 r₄ = 2.94572742 4th approximation ρ₁ = 2.00455348 ρ₄ = 1.94776713 r₁ = 2.93344165 r₄ = 2.94607627 5th approximation ρ₁ = 2.0046002 ρ₄ = 1.94781054 r₁ = 2.93348769 r₄ = 2.94611956 6th approximation ρ₁ = 2.00460599 ρ₄ = 1.94781593 r₁ = 2.9334934 r₄ = 2.94612492 7th approximation ρ₁ = 2.00460671 ρ₄ = 1.9478166 r₁ = 2.93349411 r₄ = 2.94612559 8th approximation ρ₁ = 2.0046068 ρ₄ = 1.94781668 r₁ = 2.9334942 r₄ = 2.94612567 9th approximation ρ₁ = 2.00460681 ρ₄ = 1.94781669 r₁ = 2.93349421 r₄ = 2.94612568 10th approximation (final, converged) ρ₁ = 2.00460681 ρ₄ = 1.94781669 r₁ = 2.93349421 r₄ = 2.94612568 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 The reciprocal of the speed of light ç = 0.00577551833 days/AU Aberration corrections to time çρ₁ = 0.011577643 days t₁° = 2457204.61 çρ₄ = 0.011249651 days t₄° = 2457234.61 Epoch of state vector t₀ = 2457219.61 JD The object's state vector in heliocentric ecliptic coordinates 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 The object's distance from the sun at t₀ r₀ = 2.93980995 AU Sunrelative speed V₀ = 16819.9661 m/s The sun's gravitational parameter μ = 1.32712440018ᴇ20 m³ sec⁻² The semimajor axis of the object's orbit a = 2.76694735 AU The angular momentum per unit mass in the object's orbit hx = +1.3390648ᴇ15 m²/sec hy = −2.2844842ᴇ14 m²/sec hz = +7.26435028ᴇ15 m²/sec h = 7.39026848ᴇ15 m²/sec The eccentricity of the object's orbit e = 0.076026341 The inclination of the object's orbit i = 10.5918141° The longitude of the ascending node of the object's orbit Ω = 80.3183813° The true anomaly at time t₀ θ₀ = 147.669798° The sum of the true anomaly at time t₀ and the argument of the perihelion of the object's orbit (θ₀+ω) = 220.296384° The argument of the perihelion of the object's orbit ω = 72.6265867° The eccentric anomaly of the object at time t₀ u₀ = 145.259666° The mean anomaly of the object at time t₀ M₀ = 142.777370° The period of the object's orbit P = 1681.12408 days Times of perihelion passage T₀ = 2456552.87 T₁ = 2458234.01 Summary of the Results of Calculation and Comparison with JPL's Numbers Orbital elements (as calculated) a = 2.76694735 AU e = 0.076026341 i = 10.5918141° Ω = 80.3183813° ω = 72.6265868° T₀ = JD 2456552.87 T₁ = JD 2458234.01 Orbital elements (from the JPL SmallBody Database) a = 2.76916515 AU e = 0.076009027 i = 10.5940672 Ω = 80.3055309° ω = 73.5976947° T = JD 2458238.75 Last edited by Jerry Abbott; September 13th, 2020 at 05:10 PM. 
September 14th, 2020  #2 
Senior Member
Join Date: May 2014
Posts: 9,967

The orbital mechanics and mathematics behind it are beyond me personally, but this is proof that only White people are intelligent enough to figure this out and achieve something like space travel.
__________________
"Military men are dumb, stupid animals to be used as pawns for foreign policy." Henry A. Kissinger, jewish politician and advisor 
September 14th, 2020  #3 
Sex: Male

Maybe we need a math and physics subforum.
__________________
All these ideas…are chained to the existence of men, to who[m]…they owe their existence. Precisely in this case the preservation of these definite races and men is the precondition for the existence of these ideas. Adolf Hitler 