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.

When neglecting the influence of other objects, two celestial bodies orbit one another along a conic trajectory. In the orbital plane, the radial equation is thus:

 r = l/(1 + e cos(angle)) 

The   l   (unity),   e   and   angle   are respectively called the   semi-latus rectum,   the eccentricity   and the   true anomaly.   The eccentricity and the true anomaly are two of the six so-called   orbital elements   often used to specify an orbit and the position of a point on this orbit.

The four other parameters are the   semi-major axis,   the   longitude of the ascending node,   the   inclination   and the   argument of periapsis.

The semi-major axis is half the distance between   perihelion and aphelion.   It is often noted   a,   and it's not too hard to see how it's related to the semi-latus-rectum:

 a = l/(1 - e2) 

The longitude of the ascending node, the inclination and the argument of the periapsis specify the orientation of the orbiting plane with respect to a reference plane defined with three arbitrarily chosen reference distant stars.

Those six parameters, along with dynamical considerations implying notably the   vis-viva equation,   allow for the determination of both the position and the speed of the orbiting object in   cartesian coordinates,   those two vectors constituting the so-called   orbital state vectors.

The determination of the speed also requires the knowledge of the so-called   gravitational parameter   which we shall consider equal to one here:

 µ = GM = 1 

The purpose of this task is to show how to perform this conversion from orbital elements to orbital state vectors in your programming language.

TODO:   pick an example from a reputable source, and bring the algorithm description onto this site. (Restating those pages in concise a fashion comprehensible to the coders and readers of this site will be a good exercise.)

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>sub orbital-state-vectors(

   Real :$semimajor-axis where * >= 0,
   Real :$eccentricity   where * >= 0,
   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) = .cos, .sin given $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 /= sqrt($speed**2);
   $speed *= sqrt(2/r - 1/$semimajor-axis);

   { :$position, :$speed }

}

say orbital-state-vectors

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

</lang>

Output:

{position => 0.237771283982207*e0+0.860960261697716*e1+0.110509023572076*e2, speed => -1.06193301748006*e0+0.27585002056925*e1+0.135747024865598*e2}

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