Orbital elements: Difference between revisions
Added FreeBASIC
(link to Wikipedia article about polar coordinates) |
(Added FreeBASIC) |
||
| (5 intermediate revisions by 3 users not shown) | |||
Line 189:
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>
Line 657 ⟶ 717:
<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}}==
Line 1,027 ⟶ 1,158:
<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}}==
Line 1,972 ⟶ 2,174:
=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-vector}}
<syntaxhighlight lang="wren">import "./vector" for Vector3
var orbitalStateVectors = Fn.new { |semimajorAxis, eccentricity, inclination,
longitudeOfAscendingNode, argumentOfPeriapsis, trueAnomaly|
var i =
var j =
var k =
var mulAdd = Fn.new { |v1, x1, v2, x2| v1 * x1 + v2 * x2 }
Line 2,023 ⟶ 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.
speed = speed * (2 / r - 1 / semimajorAxis).sqrt
return [position, speed]
Line 2,036 ⟶ 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>
| |||