Rosetta Code

Orbital elements: Difference between revisions

Page Discussion

Orbital elements (view source)

Revision as of 04:07, 24 March 2024

11,530 bytes added ,  1 month ago
Added FreeBASIC
m (→‎{{header|J}}: put that comma where it was originally supposed to be)
(Added FreeBASIC)
 
(8 intermediate revisions by 6 users not shown)
Line 1:
{{draft task}}
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[[w:Polar coordinate system|polar equation]] is thus:
 
&nbsp;<big> r = L/(1 + e cos(angle)) </big>
Line 35:
{{trans|Python}}
 
<langsyntaxhighlight lang="11l">F mulAdd(v1, x1, v2, x2)
R v1 * x1 + v2 * x2
 
Line 69:
V ps = orbitalStateVectors(1.0, 0.1, 0.0, 355.0 / (113.0 * 6.0), 0.0, 0.0)
print(‘Position : ’ps[0])
print(‘Speed : ’ps[1])</langsyntaxhighlight>
 
{{out}}
Line 79:
=={{header|Ada}}==
{{Trans|Kotlin}}
<langsyntaxhighlight Adalang="ada">with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Generic_Real_Arrays;
with Ada.Numerics.Generic_Elementary_Functions;
Line 184:
Put ("Position : "); Put (State.Position); New_Line;
Put ("Speed : "); Put (State.Speed); New_Line;
end Orbit;</langsyntaxhighlight>
{{out}}
<pre>
Position : ( 0.779423, 0.450000, 0.000000)
Speed : (-0.552771, 0.957427, 0.000000)
</pre>
 
=={{header|ALGOL 68}}==
{{Trans|Lua}} (which is a translation of C which is...)
<syntaxhighlight lang="algol68">
BEGIN # orbital elements #
 
MODE VECTOR = STRUCT( REAL x, y, z );
 
OP + = ( VECTOR v, w )VECTOR: ( x OF v + x OF w, y OF v + y OF w, z OF v + z OF w );
OP * = ( VECTOR v, REAL m )VECTOR: ( x OF v * m, y OF v * m, z OF v * m );
OP / = ( VECTOR v, REAL d )VECTOR: v * ( 1 / d );
OP ABS = ( VECTOR v )REAL: sqrt( x OF v * x OF v + y OF v * y OF v + z OF v * z OF v );
 
PROC muladd = ( VECTOR v1, v2, REAL x1, x2 )VECTOR: ( v1 * x1 ) + ( v2 * x2 );
PROC set = ( REF VECTOR v, w, []VECTOR ps )VOID: BEGIN v := ps[ LWB ps ]; w := ps[ UPB ps ] END;
PROC rotate = ( VECTOR i, j, REAL alpha )[]VECTOR:
( muladd( i, j, cos( alpha ), sin( alpha ) ), muladd( i, j, -sin( alpha ), cos( alpha ) ) );
PROC orbital state vectors = ( REAL semimajor axis, eccentricity, inclination
, longitude of ascending node, argument of periapsis
, true anomaly
, REF VECTOR position, speed
) VOID:
BEGIN
VECTOR i := ( 1.0, 0.0, 0.0 ), j := ( 0.0, 1.0, 0.0 ), k := ( 0.0, 0.0, 1.0 );
set( i, j, rotate( i, j, longitude of ascending node ) );
set( j, LOC VECTOR, rotate( j, k, inclination ) );
set( i, j, rotate( i, j, argument of periapsis ) );
REAL l = IF eccentricity /= 1 THEN 1 - eccentricity * eccentricity ELSE 2 FI
* semimajor axis;
REAL c = cos( true anomaly ), s = sin( true anomaly );
REAL r = l / ( 1.0 + eccentricity * c );
REAL rprime = s * r * r / l;
position := muladd( i, j, c, s ) * r;
speed := muladd( i, j, rprime * c - r * s, rprime * s + r * c );
speed := speed / ABS speed;
speed := speed * sqrt( 2 / r - 1 / semimajor axis )
END # orbital state vectors # ;
 
