Classical orbital elements Vectors

If I would be given a R and a V vector and I would have to find orbital elements like... (a=semi-major axis)(eccentricity)(inclination)(right ascension of the ascending node)(argument of perigee)(true anomaly). After I defined the vectors and found the magnitudes I was trying to write out my equations but I guess matlab didn't like it. If you could look at my code and point me in the right direction I would appreciate it. I had quite a few of errors just to keep ending my functions over and over again?

First, solve for the angular momentum: h⃗ =r⃗ ×v⃗ h

The eccentricity vector is then: e⃗ =(v2−μ/r)r⃗ −(r⃗ ⋅v⃗ )v⃗

a=−μ/2E

i=cos−1hKhi = cos−1⁡hKh

Ω=cos−1nInΩ = cos−1⁡nIn

ω=cos−1n⃗ ⋅e⃗ neω = cos−1⁡n→⋅e→ne

ν=cos−1e⃗ ⋅r⃗ er

R= [-7953.8073 - 4174.5370 - 1008.9496];

v= [3.6460035 - 4.9118820 - 4.9193608];

h=cross(R,v);

nhat=cross([0 0 1],h)

r=norm(R);

mu=3.986*10^5

energy = mag(v)^2/2-mu/mag(R)

e = mag(evec)

evec = ((mag(v)^2-mu/mag(R))*R-dot(R,v)*v)/mu

if abs(e-1.0)>eps

a = -mu/(2*energy)

p = a*(1-e^2)

else

p = mag(h)^2/mu

a = inf

end

i = acos(h(3)/mag(h))

Omega = acos(n(1)/mag(n))

if n(2)<0

Omega = 360-Omega

argp = acos(dot(n,evec)/(mag(n)*e))

end

if e(3)<0

argp = 360-argp

nu = acos(dot(evec,R)/(e*mag(R))

end

if dot(R,v)<0

nu = 360 - nu

end

Answer 1

[James Tursa]

Use the norm( ) function instead of mag( ). E.g., norm(v) instead of mag(v).

Calculate evec before you calculate e.

Your evec code appears correct, but your description of this is missing the mu term in the denominator. E.g.,

e⃗ =((v2−μ/r)r⃗ −(r⃗ ⋅v⃗ )v⃗ )/μ

[Oskar Kinat]

keep getting this error.

MATLAB_code.m Line: 36 Column: 34

Invalid expression. When calling a function or indexing a variable, use

parentheses. Otherwise, check for mismatched delimiters. I don't understand.

R= [-7953.8073 -4174.5370 -1008.9496];

V= [3.6460035 -4.9118820 -4.9193608];

h=cross(R,V);

nhat=cross([0 0 1],h);

r=norm(R);

v=norm(V);

mu=3.986*10^5;

evec = ((norm(V)^2-mu/norm(R))*R-dot(R,V)*V)/mu;

e = mag(evec);

energy = mag(V)^2/2-mu/mag(R);

if abs(e-1.0)>eps

a = -mu/(2*energy);

p = a*(1-e^2);

else

p = norm(h)^2/mu;

a = inf;

end

i = acos(h(3)/mag(h));

Omega = acos(n(1)/norm(n));

if n(2)<0

Omega = 360-Omega;

argp = acos(dot(n,evec)/(norm(n)*e));

end

if e(3)<0

argp = 360-argp;

nu = acos(dot(evec,R)/(e*norm(R)); % line 36 error

end

if dot(R,V)<0

nu = 360 - nu;

end

[James Tursa]

The error message is telling you the exact problem ... you have mismatched parentheses. So add a closing paren:

nu = acos(dot(evec,R)/(e*norm(R)));