Orbital elements: Difference between revisions
(Added D) |
|||
| Line 124: | Line 124: | ||
writeln("Position : ", res[0]); |
writeln("Position : ", res[0]); |
||
writeln("Speed : ", res[1]); |
writeln("Speed : ", res[1]); |
||
writeln("Edit source/app.d to start your project."); |
|||
}</lang> |
}</lang> |
||
{{out}} |
{{out}} |
||
Revision as of 01:35, 24 January 2018
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. An other parameter, called the gravitational parameter, along with dynamical considerations described further, also allows for the determination of the speed of the orbiting object.
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.
The gravitational parameter is the coefficent GM in Newton's gravitational force. It is sometimes noted µ and will be chosen as one here for the sake of simplicity:
µ = GM = 1
As mentioned, dynamical considerations allow for the determination of the speed. They result in the so-called vis-viva equation:
v2 = GM(2/r - 1/a)
This only gives the magnitude of the speed. The direction is easily determined since it's tangent to the conic.
Those parameters 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.
- Task
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.)
D
<lang D>import std.math; import std.stdio; import std.typecons;
struct Vector {
double x, y, z;
auto opBinary(string op : "+")(Vector rhs) {
return Vector(x+rhs.x, y+rhs.y, z+rhs.z);
}
auto opBinary(string op : "*")(double m) {
return Vector(x*m, y*m, z*m);
}
auto opOpAssign(string op : "*")(double m) {
this.x *= m;
this.y *= m;
this.z *= m;
return this;
}
auto opBinary(string op : "/")(double d) {
return Vector(x/d, y/d, z/d);
}
auto opOpAssign(string op : "/")(double m) {
this.x /= m;
this.y /= m;
this.z /= m;
return this;
}
auto abs() {
return sqrt(x * x + y * y + z * z);
}
void toString(scope void delegate(const(char)[]) sink) const {
import std.format;
sink("(");
formattedWrite(sink, "%.16f", x);
sink(", ");
formattedWrite(sink, "%.16f", y);
sink(", ");
formattedWrite(sink, "%.16f", z);
sink(")");
}
}
auto orbitalStateVectors(
double semiMajorAxis, double eccentricity, double inclination, double longitudeOfAscendingNode, double argumentOfPeriapsis, double trueAnomaly
) {
auto i = Vector(1.0, 0.0, 0.0); auto j = Vector(0.0, 1.0, 0.0); auto k = Vector(0.0, 0.0, 1.0);
auto mulAdd = (Vector v1, double x1, Vector v2, double x2) => v1 * x1 + v2 * x2;
auto rotate = (Vector i, Vector j, double alpha) =>
tuple(mulAdd(i, +cos(alpha), j, sin(alpha)),
mulAdd(i, -sin(alpha), j, cos(alpha)));
auto p = rotate(i, j, longitudeOfAscendingNode); i = p[0]; j = p[1]; p = rotate(j, k, inclination); j = p[0]; p = rotate(i, j, argumentOfPeriapsis); i = p[0]; j = p[1];
auto l = semiMajorAxis * ((eccentricity == 1.0) ? 2.0 : (1.0 - eccentricity * eccentricity)); auto c = cos(trueAnomaly); auto s = sin(trueAnomaly); auto r = l / (1.0 + eccentricity * c); auto rprime = s * r * r / l; auto position = mulAdd(i, c, j, s) * r; auto speed = mulAdd(i, rprime * c - r * s, j, rprime * s + r * c); speed /= speed.abs(); speed *= sqrt(2.0 / r - 1.0 / semiMajorAxis); return tuple(position, speed);
}
void main() {
auto res = orbitalStateVectors(1.0, 0.1, 0.0, 355.0 / (113.0 * 6.0), 0.0, 0.0);
writeln("Position : ", res[0]);
writeln("Speed : ", res[1]);
}</lang>
- Output:
Position : (0.7794228433986798, 0.4500000346536842, 0.0000000000000000) Speed : (-0.5527708409604437, 0.9574270831797614, 0.0000000000000000) Edit source/app.d to start your project.
Kotlin
<lang scala>// version 1.1.4-3
class Vector(val x: Double, val y: Double, val z: Double) {
operator fun plus(other: Vector) = Vector(x + other.x, y + other.y, z + other.z) operator fun times(m: Double) = Vector(x * m, y * m, z * m)
operator fun div(d: Double) = this * (1.0 / d)
fun abs() = Math.sqrt(x * x + y * y + z * z)
override fun toString() = "($x, $y, $z)"
}
fun orbitalStateVectors(
semimajorAxis: Double, eccentricity: Double, inclination: Double, longitudeOfAscendingNode: Double, argumentOfPeriapsis: Double, trueAnomaly: Double
): Pair<Vector, Vector> {
var i = Vector(1.0, 0.0, 0.0) var j = Vector(0.0, 1.0, 0.0) var k = Vector(0.0, 0.0, 1.0)
fun mulAdd(v1: Vector, x1: Double, v2: Vector, x2: Double) = v1 * x1 + v2 * x2
fun rotate(i: Vector, j: Vector, alpha: Double) =
Pair(mulAdd(i, +Math.cos(alpha), j, Math.sin(alpha)),
mulAdd(i, -Math.sin(alpha), j, Math.cos(alpha)))
var p = rotate(i, j, longitudeOfAscendingNode) i = p.first; j = p.second p = rotate(j, k, inclination) j = p.first p = rotate(i, j, argumentOfPeriapsis) i = p.first; j = p.second
val l = semimajorAxis * (if (eccentricity == 1.0) 2.0 else (1.0 - eccentricity * eccentricity)) val c = Math.cos(trueAnomaly) val s = Math.sin(trueAnomaly) val r = l / (1.0 + eccentricity * c) val rprime = s * r * r / l val position = mulAdd(i, c, j, s) * r var speed = mulAdd(i, rprime * c - r * s, j, rprime * s + r * c) speed /= speed.abs() speed *= Math.sqrt(2.0 / r - 1.0 / semimajorAxis) return Pair(position, speed)
}
fun main(args: Array<String>) {
val (position, speed) = orbitalStateVectors(
semimajorAxis = 1.0,
eccentricity = 0.1,
inclination = 0.0,
longitudeOfAscendingNode = 355.0 / (113.0 * 6.0),
argumentOfPeriapsis = 0.0,
trueAnomaly = 0.0
)
println("Position : $position")
println("Speed : $speed")
}</lang>
- Output:
Position : (0.7794228433986797, 0.45000003465368416, 0.0) Speed : (-0.5527708409604438, 0.9574270831797618, 0.0)
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}
Sidef
<lang ruby>func orbital_state_vectors(
semimajor_axis, eccentricity, inclination, longitude_of_ascending_node, argument_of_periapsis, true_anomaly
) {
static vec = frequire('Math::Vector::Real')
var (i, j, k) = (vec.V(1,0,0), vec.V(0,1,0), vec.V(0,0,1))
func muladd(v1, x1, v2, x2) {
v1.mul(x1).add(v2.mul(x2))
}
func rotate(Ref i, Ref j, α) {
(*i, *j) = (
muladd(*i, +cos(α), *j, sin(α)),
muladd(*i, -sin(α), *j, cos(α)),
)
}
rotate(\i, \j, longitude_of_ascending_node) rotate(\j, \k, inclination) rotate(\i, \j, argument_of_periapsis)
var l = (eccentricity == 1 ? 2*semimajor_axis
: semimajor_axis*(1 - eccentricity**2))
var (c, s) = with(true_anomaly) { (.cos, .sin) }
var r = l/(1 + eccentricity*c) var rprime = (s * r**2 / l) var position = muladd(i, c, j, s).mul(r)
var speed = muladd(i, rprime*c - r*s, j, rprime*s + r*c) speed.div!(abs(speed)) speed.mul!(sqrt(2/r - 1/semimajor_axis))
struct Result { position, speed }
Result([position,speed].map {|v| [v{:module}[]].map{Num(_)} }...)
}
var r = orbital_state_vectors(
semimajor_axis: 1, eccentricity: 0.1, inclination: 0, longitude_of_ascending_node: 355/(113*6), argument_of_periapsis: 0, true_anomaly: 0,
)
say '['+r.position.join(', ')+']' say '['+r.speed.join(', ')+']'</lang>
- Output:
[0.77942284339868, 0.450000034653684, 0] [-0.552770840960444, 0.957427083179761, 0]
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))