Orbital elements: Difference between revisions

<lang zkl>fcn orbital_state_vectors(semimajor_axis, eccentricity, inclination,

longitude_of_ascending_node, argument_of_periapsis, true_anomaly){

i,j,k:=T(1.0, 0.0, 0.0), T(0.0, 1.0, 0.0), T(0.0, 0.0, 1.0);

vdot:=fcn(c,vector){ vector.apply('*,c) };

vsum:=fcn(v1,v2) { v1.zipWith('+,v2) };

rotate:='wrap(alpha, a,b){ // a&b are vectors: (x,y,z)

return(vsum(vdot( alpha.cos(),a), vdot(alpha.sin(),b)), #cos(alpha)*a + sin(alpha)*b

vsum(vdot(-alpha.sin(),a), vdot(alpha.cos(),b)));

};

i,j=rotate(longitude_of_ascending_node,i,j);

j,k=rotate(inclination, j,k);

i,j=rotate(argument_of_periapsis, i,j);

l:=if(eccentricity==1) # PARABOLIC CASE

semimajor_axis*2 else

semimajor_axis*(1.0 - eccentricity.pow(2));;

c,s,r:=true_anomaly.cos(), true_anomaly.sin(), l/(eccentricity*c + 1);

rprime:=s*r.pow(2)/l;

position:=vdot(r,vsum(vdot(c,i), vdot(s,j))); #r*(c*i + s*j)

speed:=vsum(vdot(rprime*c - r*s,i), vdot(rprime*s + r*c,j)); #(rprime*c - r*s)*i + (rprime*s + r*c)*j

z:=speed.zipWith('*,speed).sum(0.0).sqrt(); #sqrt(speed**2)

speed=vdot(1.0/z,speed); #speed/z

speed=vdot((2.0/r - 1.0/semimajor_axis).sqrt(),speed); #speed*sqrt(2/r - 1/semimajor_axis)

return(position,speed);

}</lang>

<lang zkl>orbital_state_vectors(

1.0, # semimajor axis

0.1, # eccentricity

0.0, # inclination

(0.0).pi/6, # longitude of ascending node

0.0, # argument of periapsis

0.0 # true-anomaly

).println();</lang>

{{out}}

<pre>

L(L(0.779423,0.45,0),L(-0.552771,0.957427,0))

</pre>