<p>I'm trying to calculate the state vectors of the ISS using Gauss's method with three sets of values (right ascension, declination and sidereal time). I'm following the method described in &quot;Orbital Mechanics for engineering students&quot;. But the values i'm getting from my code are way off from the state vectors given by NASA. I used the on-date RA and Dec and mean sidereal time on 10 June 2023 from stellarium. Any errors I'm missing?</p>
<pre><code>1: RA: 7h 16m 55.04s, Dec: -1* 30’ 51.2’’, Sidereal time: 12h 03m 3.2s 
2: RA: 7h 45m 46.35s, Dec: -23* 11’ 35.5’’, Sidereal time: 12h 03m 56.7s 
3: RA: 8h 17m 38.14s, Dec: -43* 45’ 43.4’’, Sidereal time: 12h 04m 42.6s 
lat/long: 1.347580, 103.843500
</code></pre>
<pre><code># this is the changed one
import numpy as np

todeg = 180 / np.pi
torad = np.pi / 180

phi = 1.347580 * torad
f = 0.003353
h = 30 / 1000
re = 6378
mu = 389600

theta1 = (12 + 3/60 + 3.2/3600) * 15 * torad
theta2 = (12 + 3/60 + 56.7/3600) * 15 * torad
theta3 = (12 + 4/60 + 42.6/3600) * 15 * torad

ra1 = (7 + 16/60 + 55.04/3600) * 15 * torad
dec1 = (-1 - 30/60 - 51.2/3600) * torad
ra2 = (7 + 45/60 + 46.35/3600) * 15 * torad
dec2 = (-23 - 11/60 - 35.5/3600) * torad
ra3 = (8 + 17/60 + 38.14/3600) * 15 * torad
dec3 = (-43 - 45/60 - 43.4/3600) * torad

t1 = 0
t2 = 53.35
t3 = 99.4

def posroot(roots):
    # Extract positive real roots
    posroots = roots[np.logical_and(roots &gt; 0, np.isreal(roots))].real
    
    # Check if any positive real roots exist
    npositive = len(posroots)
    if npositive == 0:
        raise ValueError(&quot;There are no positive real roots.&quot;)
    
    # If there is only one positive real root, return it
    if npositive == 1:
        return posroots[0]
    
    # If there are multiple positive real roots, prompt the user to select one
    print(&quot;There are multiple positive real roots. Please select one:&quot;)
    for i, root in enumerate(posroots):
        print(f&quot;Root #{i+1}: {root}&quot;)
    
    while True:
        choice = input(&quot;Enter the root number to use: &quot;)
        try:
            choice = int(choice)
            if 1 &lt;= choice &lt;= npositive:
                return posroots[choice-1]
            else:
                print(&quot;Invalid choice. Please try again.&quot;)
        except ValueError:
            print(&quot;Invalid input. Please enter a valid root number.&quot;)

def position_vector(theta):
    a = 1 - (2 * f - f**2) * (np.sin(phi))**2
    b = np.sqrt(a)
    c = re / b
    D = c + h
    e = c * (1 - f)**2
    F = e + h

    Rx = D * np.cos(phi) * np.cos(theta)
    Ry = D * np.cos(phi) * np.sin(theta)
    Rz = F * np.sin(phi)
    R = np.array([Rx, Ry, Rz])

    return R

R1 = position_vector(theta1)
R2 = position_vector(theta2)
R3 = position_vector(theta3)

def rhohat(ra, dec):
    rx = np.cos(dec) * np.cos(ra)
    ry = np.cos(dec) * np.sin(ra)
    rz = np.sin(dec)
    rhohat = np.array([rx, ry, rz])
    return rhohat

rhohat1 = rhohat(ra1, dec1)
rhohat2 = rhohat(ra2, dec2)
rhohat3 = rhohat(ra3, dec3)

tau1 = t1 - t2
tau3 = t3 - t2
tau = tau3 - tau1

p1 = np.cross(rhohat2, rhohat3)
p2 = np.cross(rhohat1, rhohat3)
p3 = np.cross(rhohat1, rhohat2)

D0 = np.dot(rhohat1, p1)
D11 = np.dot(R1, p1)
D21 = np.dot(R2, p1)
D31 = np.dot(R3, p1)
D12 = np.dot(R1, p2)
D22 = np.dot(R2, p2)
D32 = np.dot(R3, p2)
D13 = np.dot(R1, p3)
D23 = np.dot(R2, p3)
D33 = np.dot(R3, p3)

A = 1 / D0 * (-1 * D12 * tau3/tau + D22 + D32 * tau1/tau)
B = 1 / 6 / D0 * (D12 * (tau3**2-tau*2) * tau3/tau + D32 * (tau**2 - tau1**2) * tau1 / tau)

E = np.dot(R2, rhohat2)
R2squared = np.dot(R2, R2)

a = -1 * (A**2 + 2 * A * E + R2squared)
b = -2 * mu * B * (A + E)
c = -1 * mu**2 * B**2

r2 = posroot(np.roots([1, 0, a, 0, 0, b, 0, 0, c]))

rho1 = 1 / D0 * ((6 * (D31 * tau1 / tau3 + D21 * tau / tau3) * r2**3 + mu * D31 * (tau**2 - tau1**2) * tau1 / tau3)/(6 * r2**3 + mu * (tau**2 - tau3**2))- D11)
rho2 = A + mu * B / r2**3
rho3 = 1 / D0 * ((6 * (D13 * tau3 / tau1 - D23 * tau / tau1) * r2**3 + mu * D13 * (tau**2 - tau3**2) * tau3 / tau1)/(6 * r2**3 + mu * (tau**2 - tau1**1))- D33)

r1vec = R1 + rho1 * rhohat1
r2vec = R2 + rho2 * rhohat2
r3vec = R3 + rho3 * rhohat3
print(r1vec)
print(r2vec)
print(r3vec)
</code></pre>
<p>Errors - thanks for your help :)</p>
<ol>
<li>tau3 should be t3 - t2</li>
<li>declination values should not be multiplied by 15</li>
</ol>

