Orbital elements: Difference between revisions
m (added whitespace and highlighting to the task's preamble, used a shaded window for the formulae.) |
(various rephrasing + changing 'l' to 'L' in order to avoid confusion with 1) |
||
Line 3: | Line 3: | ||
When neglecting the influence of other objects, two celestial bodies orbit one another along a [[wp:conic section|conic]] trajectory. In the orbital plane, the radial equation is thus: |
When neglecting the influence of other objects, two celestial bodies orbit one another along a [[wp:conic section|conic]] trajectory. In the orbital plane, the radial equation is thus: |
||
<big> r = |
<big> r = L/(1 + e cos(angle)) </big> |
||
<big> '''<tt>L</tt>''' </big>, <big> '''e''' </big> and <big> '''angle''' </big> are respectively called ''semi-latus rectum'', ''eccentricity'' and ''true anomaly''. The eccentricity and the true anomaly are two of the six so-called [[wp:orbital elements|orbital elements]] often used to specify an orbit and the position of a point on this orbit. |
|||
The four other parameters are the |
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 |
The semi-major axis is half the distance between [[wp:perihelion and aphelion|perihelion and aphelion]]. It is often noted <big> '''a'''</big>, and it's not too hard to see how it's related to the semi-latus-rectum: |
||
<big> a = |
<big> a = L/(1 - e<sup>2</sup>) </big> |
||
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. |
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 |
Those six parameters, along with dynamical considerations explained below, allow for the determination of both the position and the speed of the orbiting object in [[wp:cartesian coordinates|cartesian coordinates]], those two vectors constituting the so-called [[wp:orbital state vectors|orbital state vectors]]. |
||
The aforementioned dynamical considerations imply the so-called [[wp:vis-viva equation|vis-viva equation]]: |
|||
The determination of the speed also requires the knowledge of the so-called ''gravitational parameter'' which we shall consider equal to one here: |
|||
<big>v<sup>2</sup> = GM(2/r - 1/a)</big> |
|||
<tt>GM</tt> is the gravitational parameter sometimes noted <tt>µ</tt>. It will be chosen as one here for the sake of simplicity: |
|||
<big> µ = GM = 1 </big> |
<big> µ = GM = 1 </big> |
||
Line 23: | Line 27: | ||
The purpose of this task is to show how to perform this conversion from orbital elements to orbital state vectors in your programming language. |
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: |
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.) |
||
<br><br> |
<br><br> |
||
Revision as of 10:17, 4 August 2016
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))
L , e and angle are respectively called semi-latus rectum, eccentricity and 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 explained below, 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 aforementioned dynamical considerations imply the so-called vis-viva equation:
v2 = GM(2/r - 1/a)
GM is the gravitational parameter sometimes noted µ. It will be chosen as one here for the sake of simplicity:
µ = 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
<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
<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))