Orbital elements

Orbital elements is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The purpose of this task is to convert orbital elements into orbital state vectors. It will be assumed that µ = GM = 1.

TODO: pick an example from a reputable source.

Perl 6

We'll use the Clifford geometric algebra library but only for the vector operations.

<lang perl6>subset PositiveReal of Real where * >= 0; sub orbital-state-vectors(

   PositiveReal :$semimajor-axis,
   PositiveReal :$eccentricity,
   Real :$inclination,
   Real :$longitude-of-ascending-node,
   Real :$argument-of-periapsis,
   Real :$true-anomaly

) {

   use Clifford;
   my ($i, $j, $k) = @e[^3];

   sub rotate($a is rw, $b is rw, Real \α) {
       ($a, $b) = cos(α)*$a + sin(α)*$b, -sin(α)*$a + cos(α)*$b;
   }
   rotate($i, $j, $longitude-of-ascending-node);
   rotate($j, $k, $inclination);
   rotate($i, $j, $argument-of-periapsis);

   my \l = $eccentricity == 1 ?? # PARABOLIC CASE
   2*$semimajor-axis !!
   $semimajor-axis*(1 - $eccentricity**2);

   my ($c, $s) = map {.($true-anomaly)}, &cos, &sin;

   my \r = l/(1 + $eccentricity*$c);
   my \rprime = $s*r**2/l;

   my $position = r*($c*$i + $s*$j);

   my $speed = 
   (rprime*$c - r*$s)*$i + (rprime*$s + r*$c)*$j;
   $speed /= sqrt(($speed·$speed).Real);
   $speed *= sqrt(2/r - 1/$semimajor-axis);

   { :position($position X· @e[^3]), :speed($speed X· @e[^3]) }

}

say orbital-state-vectors

   semimajor-axis => 1,
   eccentricity => 0.1,
   inclination => 0,
   longitude-of-ascending-node => pi/6,
   argument-of-periapsis => 0,
   true-anomaly => 0;

</lang>

Output:

{position => (0.779422863405995 0.45 0), speed => (-0.552770798392567 0.957427107756338 0)}