CC BY-SA 4.0
2023-06-20T14:17:46.853Z
I don't know this algorithm, but I suspect that you shouldn't be multiplying your declinations in degrees by 15. (You do need that for the RA, to convert hours to degrees).
PM 2Ring
CC BY-SA 4.0
2023-06-20T19:15:45.020Z
It would help to put the NASA vectors you're comparing them to.
Greg Miller
CC BY-SA 4.0
2023-06-21T00:12:41.037Z
you may find other references to Gauss' method in [Space Exploration SE](https://space.stackexchange.com/search?q=gauss+method) and here in [Astronomy SE](https://astronomy.stackexchange.com/search?q=gauss+method) helpful.
uhoh
CC BY-SA 4.0
2023-06-22T04:29:11.577Z
@PM2Ring thanks for that, it really helped
specimen_n1l3
<orbital-mechanics><orbital-elements><celestial-mechanics>
Code for Gauss's method orbit determination
specimen_n1l3
# Code for Gauss's method orbit determination

(Asked by: username_0 on 2023-06-20)

I'm trying to calculate the state vectors of the ISS using Gauss's method with three sets of values (right ascension, declination and sidereal time). I'm following the method described in "Orbital Mechanics for engineering students". But the values i'm getting from my code are way off from the state vectors given by NASA. I used the on-date RA and Dec and mean sidereal time on 10 June 2023 from stellarium. Any errors I'm missing?
1: RA: 7h 16m 55.04s, Dec: -1* 30’ 51.2’’, Sidereal time: 12h 03m 3.2s 
2: RA: 7h 45m 46.35s, Dec: -23* 11’ 35.5’’, Sidereal time: 12h 03m 56.7s 
3: RA: 8h 17m 38.14s, Dec: -43* 45’ 43.4’’, Sidereal time: 12h 04m 42.6s 
lat/long: 1.347580, 103.843500

# this is the changed one
import numpy as np

todeg = 180 / np.pi
torad = np.pi / 180

phi = 1.347580 * torad
f = 0.003353
h = 30 / 1000
re = 6378
mu = 389600

theta1 = (12 + 3/60 + 3.2/3600) * 15 * torad
theta2 = (12 + 3/60 + 56.7/3600) * 15 * torad
theta3 = (12 + 4/60 + 42.6/3600) * 15 * torad

ra1 = (7 + 16/60 + 55.04/3600) * 15 * torad
dec1 = (-1 - 30/60 - 51.2/3600) * torad
ra2 = (7 + 45/60 + 46.35/3600) * 15 * torad
dec2 = (-23 - 11/60 - 35.5/3600) * torad
ra3 = (8 + 17/60 + 38.14/3600) * 15 * torad
dec3 = (-43 - 45/60 - 43.4/3600) * torad

t1 = 0
t2 = 53.35
t3 = 99.4