OP FMT = ( REAL v )STRING:
BEGIN
STRING result := fixed( ABS v, 0, 15 );
IF result[ LWB result ] = "." THEN "0" +=: result FI;
WHILE result[ UPB result ] = "0" DO result := result[ : UPB result - 1 ] OD;
IF result[ UPB result ] = "." THEN result := result[ : UPB result - 1 ] FI;
IF v < 0 THEN "-" + result ELSE result FI
END # FMT # ;
OP TOSTRING = ( VECTOR v )STRING: "(" + FMT x OF v + ", " + FMT y OF v + ", " + FMT z OF v + ")";
 
REAL longitude = 355 / ( 113 * 6 );
VECTOR position, speed;
orbital state vectors( 1.0, 0.1, 0.0, longitude, 0.0, 0.0, position, speed );
print( ( "Position : ", TOSTRING position, newline ) );
print( ( "Speed : ", TOSTRING speed ) )
END
</syntaxhighlight>
{{out}}
<pre>
Position : (0.77942284339868, 0.450000034653684, 0)
Speed : (-0.552770840960444, 0.957427083179762, 0)
</pre>
 
=={{header|ALGOL W}}==
{{Trans|C}} (which is a translation of Kotlin which is a translation of ...).
<langsyntaxhighlight lang="algolw">begin
% compute orbital elements %
% 3-element vector %
Line 277 ⟶ 337:
write( "Speed : " ); writeOnVector( speed )
end
end.</langsyntaxhighlight>
{{out}}
<pre>
Line 285 ⟶ 345:
=={{header|C}}==
{{trans|Kotlin}}
<langsyntaxhighlight lang="c">#include <stdio.h>
#include <math.h>
 
Line 360 ⟶ 420:
printf("Speed : %s\n", buffer);
return 0;
}</langsyntaxhighlight>
 
{{output}}
Line 370 ⟶ 430:
=={{header|C sharp|C#}}==
{{trans|D}}
<langsyntaxhighlight lang="csharp">using System;
 
namespace OrbitalElements {
Line 455 ⟶ 515:
}
}
}</langsyntaxhighlight>
{{out}}
<pre>Position : (0.77942284339868, 0.450000034653684, 0)
Line 462 ⟶ 522:
=={{header|C++}}==
{{trans|C#}}
<langsyntaxhighlight lang="cpp">#include <iostream>
#include <tuple>
 
Line 558 ⟶ 618:
 
return 0;
}</langsyntaxhighlight>
{{out}}
<pre>Position : (0.779423, 0.45, 0)
Line 565 ⟶ 625:
=={{header|D}}==
{{trans|Kotlin}}
<langsyntaxhighlight Dlang="d">import std.math;
import std.stdio;
import std.typecons;
Line 653 ⟶ 713:
writeln("Position : ", res[0]);
writeln("Speed : ", res[1]);
}</langsyntaxhighlight>
{{out}}
<pre>Position : (0.7794228433986798, 0.4500000346536842, 0.0000000000000000)
Speed : (-0.5527708409604437, 0.9574270831797614, 0.0000000000000000)</pre>
 
=={{header|FreeBASIC}}==
{{trans|Phix}}
<syntaxhighlight lang="vbnet">Sub vabs(v() As Double, Byref result As Double)
result = 0
For i As Integer = Lbound(v) To Ubound(v)
result += v(i) * v(i)
Next i
result = Sqr(result)
End Sub
 
Sub mulAdd(v1() As Double, x1 As Double, v2() As Double, x2 As Double, result() As Double)
For i As Integer = Lbound(v1) To Ubound(v1)
result(i) = v1(i) * x1 + v2(i) * x2
Next i
End Sub
 
Sub rotate(i() As Double, j() As Double, alfa As Double, result1() As Double, result2() As Double)
Dim As Double ca = Cos(alfa), sa = Sin(alfa)
mulAdd(i(), ca, j(), sa, result1())
mulAdd(i(), -sa, j(), ca, result2())
End Sub
 
