/
GoodingTests.cs
89 lines (70 loc) · 2.91 KB
/
GoodingTests.cs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
using System.Collections.Generic;
using AssertExtensions;
using MechJebLib.Core;
using Xunit;
using Xunit.Abstractions;
using static MechJebLib.Utils.Statics;
using static System.Math;
namespace MechJebLibTest.Maths
{
public class GoodingTests
{
private readonly ITestOutputHelper _testOutputHelper;
public GoodingTests(ITestOutputHelper testOutputHelper)
{
_testOutputHelper = testOutputHelper;
}
[Fact]
private void SingleRevolution()
{
double mu = 1.0;
for (int j = 0; j < 40; j++)
{
double ecc = 0.1 * j;
double sma = ecc > 1 ? -1.0 : 1.0;
_testOutputHelper.WriteLine($"{ecc}");
var elist = new List<double>(); // eccentric anomaly
var tlist = new List<double>(); // time of flight
var rlist = new List<double>(); // magnitude of r
var vlist = new List<double>(); // mangitude of v
var flist = new List<double>(); // true anomaly
for (int i = 0; i < 360; i += 4)
{
double eanom = Deg2Rad(i);
double time = MechJebLib.Core.Maths.TimeSincePeriapsisFromEccentricAnomaly(mu, sma, ecc, eanom);
double tanom = MechJebLib.Core.Maths.TrueAnomalyFromEccentricAnomaly(ecc, eanom);
double smp = ecc == 1 ? 2 * sma : sma * (1.0 - ecc * ecc);
double energy = ecc != 1 ? -1.0 / (2.0 * sma) : 0;
double r = smp / (1.0 + ecc * Cos(tanom));
double v = Sqrt(2 * (energy + 1.0 / r));
elist.Add(eanom);
tlist.Add(time);
rlist.Add(r);
vlist.Add(v);
flist.Add(tanom);
}
double diffmax = 0;
for (int n1 = 0; n1 < elist.Count; n1++)
{
for (int n2 = n1 + 1; n2 < elist.Count; n2++)
{
double VR11, VT11, VR12, VT12;
(_, VR11, VT11, VR12, VT12, _, _, _, _) =
Gooding.VLAMB(1.0, rlist[n1], rlist[n2], flist[n2] - flist[n1], tlist[n2] - tlist[n1]);
double vi = Sqrt(VR11 * VR11 + VT11 * VT11);
double vf = Sqrt(VR12 * VR12 + VT12 * VT12);
double diff1 = vlist[n1] - vi;
double diff2 = vlist[n2] - vf;
double diff = Sqrt(diff1 * diff1 + diff2 * diff2);
if (diff > diffmax)
{
diffmax = diff;
}
}
}
_testOutputHelper.WriteLine($"{diffmax}");
diffmax.ShouldBeZero(1e-10);
}
}
}
}