Skip to content
This repository was archived by the owner on Aug 5, 2024. It is now read-only.

Commit ba0dfa8

Browse files
sylvain-guilletsylvain-guillet
authored
Fix/orbitalparameters (#99)
* Define convention for special cases
1 parent 8c88897 commit ba0dfa8

2 files changed

Lines changed: 67 additions & 34 deletions

File tree

‎IO.Astrodynamics.Tests/OrbitalParameters/StateVectorTests.cs‎

Lines changed: 47 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@ public StateVectorTests()
1919
[Fact]
2020
public void Create()
2121
{
22-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
22+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
2323

2424
Vector3 pos = new Vector3(1.0, 2.0, 3.0);
2525
Vector3 vel = new Vector3(4.0, 5.0, 6.0);
@@ -35,7 +35,7 @@ public void Create()
3535
[Fact]
3636
public void Inverse()
3737
{
38-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
38+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
3939

4040
Vector3 pos = new Vector3(1.0, 2.0, 3.0);
4141
Vector3 vel = new Vector3(4.0, 5.0, 6.0);
@@ -48,7 +48,7 @@ public void Inverse()
4848
[Fact]
4949
public void Add()
5050
{
51-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
51+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
5252

5353
Vector3 pos = new Vector3(1.0, 2.0, 3.0);
5454
Vector3 vel = new Vector3(4.0, 5.0, 6.0);
@@ -65,7 +65,7 @@ public void Add()
6565
[Fact]
6666
public void AddExcept()
6767
{
68-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
68+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
6969

7070
Vector3 pos = new Vector3(1.0, 2.0, 3.0);
7171
Vector3 vel = new Vector3(4.0, 5.0, 6.0);
@@ -81,7 +81,7 @@ public void AddExcept()
8181
[Fact]
8282
public void Subtract()
8383
{
84-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
84+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
8585

8686
Vector3 pos = new Vector3(1.0, 2.0, 3.0);
8787
Vector3 vel = new Vector3(4.0, 5.0, 6.0);
@@ -95,7 +95,7 @@ public void Subtract()
9595
[Fact]
9696
public void SubtractExcept()
9797
{
98-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
98+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
9999

100100
Vector3 pos = new Vector3(1.0, 2.0, 3.0);
101101
Vector3 vel = new Vector3(4.0, 5.0, 6.0);
@@ -111,17 +111,17 @@ public void SubtractExcept()
111111
[Fact]
112112
public void Eccentricity()
113113
{
114-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
114+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
115115

116116
StateVector sv = new StateVector(new Vector3(-6.116559469556896E+06, -1.546174698676721E+06, 2.521950157430313E+06),
117117
new Vector3(-8.078523150700097E+02, -5.477647950892673E+03, -5.297615757935174E+03), earth, DateTime.UtcNow, Frames.Frame.ICRF);
118-
Assert.Equal(1.3532176446914895E-03, sv.Eccentricity(),6);
118+
Assert.Equal(1.3532176446914895E-03, sv.Eccentricity(), 6);
119119
}
120120

121121
[Fact]
122122
public void EccentricityVector()
123123
{
124-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
124+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
125125

126126
StateVector sv = new StateVector(new Vector3(6800000.0, 0.0, 0.0), new Vector3(0.0, 8000.0, 0.0), earth, DateTime.UtcNow, Frames.Frame.ICRF);
127127
Vector3 ev = sv.EccentricityVector();
@@ -134,7 +134,7 @@ public void EccentricityVector()
134134
[Fact]
135135
public void SpecificAngularMomentum()
136136
{
137-
CelestialBody sun = new CelestialBody( Stars.Sun);
137+
CelestialBody sun = new CelestialBody(Stars.Sun);
138138

139139
StateVector sv = new StateVector(new Vector3(149600000.0, 0.0, 0.0), new Vector3(0.0, 29.8, 0.0), sun, DateTime.UtcNow, Frames.Frame.ICRF);
140140

@@ -147,7 +147,7 @@ public void SpecificAngularMomentum()
147147
[Fact]
148148
public void SpecificOrbitalEnergyMomentum()
149149
{
150-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
150+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
151151

152152
StateVector sv = new StateVector(new Vector3(-6.116559469556896E+06, -1.546174698676721E+06, 2.521950157430313E+06),
153153
new Vector3(-8.078523150700097E+02, -5.477647950892673E+03, -5.297615757935174E+03), earth, DateTime.UtcNow, Frames.Frame.ICRF);
@@ -158,7 +158,7 @@ public void SpecificOrbitalEnergyMomentum()
158158
[Fact]
159159
public void Inclination()
160160
{
161-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
161+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
162162

163163
StateVector sv = new StateVector(new Vector3(6800.0, 0.0, 0.0), new Vector3(0.0, 5.0, 5.0), earth, DateTime.UtcNow, Frames.Frame.ICRF);
164164
Assert.Equal(System.Math.PI / 4.0, sv.Inclination());
@@ -167,7 +167,7 @@ public void Inclination()
167167
[Fact]
168168
public void SemiMajorAxis()
169169
{
170-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
170+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
171171

172172
StateVector sv = new StateVector(new Vector3(8000000.0, 0.0, 0.0), new Vector3(0.0, 6000.0, 6000.0), earth, DateTime.UtcNow, Frames.Frame.ICRF);
173173
Assert.Equal(14415872.186388023, sv.SemiMajorAxis());
@@ -176,7 +176,7 @@ public void SemiMajorAxis()
176176
[Fact]
177177
public void AscendingNodeVector()
178178
{
179-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
179+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
180180

181181
StateVector sv = new StateVector(new Vector3(8000000.0, 0.0, 0.0), new Vector3(0.0, 6000.0, 0.0), earth, DateTime.UtcNow, Frames.Frame.ICRF);
182182
var v = sv.AscendingNodeVector().Normalize();
@@ -188,7 +188,7 @@ public void AscendingNodeVector()
188188
[Fact]
189189
public void AscendingNode()
190190
{
191-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
191+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
192192

193193
StateVector sv = new StateVector(new Vector3(9208000.0, 3352000, 0.0), new Vector3(-1750, 4830, 5140), earth, DateTime.UtcNow, Frames.Frame.ICRF);
194194
Assert.Equal(20.00308830929978, sv.AscendingNode() * IO.Astrodynamics.Constants.Rad2Deg);
@@ -200,78 +200,78 @@ public void AscendingNode()
200200
[Fact]
201201
public void ArgumentOfPeriapis()
202202
{
203-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
203+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
204204

205205
StateVector sv = new StateVector(new Vector3(8237000.0, 17000.0, 5308000.0), new Vector3(-2000.0, 6000.0, 3000.0), earth, DateTime.UtcNow, Frames.Frame.ICRF);
206206

207-
Assert.Equal(53.621820299859436, sv.ArgumentOfPeriapsis() * IO.Astrodynamics.Constants.Rad2Deg,3);
207+
Assert.Equal(53.621820299859436, sv.ArgumentOfPeriapsis() * IO.Astrodynamics.Constants.Rad2Deg, 3);
208208

209209
sv = new StateVector(new Vector3(3973000.0, -4881000.0, -931000.0), new Vector3(5400.0, 3360.0, 5890.0), earth, DateTime.UtcNow, Frames.Frame.ICRF);
210210

211-
Assert.Equal(350.51061587951995, sv.ArgumentOfPeriapsis() * IO.Astrodynamics.Constants.Rad2Deg,3);
211+
Assert.Equal(350.51061587951995, sv.ArgumentOfPeriapsis() * IO.Astrodynamics.Constants.Rad2Deg, 3);
212212
}
213213

214214
[Fact]
215215
public void TrueAnomaly()
216216
{
217-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
217+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
218218

219219
StateVector sv = new StateVector(new Vector3(5070000.0, -2387000.0, 1430000.0), new Vector3(2450.0, 6350.0, 6440.0), earth, DateTime.UtcNow, Frames.Frame.ICRF);
220220

221-
Assert.Equal(30.57538215436613, sv.TrueAnomaly() * IO.Astrodynamics.Constants.Rad2Deg,6);
221+
Assert.Equal(30.57538215436613, sv.TrueAnomaly() * IO.Astrodynamics.Constants.Rad2Deg, 6);
222222

223223
sv = new StateVector(new Vector3(1664000.0, -4862000.0, -2655000.0), new Vector3(7520.0, 890, 5520.0), earth, DateTime.UtcNow, Frames.Frame.ICRF);
224224

225-
Assert.Equal(329.477000818209, sv.TrueAnomaly() * IO.Astrodynamics.Constants.Rad2Deg,6);
225+
Assert.Equal(329.477000818209, sv.TrueAnomaly() * IO.Astrodynamics.Constants.Rad2Deg, 6);
226226
}
227227

228228
[Fact]
229229
public void EccentricAnomaly()
230230
{
231-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
231+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
232232

233233
StateVector sv = new StateVector(new Vector3(6700000.0, 2494000.0, 0.0), new Vector3(-2150.0, 8850.0, 0.0), earth, DateTime.UtcNow, Frames.Frame.ICRF);
234234

235-
Assert.Equal(0.20776173316152752, sv.EccentricAnomaly(),6);
235+
Assert.Equal(0.20776173316152752, sv.EccentricAnomaly(), 6);
236236
}
237237

238238
[Fact]
239239
public void MeanAnomaly()
240240
{
241-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
241+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
242242

243243
StateVector sv = new StateVector(new Vector3(-5775.068936894231E+03, -3372.353197848874E+03, 651.695854037289E+03),
244244
new Vector3(-0.661469579672604E+03, -7.147573777688288E+03, -2.915719736461653E+03), earth, DateTime.UtcNow, Frames.Frame.ICRF);
245245

246-
Assert.Equal(59.99823497591805, sv.MeanAnomaly() * IO.Astrodynamics.Constants.Rad2Deg,3);
246+
Assert.Equal(59.99823497591805, sv.MeanAnomaly() * IO.Astrodynamics.Constants.Rad2Deg, 3);
247247
}
248248

249249
[Fact]
250250
public void Period()
251251
{
252-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
252+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
253253

254254
StateVector sv = new StateVector(new Vector3(-5775.068936894231E+03, -3372.353197848874E+03, 651.695854037289E+03),
255255
new Vector3(-0.661469579672604E+03, -7.147573777688288E+03, -2.915719736461653E+03), earth, DateTime.UtcNow, Frames.Frame.ICRF);
256256

257-
Assert.Equal(1.5501796754722221, sv.Period().TotalHours,6);
257+
Assert.Equal(1.5501796754722221, sv.Period().TotalHours, 6);
258258
}
259259

260260
[Fact]
261261
public void MeanMotion()
262262
{
263-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
263+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
264264

265265
StateVector sv = new StateVector(new Vector3(-5775.068936894231E+03, -3372.353197848874E+03, 651.695854037289E+03),
266266
new Vector3(-0.661469579672604E+03, -7.147573777688288E+03, -2.915719736461653E+03), earth, DateTime.UtcNow, Frames.Frame.ICRF);
267267

268-
Assert.Equal(0.0011258883596591215, sv.MeanMotion(),3);
268+
Assert.Equal(0.0011258883596591215, sv.MeanMotion(), 3);
269269
}
270270

271271
[Fact]
272272
public void ToFrame()
273273
{
274-
CelestialBody sun = new CelestialBody( Stars.Sun);
274+
CelestialBody sun = new CelestialBody(Stars.Sun);
275275

276276
//J2000->Ecliptic
277277
//Earth from sun at 0 TDB
@@ -290,7 +290,7 @@ public void ToFrame()
290290
[Fact]
291291
public void ToNonInertialFrame()
292292
{
293-
CelestialBody earth = new CelestialBody( PlanetsAndMoons.EARTH);
293+
CelestialBody earth = new CelestialBody(PlanetsAndMoons.EARTH);
294294

295295
var epoch = new DateTime(2000, 1, 1, 12, 0, 0);
296296
var earthFrame = new Frames.Frame(PlanetsAndMoons.EARTH.Frame);
@@ -355,6 +355,19 @@ public void ToKeplerian()
355355
1.77688489436688, 6.259056257653703, TestHelpers.Sun, DateTimeExtension.J2000, Frames.Frame.ICRF), ke);
356356
}
357357

358+
[Fact]
359+
public void ToKeplerianSingularity()
360+
{
361+
var svOrigin =
362+
new StateVector(new Vector3(0, 6700000, 0), new Vector3(7900, 0, 0), TestHelpers.EarthAtJ2000, DateTimeExtension.J2000, Frames.Frame.ICRF);
363+
364+
var keOrigin = svOrigin.ToKeplerianElements();
365+
Assert.Equal(new KeplerianElements(7045497.0180164594, 0.049037990809302601, 3.1415926535897931, 0.0,
366+
4.7123889803846897, 0.0, TestHelpers.EarthAtJ2000, DateTimeExtension.J2000, Frames.Frame.ICRF), keOrigin);
367+
var svFinal = keOrigin.ToStateVector();
368+
Assert.Equal(svFinal,svOrigin,TestHelpers.StateVectorComparer);
369+
}
370+
358371
[Fact]
359372
public void ToEquatorial()
360373
{
@@ -368,7 +381,9 @@ public void RelativeTo()
368381
{
369382
var originalSV = new StateVector(new Vector3(6800000.0, 0.0, 0.0), new Vector3(0.0, 8000.0, 0.0), TestHelpers.EarthAtJ2000, DateTimeExtension.J2000, Frames.Frame.ICRF);
370383
var moonSv = originalSV.RelativeTo(TestHelpers.MoonAtJ2000, Aberration.None);
371-
Assert.Equal(new StateVector(new Vector3(298408384.63343549, 266716833.39423338, 76102487.099902019), new Vector3(-643.53138771903275, 8666.0876840916299, 301.32570498227307), TestHelpers.MoonAtJ2000, DateTimeExtension.J2000, Frames.Frame.ICRF),
384+
Assert.Equal(
385+
new StateVector(new Vector3(298408384.63343549, 266716833.39423338, 76102487.099902019), new Vector3(-643.53138771903275, 8666.0876840916299, 301.32570498227307),
386+
TestHelpers.MoonAtJ2000, DateTimeExtension.J2000, Frames.Frame.ICRF),
372387
moonSv);
373388
}
374389
}

‎IO.Astrodynamics/OrbitalParameters/StateVector.cs‎

Lines changed: 20 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -107,7 +107,15 @@ public override double AscendingNode()
107107

108108
Vector3 n = AscendingNodeVector();
109109

110+
if (n.Magnitude() == 0.0)
111+
{
112+
_ascendingNode = 0.0;
113+
return _ascendingNode.Value;
114+
}
115+
110116
var omega = System.Math.Acos(n.X / n.Magnitude());
117+
118+
111119
if (n.Y < 0.0)
112120
{
113121
omega = 2 * System.Math.PI - omega;
@@ -133,8 +141,18 @@ public override double ArgumentOfPeriapsis()
133141
{
134142
return 0.0;
135143
}
136-
137-
_argumentOfPeriapsis = System.Math.Acos(System.Math.Clamp((n * e) / (n.Magnitude() * e.Magnitude()),-1.0,1.0));
144+
145+
if (n == Vector3.Zero)
146+
{
147+
_argumentOfPeriapsis = System.Math.Atan2(e.Y, e.X);
148+
if (Inclination() > Constants.PI2)
149+
{
150+
_argumentOfPeriapsis = Constants._2PI - _argumentOfPeriapsis;
151+
return _argumentOfPeriapsis!.Value;
152+
}
153+
}
154+
155+
_argumentOfPeriapsis = System.Math.Acos(System.Math.Clamp((n * e) / (n.Magnitude() * e.Magnitude()), -1.0, 1.0));
138156
if (e.Z < 0.0)
139157
{
140158
_argumentOfPeriapsis = System.Math.PI * 2.0 - _argumentOfPeriapsis;

0 commit comments

Comments
 (0)