Sub orbitalStateVectors(semimajorAxis As Double, eccentricity As Double, inclination As Double, longitudeOfAscendingNode As Double, argumentOfPeriapsis As Double, trueAnomaly As Double)
Dim As Double i(1 To 3) = {1, 0, 0}
Dim As Double j(1 To 3) = {0, 1, 0}
Dim As Double k(1 To 3) = {0, 0, 1}
Dim As Double temp1(1 To 3), temp2(1 To 3)
Dim As Integer index, t
rotate(i(), j(), longitudeOfAscendingNode, temp1(), temp2())
For index = 1 To 3
i(index) = temp1(index)
j(index) = temp2(index)
Next index
rotate(j(), k(), inclination, temp1(), temp2())
For index = 1 To 3
j(index) = temp1(index)
Next index
rotate(i(), j(), argumentOfPeriapsis, temp1(), temp2())
For index = 1 To 3
i(index) = temp1(index)
j(index) = temp2(index)
Next index
Dim As Double l = Iif(eccentricity = 1, 2, 1 - eccentricity * eccentricity) * semimajorAxis
Dim As Double c = Cos(trueAnomaly), s = Sin(trueAnomaly)
Dim As Double r = 1 / (1 + eccentricity * c)
Dim As Double rprime = s * r * r / l
Dim As Double posn(1 To 3), speed(1 To 3), vabsResult
mulAdd(i(), c, j(), s, posn())
For t = Lbound(posn) To Ubound(posn)
posn(t) *= r
Next t
mulAdd(i(), rprime * c - r * s, j(), rprime * s + r * c, speed())
vabs(speed(), vabsResult)
mulAdd(speed(), 1 / vabsResult, speed(), 0, speed())
For t = Lbound(speed) To Ubound(speed)
speed(t) *= Sqr(2 / r - 1 / semimajorAxis)
Next t
Print Using "Position : {&, &, &}"; posn(1); posn(2); posn(3)
Print Using "Speed : {&, &, &}"; speed(1); speed(2); speed(3)
End Sub
 
orbitalStateVectors(1.0, 0.1, 0.0, 355.0 / (113.0 * 6.0), 0.0, 0.0)
 
Sleep</syntaxhighlight>
{{out}}
<pre>Position : {0.7872958014128079, 0.454545489549176, 0}
Speed : {-0.5477225996842874, 0.9486832736983857, 0}</pre>
 
=={{header|Go}}==
{{trans|Kotlin}}
<langsyntaxhighlight lang="go">package main
 
import (
Line 731 ⟶ 862:
fmt.Println("Position :", position)
fmt.Println("Speed :", speed)
}</langsyntaxhighlight>
 
{{out}}
Line 740 ⟶ 871:
 