def posroot(roots):
    # Extract positive real roots
    posroots = roots[np.logical_and(roots > 0, np.isreal(roots))].real
    
    # Check if any positive real roots exist
    npositive = len(posroots)
    if npositive == 0:
        raise ValueError("There are no positive real roots.")
    
    # If there is only one positive real root, return it
    if npositive == 1:
        return posroots[0]
    
    # If there are multiple positive real roots, prompt the user to select one
    print("There are multiple positive real roots. Please select one:")
    for i, root in enumerate(posroots):
        print(f"Root #{i+1}: {root}")
    
    while True:
        choice = input("Enter the root number to use: ")
        try:
            choice = int(choice)
            if 1 <= choice <= npositive:
                return posroots[choice-1]
            else:
                print("Invalid choice. Please try again.")
        except ValueError:
            print("Invalid input. Please enter a valid root number.")

def position_vector(theta):
    a = 1 - (2 * f - f**2) * (np.sin(phi))**2
    b = np.sqrt(a)
    c = re / b
    D = c + h
    e = c * (1 - f)**2
    F = e + h

    Rx = D * np.cos(phi) * np.cos(theta)
    Ry = D * np.cos(phi) * np.sin(theta)
    Rz = F * np.sin(phi)
    R = np.array([Rx, Ry, Rz])

    return R

R1 = position_vector(theta1)
R2 = position_vector(theta2)
R3 = position_vector(theta3)

def rhohat(ra, dec):
    rx = np.cos(dec) * np.cos(ra)
    ry = np.cos(dec) * np.sin(ra)
    rz = np.sin(dec)
    rhohat = np.array([rx, ry, rz])
    return rhohat

rhohat1 = rhohat(ra1, dec1)
rhohat2 = rhohat(ra2, dec2)
rhohat3 = rhohat(ra3, dec3)

tau1 = t1 - t2
tau3 = t3 - t2
tau = tau3 - tau1

p1 = np.cross(rhohat2, rhohat3)
p2 = np.cross(rhohat1, rhohat3)
p3 = np.cross(rhohat1, rhohat2)

D0 = np.dot(rhohat1, p1)
D11 = np.dot(R1, p1)
D21 = np.dot(R2, p1)
D31 = np.dot(R3, p1)
D12 = np.dot(R1, p2)
D22 = np.dot(R2, p2)
D32 = np.dot(R3, p2)
D13 = np.dot(R1, p3)
D23 = np.dot(R2, p3)
D33 = np.dot(R3, p3)

A = 1 / D0 * (-1 * D12 * tau3/tau + D22 + D32 * tau1/tau)
B = 1 / 6 / D0 * (D12 * (tau3**2-tau*2) * tau3/tau + D32 * (tau**2 - tau1**2) * tau1 / tau)

E = np.dot(R2, rhohat2)
R2squared = np.dot(R2, R2)

a = -1 * (A**2 + 2 * A * E + R2squared)
b = -2 * mu * B * (A + E)
c = -1 * mu**2 * B**2

r2 = posroot(np.roots([1, 0, a, 0, 0, b, 0, 0, c]))

rho1 = 1 / D0 * ((6 * (D31 * tau1 / tau3 + D21 * tau / tau3) * r2**3 + mu * D31 * (tau**2 - tau1**2) * tau1 / tau3)/(6 * r2**3 + mu * (tau**2 - tau3**2))- D11)
rho2 = A + mu * B / r2**3
rho3 = 1 / D0 * ((6 * (D13 * tau3 / tau1 - D23 * tau / tau1) * r2**3 + mu * D13 * (tau**2 - tau3**2) * tau3 / tau1)/(6 * r2**3 + mu * (tau**2 - tau1**1))- D33)

r1vec = R1 + rho1 * rhohat1
r2vec = R2 + rho2 * rhohat2
r3vec = R3 + rho3 * rhohat3
print(r1vec)
print(r2vec)
print(r3vec)

Errors - thanks for your help :)

tau3 should be t3 - t2
declination values should not be multiplied by 15

**Tags:** orbital-mechanics, orbital-elements, celestial-mechanics

> **Comments:**
> **username_1**: I don't know this algorithm, but I suspect that you shouldn't be multiplying your declinations in degrees by 15. (You do need that for the RA, to convert hours to degrees).  
> **username_2**: It would help to put the NASA vectors you're comparing them to.  
> **username_3**: you may find other references to Gauss' method in [Space Exploration SE](https://space.stackexchange.com/search?q=gauss+method) and here in [Astronomy SE](https://astronomy.stackexchange.com/search?q=gauss+method) helpful.  
> **username_0**: @PM2Ring thanks for that, it really helped  

---

## 0 answers

