The heliocentric position of Voyager 1 at 12h UT on 13 March 2009 (JD 2454904) was

r = 109.17 AU
hlong = 254.146 degrees
hlat = +34.858 degrees

Here's how to do the calculation.

Voyager 1, orbital elements.
a = -3.22185 AU
e = 3.72466
i = 35.777 degrees
Ω = 179.001 degrees
ω = 338.272 degrees
T = JD 2444227.38

Hyperbolic mean anomaly.
M = 0.01720209895 (t - T) sqrt[1/(-a)³]

Hyperbolic eccentric anomaly. (Danby's method shown.)

u = 0
Repeat...
. u₀ = u
. f₀ = e sinh u₀ - u₀ - M
. f₁ = e cosh u₀ - 1
. f₂ = e sinh u₀
. f3 = e cosh u₀
. d₁ = -f₀ / f₁
. d₂ = -f₀ / [ f₁ + (d₁)(f₂)/2 ]
. d₃ = -f₀ / [ f₁ + (d₁)(f₂)/2 + (d₂)² (f₃)/6 ]
. u = u₀ + d₃
Until |u - u₀| 0 and y>0 then δ=0
if x>0 and y1.e-14
:u→u0
:e\*sinh(u0)-u0-m→f0
:e\*cosh(u0)-1→f1
:e\*sinh(u0)→f2
:e\*cosh(u0)→f3
:-f0/f1→d1
:-f0/(f1+d1\*f2)→d2
:-f0/(f1+d1\*f2/2+d2^2\*f3/6)→d3
:u0+d3→u
:EndWhile
:cos⁻¹((e-cosh(u))/(e\*cosh(u)-1))→q
:a\*(1-e\*cosh(u))→r
:r\*cos(q)→xp
:r\*sin(q)→yp
:xp\*cos(w)-yp\*sin(w)→x
:xp\*sin(w)+yp\*cos(w)→y
:x→xp
:y\*cos(i)→yp
:y\*sin(i)→zp
:xp\*cos(L)-yp\*sin(L)→x
:xp\*sin(L)+yp\*cos(L)→y
:zp→z
:(180/π)\*tan⁻¹(y/x)→lh
:If x0 and y<0 Then
:lh+360→lh
:EndIf
:(180/π)\*sin⁻¹(z/r)→bh
:Output 1,1,nam
:Output 11,20,x
:Output 21.20,y
:Output 31,20,z
:Output 41,1,r
:Output 51,1,lh
:Output 61,1,bh
:EndPrgm


THE DATA FILE.

TI-89 Titanium data file "probes" - used in program probe(jd)

r1c1=-5.663407
r1c2=1.39939
r1c3=2.346188
r1c4=229.543416
r1c5=287.930748
r1c6=2453772.55
r1c7=newhoriz
r2c1=-3.22185
r2c2=3.72466
r2c3=35.777
r2c4=179.001
r2c5=338.272
r2c6=2444227.38
r2c7=voyager1
r3c1=-4.02568
r3c2=6.27722
r3c3=78.7783
r3c4=101.6368
r3c5=130.06246
r3c6=2445454.94
r3c7=voyager2


EXECUTION COMMAND (from the TI-89T HOME screen).

probe(2,2454904)


OUTPUT (on the IO screen).

voyager1
-24.473811
-86.175767
62.397398
109.172641
254.145519
34.858239


The first three numbers are the rectangular position in heliocentric ecliptic coordinates, x,y,z. The last three numbers are the position in heliocentric spherical coordinates, r,hlong,hlat.

Jerry Abbott

By accident, some of the text I wrote got deleted from the post right before I sent it. Here is a revision.


The heliocentric position of Voyager 1 at 12h UT on 13 March 2009 (JD 2454904) was

r = 109.17 AU
hlong = 254.146 degrees
hlat = +34.858 degrees

Here's how to do the calculation.

Voyager 1, orbital elements.
a = -3.22185 AU
e = 3.72466
i = 35.777 degrees
Ω = 179.001 degrees
ω = 338.272 degrees
T = JD 2444227.38

Hyperbolic mean anomaly.
M = 0.01720209895 (t - T) sqrt[1/(-a)³]

Hyperbolic eccentric anomaly. (Danby's method shown.)

u = 0
Repeat...
. u₀ = u
. f₀ = e sinh u₀ - u₀ - M
. f₁ = e cosh u₀ - 1
. f₂ = e sinh u₀
. f3 = e cosh u₀
. d₁ = -f₀ / f₁
. d₂ = -f₀ / [ f₁ + (d₁)(f₂)/2 ]
. d₃ = -f₀ / [ f₁ + (d₁)(f₂)/2 + (d₂)² (f₃)/6 ]
. u = u₀ + d₃
Until |u - u₀| 0 and y>0 then δ=0
if x>0 and y1.e-14

True anomaly.
Q = arccos { (e - cosh u) / (e cosh u - 1) }

Heliocentric distance.
r = a (1 - e cosh u)

Canonical position vector.
x''' = r cos Q
y''' = r sin Q

Rotate by argument of the perihelion.
x'' = x''' cos ω - y''' sin ω
y'' = x''' sin ω + y''' cos ω

Rotate by inclination.
x' = x''
y' = y'' cos i
z' = y'' sin i

Rotate by longitude of ascending node to get...
The heliocentric position in ecliptic coordinates.
x = x' cos Ω - y' sin Ω
y = x' sin Ω + y' cos Ω
z = z'

Heliocentric longitude.
hlong = (180/π) { arctan( y / x ) + δ } degrees
if x>0 and y>0 then δ=0
if x>0 and y<0 then δ=2π
if x<0 then δ=π

Heliocentric latitude.
hlat = (180/π) { arcsin( z / r ) } degrees

I won't repeat the TI-89T calculator program that I included in the original post because of the way the posting software treats asterisks.

Jerry Abbott