=={{header|J}}==
{{trans|Raku}}<langsyntaxhighlight Jlang="j">NB. euler rotation matrix, left hand rule
NB. x: axis (0, 1 or 2), y: angle in radians
R=: {{ ((2 1,:1 2) o.(,-)y*_1^2|x)(,&.>/~0 1 2-.x)} =i.3 }}
Line 752 ⟶ 883:
NB. Om: Longitude of the ascending node
NB. w: argument of Periapsis (the other "omega")
NB. f: true anomaly (varies with time)
L=. a*2:`]@.*1-*:e
'c s'=. 2{.,F=. 2 R f
Line 761 ⟶ 892:
speed=. (%:(2%ra)-%a)*norm(rp,ra,0) X ijk
position,:speed
}}</langsyntaxhighlight>
 
The true anomaly, argument of Periapsis, Longitude of the ascending node and inclination are all angles. And we use the dot product of their rotation matrices (in that order) to find the orientation of the orbit and the object's position in that orbit. Here, <code>R</code> finds the rotation matrix for a given angle around a given axis. Here's an example of what R gives us for a sixty degree angle:
 
<syntaxhighlight lang="j"> 0 1 2 R&.> 60r180p1 NB. rotate around first, second or third axis
┌────────────────────┬────────────────────┬────────────────────┐
│1 0 0│ 0.5 0 _0.866025│ 0.5 0.866025 0│
│0 0.5 0.866025│ 0 1 0│_0.866025 0.5 0│
│0 _0.866025 0.5│0.866025 0 0.5│ 0 0 1│
└────────────────────┴────────────────────┴────────────────────┘</syntaxhighlight>
 
Task example:<langsyntaxhighlight Jlang="j"> orbitalStateVectors 1 0.1 0 355r678 0 0
0.779423 0.45 0
_0.552771 0.957427 0</langsyntaxhighlight>
 
=={{header|Java}}==
{{trans|Kotlin}}
<langsyntaxhighlight Javalang="java">public class OrbitalElements {
private static class Vector {
private double x, y, z;
Line 849 ⟶ 989:
System.out.printf("Speed : %s\n", ps[1]);
}
}</langsyntaxhighlight>
 
{{out}}
Line 859 ⟶ 999:
{{works with|jq}}
'''Works with gojq, the Go implementation of jq'''
<langsyntaxhighlight lang="jq"># Array/vector operations
def addVectors: transpose | map(add);
 
Line 895 ⟶ 1,035:
| divide(abs)
| multiply( ((2 / $r) - (1 / semimajorAxis))|sqrt) as $speed
| [$position, $speed] ;</langsyntaxhighlight>
'''The Task'''
<langsyntaxhighlight lang="jq">orbitalStateVectors(1; 0.1; 0; 355 / (113 * 6); 0; 0)
| "Position : \(.[0])",
"Speed : \(.[1])"</langsyntaxhighlight>
{{out}}
<pre>
Line 908 ⟶ 1,048:
=={{header|Julia}}==
{{trans|Kotlin}}
<langsyntaxhighlight lang="julia">using GeometryTypes
import Base.abs, Base.print
 
Line 943 ⟶ 1,083:
 
testorbitalmath()
</langsyntaxhighlight>{{out}}
<pre>
Position : (0.7794228433986797, 0.45000003465368416, 0.0)
Line 951 ⟶ 1,091:
=={{header|Kotlin}}==
{{trans|Sidef}}
<langsyntaxhighlight lang="scala">// version 1.1.4-3
 
class Vector(val x: Double, val y: Double, val z: Double) {
Line 1,014 ⟶ 1,154:
println("Position : $position")
println("Speed : $speed")
}</langsyntaxhighlight>
{{out}}
<pre>Position : (0.7794228433986797, 0.45000003465368416, 0.0)
Speed : (-0.5527708409604438, 0.9574270831797618, 0.0)</pre>
 
=={{header|Lua}}==
{{Trans|C}}...which is translation of Kotlin which is ...
<syntaxhighlight lang="lua">
do -- orbital elements
 
local function Vector( x, y, z )
return { x = x, y= y, z = z }
end
 
local function add( v, w )
return Vector( v.x + w.x, v.y + w.y, v.z + w.z )
end
 
local function mul( v, m )
return Vector( v.x * m, v.y * m, v.z * m )
end
 
local function div( v, d )
return mul( v, 1.0 / d )
end
 
local function vabs( v )
return math.sqrt( v.x * v.x + v.y * v.y + v.z * v.z )
end
 
local function mulAdd( v1, v2, x1, x2 )
return add( mul( v1, x1 ), mul( v2, x2 ) )
end
 
local function vecAsStr( v )
return string.format( "(%.17g", v.x )..string.format( ", %.17g", v.y )..string.format( ", %.17g)", v.z )
end
 
local function rotate( i, j, alpha )
return mulAdd( i, j, math.cos( alpha ), math.sin( alpha ) )
, mulAdd( i, j, -math.sin( alpha ), math.cos( alpha ) )
end
local function orbitalStateVectors( semimajorAxis, eccentricity, inclination
, longitudeOfAscendingNode, argumentOfPeriapsis
, trueAnomaly
)
 
local i, j, k = Vector( 1.0, 0.0, 0.0 ), Vector( 0.0, 1.0, 0.0 ), Vector( 0.0, 0.0, 1.0 )
local L = 2.0
i, j = rotate( i, j, longitudeOfAscendingNode )
j, _ = rotate( j, k, inclination )
i, j = rotate( i, j, argumentOfPeriapsis )
if eccentricity ~= 1.0 then L = 1.0 - eccentricity * eccentricity end
L = L * semimajorAxis
local c, s = math.cos( trueAnomaly ), math.sin( trueAnomaly )
local r = L / ( 1.0 + eccentricity * c )
local rprime = s * r * r / L;
local position = mul( mulAdd( i, j, c, s ), r )
local speed = mulAdd( i, j, rprime * c - r * s, rprime * s + r * c )
speed = div( speed, vabs( speed ) )
speed = mul( speed, math.sqrt( 2.0 / r - 1.0 / semimajorAxis ) )
return position, speed
end
 
local longitude = 355.0 / ( 113.0 * 6.0 )
local position, speed = orbitalStateVectors( 1.0, 0.1, 0.0, longitude, 0.0, 0.0 )
print( "Position : "..vecAsStr( position ) )
print( "Speed : "..vecAsStr( speed ) )
end</syntaxhighlight>
{{out}}
<pre>
Position : (0.77942284339867973, 0.45000003465368416, 0)
Speed : (-0.55277084096044382, 0.95742708317976177, 0)
</pre>
 
=={{header|Nim}}==
{{trans|Kotlin}}
<langsyntaxhighlight Nimlang="nim">import math, strformat
 
type Vector = tuple[x, y, z: float]
Line 1,077 ⟶ 1,288:
trueAnomaly = 0.0)
echo "Position: ", position
echo "Speed: ", speed</langsyntaxhighlight>
 
{{out}}
Line 1,085 ⟶ 1,296:
=={{header|ooRexx}}==
{{trans|Java}}
<langsyntaxhighlight lang="oorexx">/* REXX */
Numeric Digits 16
ps = orbitalStateVectors(1.0, 0.1, 0.0, 355.0 / (113.0 * 6.0), 0.0, 0.0)
Line 1,191 ⟶ 1,402:
Return res
 
::requires 'rxmath' LIBRARY</langsyntaxhighlight>
 
{{out}}
Line 1,202 ⟶ 1,413:
=={{header|Perl}}==
{{trans|Raku}}
<langsyntaxhighlight lang="perl">use strict;
use warnings;
use Math::Vector::Real;
Line 1,261 ⟶ 1,472:
0, # argument of periapsis
0 # true-anomaly
;</langsyntaxhighlight>
{{out}}
<pre>$VAR1 = {
Line 1,278 ⟶ 1,489:
=={{header|Phix}}==
{{trans|Python}}
<!--<langsyntaxhighlight Phixlang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">vabs</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span>
Line 1,318 ⟶ 1,529:
<span style="color: #000000;">orbitalStateVectors</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">355.0</span> <span style="color: #0000FF;">/</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">113.0</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">6.0</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">0.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.0</span><span style="color: #0000FF;">)</span>
<!--</langsyntaxhighlight>-->
{{out}}
<pre>
Line 1,330 ⟶ 1,541:
This implementation uses the CLP/R library of swi-prolog, but doesn't have to. This removes the need for a vector divide and has limited capability to reverse the functionality (eg: given the position/speed find some orbital elements).
 
<langsyntaxhighlight Prologlang="prolog">:- use_module(library(clpr)).
 
v3_add(v(X1,Y1,Z1),v(X2,Y2,Z2),v(X,Y,Z)) :-
Line 1,388 ⟶ 1,599:
find_l(Ecc, SemiMajor, L) :-
dif(Ecc,1.0),
{ L = SemiMajor * (1.0 - Ecc * Ecc) }.</langsyntaxhighlight>
{{out}}
<pre>
Line 1,399 ⟶ 1,610:
 
=={{header|Python}}==
<langsyntaxhighlight lang="python">import math
 
class Vector:
Line 1,457 ⟶ 1,668:
ps = orbitalStateVectors(1.0, 0.1, 0.0, 355.0 / (113.0 * 6.0), 0.0, 0.0)
print "Position :", ps[0]
print "Speed :", ps[1]</langsyntaxhighlight>
{{out}}
<pre>Position : (0.787295801413, 0.454545489549, 0.0)
Line 1,465 ⟶ 1,676:
(formerly Perl 6)
We'll use the [https://github.com/grondilu/clifford Clifford geometric algebra library] but only for the vector operations.
<syntaxhighlight lang="raku" perl6line>sub orbital-state-vectors(
Real :$semimajor-axis where * >= 0,
Real :$eccentricity where * >= 0,
Line 1,508 ⟶ 1,719:
longitude-of-ascending-node => pi/6,
argument-of-periapsis => pi/4,
true-anomaly => 0;</langsyntaxhighlight>
{{out}}
<pre>{position => 0.237771283982207*e0+0.860960261697716*e1+0.110509023572076*e2, speed => -1.06193301748006*e0+0.27585002056925*e1+0.135747024865598*e2}</pre>
Line 1,516 ⟶ 1,727:
{{trans|Java}}
Vectors are represented by strings: 'x/y/z'
<langsyntaxhighlight lang="rexx">/* REXX */
Numeric Digits 16
Parse Value orbitalStateVectors(1.0,0.1,0.0,355.0/(113.0*6.0),0.0,0.0),
Line 1,651 ⟶ 1,862:
End
Numeric Digits prec
Return r+0</langsyntaxhighlight>
{{out}}
<pre>Position : (0.7794228433986798,0.4500000346536842,0)
Line 1,658 ⟶ 1,869:
===version 2===
Re-coding of REXX version 1, &nbsp; but with greater decimal digits precision.
<langsyntaxhighlight lang="rexx">/*REXX pgm converts orbital elements ──► orbital state vectors (angles are in radians).*/
numeric digits length( pi() ) - length(.) /*limited to pi len, but show 1/3 digs.*/
call orbV 1, .1, 0, 355/113/6, 0, 0 /*orbital elements taken from: Java */
Line 1,709 ⟶ 1,920:
numeric digits; parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g= g *.5'e'_ % 2
do j=0 while h>9; m.j= h; h= h % 2 + 1; end
do k=j+5 to 0 by '-1'; numeric digits m.k; g= (g+x/g) * .5; end; return g</langsyntaxhighlight>
{{out|output|text=&nbsp; when using the default internal inputs:}}
<pre>
Line 1,734 ⟶ 1,945:
 
=={{header|Scala}}==
<langsyntaxhighlight Scalalang="scala">import scala.language.existentials
 
object OrbitalElements extends App {
Line 1,776 ⟶ 1,987:
}
 
}</langsyntaxhighlight>
{{Out}}Best seen running in your browser either by [https://scalafiddle.io/sf/ac17jh2/0 ScalaFiddle (ES aka JavaScript, non JVM)] or [https://scastie.scala-lang.org/2NQNgj4OQkazxZNvSzcexQ Scastie (remote JVM)].
 
=={{header|Sidef}}==
{{trans|Perl}}
<langsyntaxhighlight lang="ruby">func orbital_state_vectors(
semimajor_axis,
eccentricity,
Line 1,837 ⟶ 2,048:
say "Position : #{r.position}"
say "Speed : #{r.speed}\n"
}</langsyntaxhighlight>
{{out}}
<pre>
Line 1,853 ⟶ 2,064:
{{trans|Kotlin}}
 
<langsyntaxhighlight lang="swift">import Foundation
 
public struct Vector {
Line 1,955 ⟶ 2,166:
)
 
print("Position: \(position); Speed: \(speed)")</langsyntaxhighlight>
 
{{out}}
Line 1,963 ⟶ 2,174:
=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-vector}}
<lang ecmascript>class Vector {
<syntaxhighlight lang="wren">import "./vector" for Vector3
construct new(x, y, z) {
_x = x
_y = y
_z = z
}
 
x { _x }
y { _y }
z { _z }
 
+(other) { Vector.new(_x + other.x, _y + other.y, _z + other.z) }
*(m) { Vector.new(_x * m, _y * m, _z * m) }
 
/(d) { this * (1/d) }
 
abs { (_x *_x + _y *_y + _z * _z).sqrt }
 
toString { "(%(_x), %(_y), %(_z))" }
}
 
var orbitalStateVectors = Fn.new { |semimajorAxis, eccentricity, inclination,
longitudeOfAscendingNode, argumentOfPeriapsis, trueAnomaly|
var i = VectorVector3.new(1, 0, 0)
var j = VectorVector3.new(0, 1, 0)
var k = VectorVector3.new(0, 0, 1)
 
var mulAdd = Fn.new { |v1, x1, v2, x2| v1 * x1 + v2 * x2 }
Line 2,014 ⟶ 2,206:
var position = mulAdd.call(i, c, j, s) * r
var speed = mulAdd.call(i, rprime * c - r * s, j, rprime * s + r * c)
speed = speed / speed.abslength
speed = speed * (2 / r - 1 / semimajorAxis).sqrt
return [position, speed]
Line 2,021 ⟶ 2,213:
var ps = orbitalStateVectors.call(1, 0.1, 0, 355 / (113 * 6), 0, 0)
System.print("Position : %(ps[0])")
System.print("Speed : %(ps[1])")</langsyntaxhighlight>
 
{{out}}
Line 2,027 ⟶ 2,219:
Position : (0.77942284339868, 0.45000003465368, 0)
Speed : (-0.55277084096044, 0.95742708317976, 0)
</pre>
 
=={{header|XPL0}}==
{{trans|C}}
<syntaxhighlight lang "XPL0">include xpllib; \for VCopy, VAdd, VMul, VDiv, VMag and Print
 
proc VShow(A);
real A;
Print("(%1.6f, %1.6f, %1.6f)\n", A(0), A(1), A(2));
 
func real VMulAdd(V0, V1, V2, X1, X2);
real V0, V1, V2, X1, X2, V3(3), V4(3);
return VAdd(V0, VMul(V3, V1, X1), VMul(V4, V2, X2));
 
proc Rotate(I, J, Alpha, PS);
real I, J, Alpha, PS;
[VMulAdd(PS(0), I, J, Cos(Alpha), Sin(Alpha));
VMulAdd(PS(1), I, J, -Sin(Alpha), Cos(Alpha));
];
 
proc OrbitalStateVectors; real SemiMajorAxis, Eccentricity, Inclination,
LongitudeOfAscendingNode, ArgumentOfPeriapsis, TrueAnomaly, PS;
real I, J, K, L, QS(2,3), C, S, R, RPrime;
[I:= [1.0, 0.0, 0.0];
J:= [0.0, 1.0, 0.0];
K:= [0.0, 0.0, 1.0];
L:= 2.0;
Rotate(I, J, LongitudeOfAscendingNode, QS);
VCopy(I, QS(0)); VCopy(J, QS(1));
Rotate(J, K, Inclination, QS);
VCopy(J, QS(0));
Rotate(I, J, ArgumentOfPeriapsis, QS);
VCopy(I, QS(0)); VCopy(J, QS(1));
if Eccentricity # 1.0 then L:= 1.0 - Eccentricity*Eccentricity;
L:= L * SemiMajorAxis;
C:= Cos(TrueAnomaly);
S:= Sin(TrueAnomaly);
R:= L / (1.0 + Eccentricity*C);
RPrime:= S * R * R / L;
VMulAdd(PS(0), I, J, C, S);
VMul(PS(0), PS(0), R);
VMulAdd(PS(1), I, J, RPrime*C - R*S, RPrime*S + R*C);
VDiv(PS(1), PS(1), VMag(PS(1)));
VMul(PS(1), PS(1), sqrt(2.0/R - 1.0/SemiMajorAxis));
];
 
def Longitude = 355.0 / (113.0 * 6.0);
real PS(2,3);
[OrbitalStateVectors(1.0, 0.1, 0.0, Longitude, 0.0, 0.0, PS);
Print("Position : "); VShow(PS(0));
Print("Speed : "); VShow(PS(1));
]</syntaxhighlight>
{{out}}
<pre>
Position : (0.779423, 0.450000, 0.000000)
Speed : (-0.552771, 0.957427, 0.000000)
</pre>
 
=={{header|zkl}}==
{{trans|Perl}}
<langsyntaxhighlight 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);
Line 2,060 ⟶ 2,308:
return(position,speed);
}</langsyntaxhighlight>
<langsyntaxhighlight lang="zkl">orbital_state_vectors(
1.0, # semimajor axis
0.1, # eccentricity
Line 2,068 ⟶ 2,316:
0.0, # argument of periapsis
0.0 # true-anomaly
).println();</langsyntaxhighlight>
{{out}}
<pre>L(L(0.779423,0.45,0),L(-0.552771,0.957427,0))</pre>
Jjuanhdez
2,123

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/268822...362722"