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

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)}

zkl

Translation of: Perl

<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>

Output:

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