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

Old March 27th, 2017 #1
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default A funny thing happened on the Physics Forums

I was censored on the Physics Forums for a reason so persnickety that it seems suspicious. Read the contents of the thread that I posted (to follow), only to have it removed just as I was finishing up.
 
Old March 27th, 2017 #2
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

The observations in this exposition of the method of Gauss (for finding preliminary Keplerian elements of a planet's orbit) are fictional. The data were taken from JPL's Horizons interface program at

https://ssd.jpl.nasa.gov/horizons.cgi

Quote:
Current Settings
Ephemeris Type [change] : OBSERVER
Target Body [change] : Mars Barycenter [4]
Observer Location [change] : Geocentric [500]
Time Span [change] : Start=2018-07-20, Stop=2018-08-03, Step=7 d
Table Settings [change]: defaults
Display/Output [change] : default (formatted HTML)
Quote:
Current Settings
Ephemeris Type [change] : VECTORS
Target Body [change] : Earth-Moon Barycenter [EMB] [3]
Coordinate Origin [change] : Solar System Barycenter (SSB) [500@0]
Time Span [change] : Start=2018-07-20, Stop=2018-08-03, Step=7 d
Table Settings [change] : defaults
Display/Output [change] : default (formatted HTML)
 
Old March 27th, 2017 #3
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

First Observation

At 0h UT on 20 July 2018, an observer on Earth sees Mars at 20h 39m 21.03s right ascension and −24° 50' 00.0" declination.

The time, t₁, converted to Julian Date:
t₁ = JD 2458319.5

The obliquity of the ecliptic at time t₁
ε₁ = 0.40905071279426947 radians

Earth's heliocentric position in ecliptic coordinates at time t₁

x⊕₁ = +0.4618841708199115 AU
y⊕₁ = −0.8984126365037669 AU
z⊕₁ = −0.00005165456609240359 AU

The sun's position in geocentric celestial coordinates at time t₁ is calculated by a rotation and a translation of the coordinate system.

X๏₁ = −0.4618841708199115 AU
Y๏₁ = +0.8242719747171781 AU
Z๏₁ = +0.3573807210716430 AU

A unit vector in the apparent direction of Mars in geocentric celestial coordinates at time t₁

a₁ = +0.5813786066519955
b₁ = −0.6968612501535729
c₁ = −0.4199801349609095
 
Old March 27th, 2017 #4
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

Second Observation

At 0h UT on 27 July 2018, Mars can be seen at 20h 31m 47.37s right ascension and −25° 32' 35.7" declination.

The time, t₂, converted to Julian Date:
t₂ = JD 2458326.5

The obliquity of the ecliptic at time t₂
ε₂ = 0.4090506693018622 radians

Earth's heliocentric position in ecliptic coordinates at time t₂

x⊕₂ = +0.5637478350603771 AU
y⊕₂ = −0.8380340129569841 AU
z⊕₂ = −0.00005336797521227511 AU

The sun's position in geocentric celestial coordinates at time t₂, therefore, is

X๏₂ = −0.5637478350603771 AU
Y๏₂ = +0.7688739928199521 AU
Z๏₂ = +0.33336735425957914 AU

A unit vector in the apparent direction of Mars in geocentric celestial coordinates at time t₂

a₂ = +0.5548334913212998
b₂ = −0.7115005281989506
c₂ = −0.43119229501561324
 
Old March 27th, 2017 #5
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

Third Observation

At 0h UT on 03 August 2018, Mars can be seen at 20h 23m 48.06s right ascension and −26° 06' 15.8" declination.

The time, t₃, converted to Julian Date:
t₃ = JD 2458333.5

The obliquity of the ecliptic at time t₃
ε₃ = 0.4090506258094553 radians

Earth's heliocentric position in ecliptic coordinates at time t₃

x⊕₃ = +0.6578217701347128 AU
y⊕₃ = −0.7659716360160346 AU
z⊕₃ = −0.00005554489987022230 AU

The sun's position in geocentric celestial coordinates at time t₃, therefore, is

X๏₃ = −0.6578217701347128 AU
Y๏₃ = +0.702755994303274 AU
Z๏₃ = +0.3047073394868152 AU

A unit vector in the apparent direction of Mars in geocentric celestial coordinates at time t₃

a₃ = +0.5271965424026115
b₃ = −0.7269503439654071
c₃ = −0.44000795798179343
 
Old March 27th, 2017 #6
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

Useful constants and conversions

GM = 1.32712440018e+20 m³ sec⁻²
k = 0.01720209895
A = 1/c = 0.00577551833109 days/AU
1 AU = 1.49597870691e+11 meters
1 AU/day = 1731456.8368056 m/s
 
Old March 27th, 2017 #7
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

First Approximation for the heliocentric and geocentric distances (equations)

τ₁ = k(t₃−t₂)
τ₂ = k(t₃−t₁)
τ₃ = k(t₂−t₁)

n₁° = τ₁/τ₂
n₃° = τ₃/τ₂

ν₁ = τ₁τ₃ (1+n₁°) / 6
ν₃ = τ₁τ₃ (1+n₃°) / 6

D = a₂ (b₃c₁−b₁c₃) + b₂ (c₃a₁−c₁a₃) + c₂ (a₃b₁−a₁b₃)

For i = 1 to 3
dᵢ = X๏ᵢ (c₁b₃−c₃b₁) + Y๏ᵢ (a₁c₃−a₃c₁) + Z๏ᵢ (b₁a₃−b₃a₁)
R๏ᵢ = √(X๏ᵢ² + Y๏ᵢ² + Z๏ᵢ²)
qᵢ = −2(aᵢX๏ᵢ + bᵢY๏ᵢ + cᵢZ๏ᵢ)

K₀ = (d₂ − d₁n₁° − d₃n₃°)/D
L₀ = (d₁ν₁ + d₃ν₃)/D

r₂ = the distance from the sun to the planet at time 2.
ρ₂ = the distance from Earth to the planet at time 2.

We must solve these two equations simultaneously for r₂ and ρ₂.

ρ₂ = K₀ − L₀/r₂³
r₂ = √(R๏₂² + q₂ρ₂ + ρ₂²)

This leads to an 8th degree polynomial in r₂.

F{r₂} = r₂⁸ − (R๏₂² + q₂K₀ + K₀²) r₂⁶ + L₀ (q₂ + 2K₀)r₂³ − L₀²
F{r₂} = 0

We will use Newton's method to get the answer.

dF/dr₂ = 8 r₂⁷ − 6 (R๏₂² + q₂K₀ + K₀²) r₂⁵ + 3 L₀ (q₂ + 2K₀) r₂²

As an initial guess at r₂,

r₂(j=0) = 5.0

Then...

Repeat over j
r₂(j+1) = r₂(j) − F{r₂(j)} / dF/dr₂(j)
Until | r₂(j+1) − r₂(j) | < 1.0e−15

Assign to r₂ the converged value of r₂ from the loop.

r₂ = r₂(j+1)

ρ₂ = K₀ − L₀/r₂³

n₁ = n₁° + ν₁/r₂³
n₃ = n₃° + ν₃/r₂³

Q₁ = n₁X๏₁ − X๏₂ + n₃X๏₃ + a₂ρ₂
Q₂ = n₁Y๏₁ − Y๏₂ + n₃Y๏₃ + b₂ρ₂
Q₃ = n₁Z๏₁ − Z๏₂ + n₃Z๏₃ + c₂ρ₂
Q₄ = a₁ − a₃c₁/c₃
Q₅ = b₁ − b₃a₁/a₃
Q₆ = b₃ − b₁a₃/a₁
Q₇ = a₃ − a₁c₃/c₁

ρ₁ = ½ { (Q₁−a₃Q₃/c₃) / (n₁Q₄) + (Q₂−b₃Q₁/a₃) / (n₁Q₅) }
ρ₃ = ½ { (Q₂−b₁Q₁/a₁) / (n₃Q₆) + (Q₁−a₁Q₃/c₁) / (n₃Q₇) }

r₁ = √(R๏₁² + q₁ρ₁ + ρ₁²)
r₃ = √(R๏₃² + q₃ρ₃ + ρ₃²)
 
Old March 27th, 2017 #8
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

First Approximation for the heliocentric and geocentric distances (example)

τ₁ = 0.12041469265
τ₂ = 0.24082938530
τ₃ = 0.12041469265

n₁° = 0.5
n₃° = 0.5

ν₁ = 0.003624924551498492
ν₃ = 0.003624924551498492

D = −0.00008086527122733167

d₁ = −0.007999178173536615
d₂ = −0.007285847860243794
d₃ = −0.006470768826225836

R๏₁ = 1.0101892175982878
R๏₂ = 1.0100062530777665
R๏₃ = 1.0096742204936617

q₁ = 1.9860511558039105
q₂ = 2.0071717322553524
q₃ = 1.983487457488723

K₀ = 0.6291249579754575
L₀ = 0.648640205396262

F{r₂} = r₂⁸ − 2.6786726757084853 r₂⁶ + 2.1180837685979137 r₂³ − 0.42073411605650496

r₂(j=0) = 5.0
r₂(j=1) = 4.3929125172136540
r₂(j=2) = 3.8645627337430860
r₂(j=3) = 3.4056413590103425
r₂(j=4) = 3.0081300591803277
r₂(j=5) = 2.6651781699814600
r₂(j=6) = 2.3710083783224194
r₂(j=7) = 2.1208545400998022
r₂(j=8) = 1.9109403021959235
r₂(j=9) = 1.7385176091531582
r₂(j=10) = 1.6019958403129406
r₂(j=11) = 1.5011460045069969
r₂(j=12) = 1.4368768296650747
r₂(j=13) = 1.4079301074280257
r₂(j=14) = 1.4021251234310480
r₂(j=15) = 1.4019087711519191
r₂(j=16) = 1.4019084781470565
r₂(j=17) = 1.4019084781465196
r₂(j=18) = 1.4019084781465194
r₂(j=19) = 1.4019084781465194

r₂ = 1.40190847814652 AU

ρ₂ = 0.39370413259633 AU

n₁ = 0.5013156488339009
n₃ = 0.5013156488339009

Q₁ = +0.22086196320358292
Q₂ = −0.28347167403935020
Q₃ = −0.17121443708688030
Q₄ = +0.07817844626231640
Q₅ = +0.10480064268380263
Q₆ = −0.09503365936122998
Q₇ = −0.08190658470375745

ρ₁ = 0.40112680805014 AU
ρ₃ = 0.39332223978467 AU

r₁ = 1.40642928447905 AU
r₃ = 1.39796070946325 AU
 
Old March 27th, 2017 #9
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

Correction for Planetary Aberration (equations)

Light doesn't travel instantaneously. We must correct the initial times for the amount of time it took light to travel from the planet to our eyes. This is done only once, after the first approximation, but before the 2nd approximation. It will not be done between any other successive approximations.

For i = 1 to 3
tᵢ° = tᵢ − A ρᵢ

Where A is the time (in days) required for light to travel 1 astronomical unit.

τ₁ = k(t₃°−t₂°)
τ₂ = k(t₃°−t₁°)
τ₃ = k(t₂°−t₁°)

Where k is the mean motion of Earth in radians per day.

n₁° = τ₁/τ₂
n₃° = τ₃/τ₂
 
Old March 27th, 2017 #10
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

Correction for Planetary Aberration (example)

t₁° = JD 2458319.4976832850
t₂° = JD 2458326.4977261545
t₃° = JD 2458333.4977283604

τ₁ = 0.12041473059503525
τ₂ = 0.24083016069401889
τ₃ = 0.12041543009898362

n₁° = 0.4999985477235360
n₃° = 0.5000014522764639
 
Old March 27th, 2017 #11
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

Recursive Procedure for Successive Approximations (equations)

For i = 1 to 3
xᵢ = aᵢρᵢ − X๏ᵢ
yᵢ = bᵢρᵢ − Y๏ᵢ
zᵢ = cᵢρᵢ − Z๏ᵢ

K₁ = √{ 2 (r₂r₃ + x₂x₃ + y₂y₃ + z₂z₃) }
K₂ = √{ 2 (r₁r₃ + x₁x₃ + y₁y₃ + z₁z₃) }
K₃ = √{ 2 (r₁r₂ + x₁x₂ + y₁y₂ + z₁z₂) }

h₁ = τ₁² / { K₁² [K₁/3 + (r₂+r₃)/2] }
h₂ = τ₂² / { K₂² [K₂/3 + (r₁+r₃)/2] }
h₃ = τ₃² / { K₃² [K₃/3 + (r₁+r₂)/2] }

For i = 1 to 3
Begin
ζ₁ = (11/9) hᵢ
ζ₂ = ζ₁
Repeat
ζ₃ = ζ₂
ζ₂ = ζ₁ / ( 1 + ζ₂ )
Until |ζ₂ − ζ₃| / ζ₃ < 1e−15
Ψᵢ = 1 + (10/11) ζ₂
End

n₁ = n₁° Ψ₂/Ψ₁
n₃ = n₃° Ψ₂/Ψ₃

ν₁ = n₁° r₂³ (Ψ₂/Ψ₁−1)
ν₃ = n₃° r₂³ (Ψ₂/Ψ₃−1)

K₀ = (d₂ − d₁n₁° − d₃n₃°) / D
L₀ = (d₁ν₁ + d₃ν₃) / D

ρ₂ = K₀ − L₀/r₂³
r₂ = √(R๏₂² + q₂ρ₂ + ρ₂²)

F{r₂} = r₂⁸ − (R๏₂² + q₂K₀ + K₀²) r₂⁶ + L₀ (q₂ + 2K₀) r₂³ − L₀²
F{r₂} = 0

dF/dr₂ = 8r₂⁷ − 6 (R๏₂² + q₂K₀ + K₀²) r₂⁵ + 3 L₀ (q₂ + 2K₀) r₂²

An initial guess at r₂ is made.

r₂(j=0) = 5.0

Repeat over j
r₂(j+1) = r₂(j) − F{r₂(j)} / dF/dr₂(j)
Until | r₂(j+1) − r₂(j) | < 1.0e−15

Assign to r₂ the converged value from the loop.

r₂ = r₂(j+1)

ρ₂ = K₀ − L₀/r₂³

Q₁ = n₁X๏₁ − X๏₂ + n₃X๏₃ + a₂p₂
Q₂ = n₁Y๏₁ − Y๏₂ + n₃Y๏₃ + b₂p₂
Q₃ = n₁Z๏₁ − Z๏₂ + n₃Z๏₃ + c₂p₂
Q₄ = a₁ − a₃c₁/c₃
Q₅ = b₁ − b₃a₁/a₃
Q₆ = b₃ − b₁a₃/a₁
Q₇ = a₃ − a₁c₃/c₁

ρ₁ = (1/2) { (Q₁−a₃Q₃/c₃) / (n₁Q₄) + (Q₂−b₃Q₁/a₃) / (n₁Q₅) }
ρ₃ = (1/2) { (Q₂−b₁Q₁/a₁) / (n₃Q₆) + (Q₁−a₁Q₃/c₁) / (n₃Q₇) }

r₁ = √(R๏₁² + q₁ρ₁ + ρ₁²)
r₃ = √(R๏₃² + q₃ρ₃ + ρ₃²)

Repeat this whole section (Recursive Procedure for Successive Approximations) until your values for r converge.
 
Old March 27th, 2017 #12
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

2nd Approximation

x₁ = +0.6950907155748627 AU
y₁ = −1.1038017036451089 AU
z₁ = −0.5258460120529784 AU

x₂ = +0.7821880734964239 AU
y₂ = −1.0489946911163521 AU
z₂ = −0.5031295427509228 AU

x₃ = +0.8651798949992420 AU
y₃ = −0.9886817318039844 AU
z₃ = −0.4777722550432934 AU

K₁ = 2.7978740465198397
K₂ = 2.7964458183802980
K₃ = 2.8063597760501007

h₁ = 0.0007940909165695694
h₂ = 0.0031771970467915564
h₃ = 0.0007869223244419950

Ψ₁ = 1.0008814685551377
Ψ₂ = 1.0035166156963282
Ψ₃ = 1.0008735187996662

n₁ = 0.5013149570937037
n₃ = 0.5013218511700457

ν₁ = 0.0036270200092790930
ν₃ = 0.0036380120924249812

K₀ = 0.6291524070013451
L₀ = 0.6498947413724790

[Run the loop to converge r₂ with these improved values for K₀ and L₀.]

r₂ = 1.40104915810870 AU

ρ₂ = 0.39284197048693 AU

Q₁ = +0.22037984626244440
Q₂ = −0.28285445669627170
Q₃ = −0.17084103674551712
Q₄ = +0.07817844626231640
Q₅ = +0.10480064268380263
Q₆ = −0.09503365936122998
Q₇ = −0.08190658470375745

ρ₁ = 0.40023158610550 AU
ρ₃ = 0.39248103668010 AU

r₁ = 1.40554188199514 AU
r₃ = 1.39712727023708 AU
 
Old March 27th, 2017 #13
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

3rd Approximation

x₁ = +0.6945702526880453 AU
y₁ = −1.1031778581616043 AU
z₁ = −0.5254700366198498 AU

x₂ = +0.7817097170831797 AU
y₂ = −1.0483812623201196 AU
z₂ = −0.5027577850922944 AU

x₃ = +0.8647364156310531 AU
y₃ = −0.9880702189177708 AU
z₃ = −0.4774021189830026 AU

K₁ = 2.7961791000789040
K₂ = 2.7947160071365147
K₃ = 2.8046107513968490

h₁ = 0.0007955352634619711
h₂ = 0.0031830908801306914
h₃ = 0.0007883948794899752

Ψ₁ = 1.0008830702760896
Ψ₂ = 1.0035231140784986
Ψ₃ = 1.0008751518307448

n₁ = 0.5013174011504072
n₃ = 0.5013242795711852

ν₁ = 0.0036270759743293342
ν₃ = 0.0036380048010743680

K₀ = 0.6291524070013451
L₀ = 0.6498996939776956

r₂ = 1.40104553690590 AU

ρ₂ = 0.39283833730013 AU

[Skipping the Q values to save space.]

ρ₁ = 0.40021552948029 AU
ρ₃ = 0.39247674784955 AU

r₁ = 1.40552596570965 AU
r₃ = 1.39712302101410 AU
 
Old March 27th, 2017 #14
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

4th Approximation

[Hereafter showing only the resulting values for r and ρ]

r₁ = 1.40551810179968 AU
r₂ = 1.40103775126222 AU
r₃ = 1.39711525882522 AU

ρ₁ = 0.40020759623016 AU
ρ₂ = 0.39283052589046 AU
ρ₃ = 0.39246891330763 AU

5th Approximation

r₁ = 1.40551793902222 AU
r₂ = 1.40103770073617 AU
r₃ = 1.39711520326079 AU

ρ₁ = 0.40020743201740 AU
ρ₂ = 0.39283047519720 AU
ρ₃ = 0.39246885722527 AU

6th Approximation

r₁ = 1.40551786734425 AU
r₂ = 1.40103763002107 AU
r₃ = 1.39711513274520 AU

ρ₁ = 0.40020735970740 AU
ρ₂ = 0.39283040424807 AU
ρ₃ = 0.39246878605239 AU

7th Approximation

r₁ = 1.40551786570544 AU
r₂ = 1.40103762940236 AU
r₃ = 1.39711513208127 AU

ρ₁ = 0.40020735805414 AU
ρ₂ = 0.39283040362731 AU
ρ₃ = 0.39246878538226 AU

8th Approximation

r₁ = 1.40551786505173 AU
r₂ = 1.40103762875969 AU
r₃ = 1.39711513144030 AU

ρ₁ = 0.40020735739466 AU
ρ₂ = 0.39283040298252 AU
ρ₃ = 0.39246878473532 AU

9th Approximation

r₁ = 1.40551786503536 AU
r₂ = 1.40103762875261 AU
r₃ = 1.39711513143281 AU

ρ₁ = 0.40020735737816 AU
ρ₂ = 0.39283040297541 AU
ρ₃ = 0.39246878472776 AU

10th Approximation

r₁ = 1.40551786502943 AU
r₂ = 1.40103762874680 AU
r₃ = 1.39711513142702 AU

ρ₁ = 0.40020735737217 AU
ρ₂ = 0.39283040296958 AU
ρ₃ = 0.39246878472191 AU

11th Approximation

r₁ = 1.40551786502925 AU
r₂ = 1.40103762874670 AU
r₃ = 1.39711513142691 AU

ρ₁ = 0.40020735737199 AU
ρ₂ = 0.39283040296948 AU
ρ₃ = 0.39246878472181 AU

12th Approximation

r₁ = 1.40551786502921 AU
r₂ = 1.40103762874666 AU
r₃ = 1.39711513142687 AU

ρ₁ = 0.40020735737195 AU
ρ₂ = 0.39283040296944 AU
ρ₃ = 0.39246878472177 AU

13th Approximation

r₁ = 1.40551786502920 AU
r₂ = 1.40103762874665 AU
r₃ = 1.39711513142687 AU

ρ₁ = 0.40020735737195 AU
ρ₂ = 0.39283040296944 AU
ρ₃ = 0.39246878472177 AU

The geocentric and heliocentric distances of Mars have converged.
 
Old March 27th, 2017 #15
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

The positions of the object in heliocentric celestial coordinates, and the velocity at time t₂

For i=1 to 3
xᵢ = aᵢρᵢ − X๏ᵢ
yᵢ = bᵢρᵢ − Y๏ᵢ
zᵢ = cᵢρᵢ − Z๏ᵢ

x₁ = +0.6945561666206916 AU
y₁ = −1.1031609740960513 AU
z₁ = −0.5254598610330625 AU

x₂ = +0.7817032990370643 AU
y₂ = −1.0483730320253148 AU
z₂ = −0.5027527972678799 AU

x₃ = +0.8647299564409844 AU
y₃ = −0.9880613123524494 AU
z₃ = −0.47739672802383704 AU

The obliquity of the ecliptic at the times corrected for planetary aberration

ţᵢ = (tᵢ°−2451545)/3652500

ţ₁ = 0.001854756381460590
ţ₂ = 0.001856672888748669
ţ₃ = 0.001858589384903587

εᵢ" = 84381.448"
− 4680.93" ţ
− 1.55" ţ²
+ 1999.25" ţ³
− 51.38" ţ⁴
− 249.67" ţ⁵
− 39.05" ţ⁶
+ 7.12" ţ⁷
+ 27.87" ţ⁸
+ 5.79" ţ⁹
+ 2.45" ţ¹⁰

[Source: J. Laskar via Wikipedia.]

εᵢ = 4.84813681109536e-6 ε"

ε₁ = 0.4090507128082722 radians
ε₂ = 0.4090506693155986 radians
ε₃ = 0.4090506258231779 radians

The positions of the object in heliocentric ecliptic coordinates

Xᵢ = xᵢ
Yᵢ = zᵢ sin qᵢ + yᵢ cos qᵢ
Zᵢ = zᵢ cos qᵢ − yᵢ sin qᵢ

X₁ = +0.6945561666206916 AU
Y₁ = −1.2211445112611437 AU
Z₁ = −0.043339161761816236 AU

X₂ = +0.7817032990370643 AU
Y₂ = −1.1618451636648004 AU
Z₂ = −0.04429678439054979 AU

X₃ = +0.8647299564409844 AU
Y₃ = −1.0964241450956385 AU
Z₃ = −0.04502096119374227 AU

Estimating the velocity at time 2 using Lagrange interpolation

w₁ = (t₂°−t₃°) / { (t₁°−t₂°) (t₁°−t₃°) }
w₂ = (2t₂°−t₁°−t₃°) / { (t₂°−t₁°) (t₂°−t₃°) }
w₃ = (t₂°−t₁°) / { (t₃°−t₁°) (t₃°−t₂°) }

w₁ = −0.07142792651844013
w₂ = −0.0000008298695931895751
w₃ = +0.07142875638803331

Here's a conversion factor that converts a velocity from AU/day to meters per second.

λ = 1731456.8368056

Vx₂ = λ {w₁X₁ + w₂X₂ + w₃X₃}
Vy₂ = λ {w₁Y₁ + w₂Y₂ + w₃Y₃}
Vz₂ = λ {w₁Z₁ + w₂Z₂ + w₃Z₃}

Vx₂ = +21046.2558357560 m/s
Vy₂ = +15424.8069263065 m/s
Vz₂ = −207.9965286294 m/s
 
Old March 27th, 2017 #16
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

From the state vector to the orbital elements (equations)

We remove the subscripts and wingdings from the position and velocity of the object at time t₂.

The distance of the object from the sun
r = √( x² + y² + z² )

The speed of the object, relative to the sun
V = √[ (Vx)² + (Vy)² + (Vz)² ]

The semimajor axis of the object's orbit around the sun

a = 1 / { 2/r − V²/(GM) }

The angular momentum, h, per unit mass (m²/sec), in the 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 around the sun

e = √[ 1 − h² / (GMa) ]

The inclination of the object's orbit, relative to the ecliptic

i = arccos( hz / h )

The longitude of the ascending node, Ω, of the object's orbit

Ω = Arctan ( hx , −hy )

The true anomaly, θ, of the object in the orbit at time t₂

cos θ = h²/(rGM) − 1
sin θ = h (x Vx + y Vy + z Vz) / (rGM)
θ = Arctan ( sin θ , cos θ )

(Arctan is the two-dimensional arctangent function.)

Argument of the perihelion, ω, of the orbit

cos ω'' = (x cos Ω + y sin Ω) / r
If sin i = 0 then sin ω'' = (y cos Ω − x sin Ω) / r
If sin i ≠ 0 then sin ω'' = z / (r sin i)
ω'' = Arctan( sin ω'' , cos ω'' )
ω' = ω'' − θ
If ω' > 0 then ω = ω'
If ω' < 0 then ω = ω' + 2π

The eccentric anomaly, u, of the object in the orbit at time t₂

sin u = (r/a) sin θ / √(1−e²)
cos u = (r/a) cos θ + e
u = Arctan( sin u , cos u )

The mean anomaly, m, of the object in the orbit at time t₂

m = u − e sin u

(Note: u must be entered in radians, and m will return in radians.)

The period of the orbit, P, days

P = (π / 43200) √[ a³/(GM) ]

The mean motion, μ, in the orbit, radians per day

μ = 2π / P

The time of perihelion passage, T, of the object in the orbit, Julian Date

T = t − m/μ
 
Old March 27th, 2017 #17
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

From the state vector to the orbital elements (example)

t = JD 2458326.4977261545
x = +1.16941149055e+11 meters
y = −1.73809562567e+11 meters
z = −6.626704624e+9 meters
Vx = +21046.2558357560 m/s
Vy = +15424.8069263065 m/s
Vz = −207.9965286294 m/s

r = 2.09592246031e+11 meters = 1.4010376287467 AU
V = 26094.3061983615 m/s

The estimated semimajor axis of Mars' orbit

a = 1.51523 AU

The estimated angular momentum per unit mass for Mars in its orbit

hx = +1.3836742502626161e+14 m²/sec
hy = −1.1514396779369656e+14 m²/sec
hz = +5.4618351660801300e+15 m²/sec

The estimated eccentricity of Mars' orbit

e = 0.085238

The estimated inclination of Mars' orbit to the ecliptic

i = 1.88766°

The estimated longitude of the ascending node of Mars' orbit

Ω = 50.23409°

The estimated true anomaly of Mars in its orbit (at time t₂)

θ = 329.77123°

The estimated argument of the perihelion for Mars' orbit

ω = 283.93619°

The estimated eccentric anomaly of Mars in its orbit (at time t₂)

u = 332.14649°

The estimated mean anomaly of Mars in its orbit (at time t₂)

m = 334.42826°

The estimated time of perihelion passage for Mars

T = JD 2458374.890
 
Old March 27th, 2017 #18
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

The official NASA values of the Keplerian elements for Mars' orbit are

a = 1.52365 AU
e = 0.093327
i = 1.84815°
Ω = 49.50637°
ω = 286.65831°
T = JD 2458378.005

Using the Keplerian elements that I've calculated (shown in the previous post), an observer on Earth would find Mars, at 0h UT on Christmas 2018, at

Right ascension: 23h 40m 54s
Declination −2° 38' 2"

The position for this same object at the same time, from JPL's Horizons interface is

Right ascension: 23h 42m 19.67s
Declination: −2° 25' 10.2"

Last edited by Jerry Abbott; March 27th, 2017 at 10:28 AM.
 
Old March 27th, 2017 #19
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

The notation used in this demonstration of the method of Gauss closely follows that used by A.D. Dubyago in chapter five of his book The Determination of Orbits, until the post headed "From the state vector to the orbital elements (equations)." From that point forward, the procedure is the one I learned as an undergraduate.
 
Old March 27th, 2017 #20
Jerry Abbott
Senior Member
 
Join Date: Nov 2007
Location: In the hills north of Hillsboro WV
Posts: 1,048
Default

Physics Forums moderator fresh_42 wrote:

Thread closed for moderation.
 
Reply

Share


Thread
Display Modes


All times are GMT -5. The time now is 06:59 AM.
Page generated in 0.90697 seconds.