Orbital elements

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

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.

TODO: pick an example from a reputable source.

Perl

Translation of: Perl 6

<lang perl>use strict; use warnings; use Math::Vector::Real;

sub orbital_state_vectors {

   my (
       $semimajor_axis,
       $eccentricity,
       $inclination,
       $longitude_of_ascending_node,
       $argument_of_periapsis,
       $true_anomaly
   ) = @_[0..5];

   my ($i, $j, $k) = (V(1,0,0), V(0,1,0), V(0,0,1));
   
   sub rotate {
       my $alpha = shift;
       @_[0,1] = (
           +cos($alpha)*$_[0] + sin($alpha)*$_[1],
           -sin($alpha)*$_[0] + cos($alpha)*$_[1]
       );
   }

   rotate $longitude_of_ascending_node, $i, $j;
   rotate $inclination,                 $j, $k;
   rotate $argument_of_periapsis,       $i, $j;

   my $l = $eccentricity == 1 ? # PARABOLIC CASE
       2*$semimajor_axis :
       $semimajor_axis*(1 - $eccentricity**2);

   my ($c, $s) = (cos($true_anomaly), sin($true_anomaly));

   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 /= abs($speed);
   $speed *= sqrt(2/$r - 1/$semimajor_axis);

   {
       position => $position,
       speed    => $speed
   }

}

use Data::Dumper;

print Dumper orbital_state_vectors

   1,                             # semimajor axis
   0.1,                           # eccentricity
   0,                             # inclination
   355/113/6,                     # longitude of ascending node
   0,                             # argument of periapsis
   0                              # true-anomaly
   ;

</lang>

Output:

$VAR1 = {
          'position' => bless( [
                                 '0.77942284339868',
                                 '0.450000034653684',
                                 '0'
                               ], 'Math::Vector::Real' ),
          'speed' => bless( [
                              '-0.552770840960444',
                              '0.957427083179762',
                              '0'
                            ], 'Math::Vector::Real' )
        };

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**2);
   $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)}