Possible bug when calculating the effect of AdditionalCelestialBody

[szabolcsvelkei]

Describe the bug JPL Spice kernels provide precisely calculated positions of celestial bodies, where their influence on each other is considered. If we want to calculate the movement of a spacecraft, we must do so either in the frame in which we are interested in the primary effects, or we must take into account not the full value of the effects of the additional celestial bodies, but the delta of the effect compared to the reference celestial body.

To Reproduce Steps to reproduce the behavior: In version 2.1.1, simply replace the Moon with a test satellite (use the moon's StateVector) and see how it behaves if you take into account the attraction of the earth and the sun. Example code:

var frame = Frame.ICRF;
var start = DateTimeExtension.J2000;
var end = DateTimeExtension.J2000.AddDays(3);
var sun = new CelestialBody(Stars.Sun);
var earth = new CelestialBody(PlanetsAndMoons.EARTH, frame, start);
var moon = new CelestialBody(PlanetsAndMoons.MOON, frame, start);
Mission mission = new Mission("mission01");
Scenario scenario = new Scenario("scn01", mission, new IO.Astrodynamics.Time.Window(start, end));
scenario.AddAdditionalCelestialBody(sun);
scenario.AddAdditionalCelestialBody(earth);
//Define test orbit
StateVector testOrbit = moon.GetEphemeris(start, earth, frame, Aberration.None).ToStateVector();
Clock clk = new Clock("My clock", 256);
Spacecraft spc = new Spacecraft(-1001, "MySpacecraft", 100.0, 10000.0, clk, testOrbit);
scenario.AddSpacecraft(spc);
var summary = await scenario.SimulateAsync(new DirectoryInfo("Simulation"), true, true, TimeSpan.FromSeconds(10.0));

var earthEphem = earth.GetEphemeris(new IO.Astrodynamics.Time.Window(start, end), earth, frame, Aberration.None, TimeSpan.FromMinutes(10)).ToArray();
var sunEphem = sun.GetEphemeris(new IO.Astrodynamics.Time.Window(start, end), earth, frame, Aberration.None, TimeSpan.FromMinutes(10)).ToArray();
var spcEphem = spc.GetEphemeris(new IO.Astrodynamics.Time.Window(start, end), earth, frame, Aberration.None, TimeSpan.FromMinutes(10)).ToArray();
var moonEphem = moon.GetEphemeris(new IO.Astrodynamics.Time.Window(start, end), earth, frame, Aberration.None, TimeSpan.FromMinutes(10)).ToArray();

Mat img = new Mat(1024, 1024, MatType.CV_8UC3);
img.SetTo(0);
double d = 1000000.0;
for (int n = 1; n < earthEphem.Length; n++)
{
    img.Line(new Point(img.Width / 2, img.Height / 2), new Point(img.Width / 2 + sunEphem[n - 1].ToStateVector().Position.X / d, img.Height / 2 + sunEphem[n - 1].ToStateVector().Position.Y / d), Scalar.Black, 1);
    // Earth at center
    img.Circle(new Point(img.Width / 2, img.Height / 2), (int)(earth.EquatorialRadius / d) + 1, Scalar.Blue, -1);
    // Moon path
    img.Circle(new Point(img.Width / 2 + moonEphem[n].ToStateVector().Position.X / d, img.Height / 2 + moonEphem[n].ToStateVector().Position.Y / d), (int)(moon.EquatorialRadius / d), Scalar.Gray, -1);
    // Calculated path
    img.Circle(new Point(img.Width / 2 + spcEphem[n].ToStateVector().Position.X / d, img.Height / 2 + spcEphem[n].ToStateVector().Position.Y / d), 1, Scalar.OrangeRed, -1);
    // Sun orientation
    img.Line(new Point(img.Width / 2, img.Height / 2), new Point(img.Width / 2 + sunEphem[n].ToStateVector().Position.X / d, img.Height / 2 + sunEphem[n].ToStateVector().Position.Y / d), Scalar.Yellow, 1);
    Cv2.ImShow("sim", img);
    Cv2.WaitKey(10);}
img.SaveImage("test.png");

Expected behavior We would expect the test satellite to follow (nearly) the same orbit as the moon nevertheless the attraction of the Sun dominates and immediately grabs the satellite from the Earth.

Screenshots

Desktop (please complete the following information):

• OS: Windows
• Version: 11

Additional context Regardless, it is very gratifying that such a solution is finally being prepared natively for .net. I hope you continue to work on the project with such intensity!

[sylvain-guillet]

Hello @szabolcsvelkei, I'am sorry for this inconvenience. I will analyze your problem and come back to you with an explanation. Thank you for your feedback and your encouragement.

[sylvain-guillet]

Hello @szabolcsvelkei, I've analyzed your scenario. To start my analysis I exported the scenario into Cosmographia to take a picture of initial state :

The framework proceeds a vector summation from Earth and sun gravity. At this initial position the influence of the sun is 3 times greater than the earth. So I put the sun as center motion and now I obtained an orbit close to the moon's orbit during 300 days, obviously orbit is not exactly the same due to mass difference, tidal effect, ... : .

This case works because the perturbing body is now the earth and the delta effect relative to the primary body(the sun) is lower. I think this test confirms your theory and I've to differentiate perturbing forces included into Spice kernels.

May you have resources to share with me about the computation of perturbation forces from Spice kernels ? It seems to be a simple vector subtraction from perturbing body to center of motion but I would like to explore this subject.

Thank you very much for your relevance.

[szabolcsvelkei]

Interesting, if I choose the Sun as the "observer", the simulation will be different but not good, again. So... I think that I found the source of the problem:

In IO.Astrodynamics/Propagator/Integrator.cs where the ComputeAcceleration() function computes the acceleration, the frame of the observer is used as reference. Usually, the frame of the observer is not an inertial frame. The correct computation uses an inertial frame such as the barycenter of the solar system.

So I made a little "hack" in the above function when using the previously mentioned test scenario with only the Sun and Earth as celestial bodies (obviously these two are the most significant). I still use the non-inertial frame but subtract the acceleration of Earth towards the Sun.

public Vector3 ComputeAcceleration(StateVector stateVector)
{
    Vector3 res = Vector3.Zero;
    foreach (var force in Forces)
    {
        if (force is GravitationalAcceleration)
        {
            if (((GravitationalAcceleration)force).CelestialBody.Name == "SUN")
            {

                var earth = new CelestialBody(PlanetsAndMoons.EARTH, Frame.ICRF, DateTimeExtension.J2000);
                var earthEphem = earth.GetEphemeris(stateVector.Epoch, earth, stateVector.Frame, Aberration.None);
                res -= force.Apply(earthEphem.ToStateVector());
            }
        }
        res += force.Apply(stateVector);
    }

    return res;
}

and the result was what I expected, the satellite trajectory was finally nearly the same as the Moon's trajectory.

Of course, there are still little differences because there are other celestial bodies in the solar system that we have to take into account. At least the Jupiter and the Saturn.

I am not in a position to dictate the development, but I think when using the SPICE data, the correct solution is to use an inertial frame for the calculation of the forces and take into account the perturbing effects of as many celestial bodies as possible since they were also taken into account during the calculation of de440s.bsp

[sylvain-guillet]

I understood the idea and I'am going to implement a more generic code to handle correctly perturbing bodies. I'll try to provide a new release this week.

[sylvain-guillet]

@szabolcsvelkei a new version(2.2.0) has been released. Now after 3 days of integration, the relative speed to the moon is lower than 10m/s with the sun as perturbing body and the earth as center of motion.

[szabolcsvelkei]

Hi @sylvain-guillet, I apologize for guiding you through the use of delta accelerations. Unfortunately, it is not sufficient, because it should also be considered that the arbitrary choice of observer gives a noninertial accelerating coordinate system. (10 m / sec is still exceptionally large after 3 days of integration.)

If you use barycenter as observer and revert the changes in EvaluateGravitationalAcceleration(...) just use Aberration.None like this:

public Vector3 EvaluateGravitationalAcceleration(OrbitalParameters.OrbitalParameters orbitalParameters)
{
    var sv = orbitalParameters.Observer as CelestialBody != this ? orbitalParameters.RelativeTo(this, Aberration.None).ToStateVector() : orbitalParameters.ToStateVector();

    return GravitationalField.ComputeGravitationalAcceleration(sv);
}

you'll get a deltaV of 7 cm/s (!) after 3 days of integration.

Even further, if you add Jupiter and Saturn to the integration, the error will be much lower after 3 days:

• 3mm /s ! (Of course, I used de440.bsp and the actual date because of the availability of the required spice kernels for this scenario)

So, I think you are very close to a precise solution: just have to use the center of mass of the solar system for the calculations, regardless of the user's chosen observer, and then transform the result into the desired coordinate system.

[szabolcsvelkei]

My annoying philosophical level problem is why it works so accurately without aberration and very poorly with LT? It's ruining me.

[sylvain-guillet]

I apologize for guiding you through the use of delta accelerations...

Hello @szabolcsvelkei, do not be sorry, you've identified an issue and suggest a solution it's all in your honor. It's my duty to offer a robust technical solution, it's my fault. I haven't done enough tests with scenarios where disturbing bodies have non negligible effect. But now thanks to you I'll be able to correct this point.

So I've adapted the propagator on my computer to use the sun as center of motion and I reached a differential velocity about 2cm/s after three days with earth and sun only.

My annoying philosophical level problem is why it works so accurately without aberration and very poorly with LT? It's ruining me.

My suggestion is now the computation is made in ICRF frame with sun as center and time relative to barycentric dynamic time. In other words with aberration set as "None" we have a geometric state of celestial bodies at barycentric epoch in the same referential. Now that we use barycentric+TDB approach bodies state already include gravitational delay and the LT will over correct bodies state. LT correction must be used only for observation purposes.

I will create some automated tests to check if there are no other failures with perturbing bodies and I will produce another version.

[sylvain-guillet]

Hello @szabolcsvelkei , I'am fixing the framework to improve the accuracy. I tried to reproduce the scenario when you reach a differential velocity about 3mm/s but I'am stuck at 6cm/s... To be sure there is no ambiguity in scenario could you please share with me the code of your scenario ? About data, can you give me the list of kernels you used with your scenario ? Thanks for your help !

[Edit] Ok I understood what you have done. I tried to reproduce your scenario with the code of 2.1.1 but I think you've adapted the code a little to support Barycenter as center of motion (else the solar system barycenter mass is used in the computation and it's wrong, we want only the relative position). However use barycenter as perturbing object can be a good idea for distant system and this is planned for next version... So i'am going to prioritize this feature, the first tests give me an accuracy of 2.8mm/s.

[sylvain-guillet]

Hi @szabolcsvelkei, The bug has been fixed in the new release (3.0.0), now after three days with a step of 100s I reached a difference in position about 35m and 0.27mm/s for the velocity. This scenario should be simulated in less than 300ms. Here the code for your project adapted to the V 3.0.0 :

var frame = Frames.Frame.ICRF;
var start = DateTimeExtension.J2000;
var end = start.AddDays(3);
var earth = new CelestialBody(PlanetsAndMoons.EARTH, frame, start);
var moon = new CelestialBody(PlanetsAndMoons.MOON, frame, start);
Astrodynamics.Mission.Mission mission = new Astrodynamics.Mission.Mission("mission01");
Scenario scenario = new Scenario("scn01", mission, new Window(start, end));

StateVector testOrbit = moon.GetEphemeris(start, new CelestialBody(399), frame, Aberration.None).ToStateVector();
Clock clk = new Clock("My clock", 256);
Spacecraft spc = new Spacecraft(-1001, "MySpacecraft", 100.0, 10000.0, clk, testOrbit,0.5);
scenario.AddSpacecraft(spc);
scenario.AddCelestialItem(new CelestialBody(10));
scenario.AddCelestialItem(new CelestialBody(399));
scenario.AddCelestialItem(new Barycenter(1));
scenario.AddCelestialItem(new Barycenter(2));
scenario.AddCelestialItem(new Barycenter(4));
scenario.AddCelestialItem(new Barycenter(5));
scenario.AddCelestialItem(new Barycenter(6));
scenario.AddCelestialItem(new Barycenter(7));
scenario.AddCelestialItem(new Barycenter(8));
var summary = await scenario.SimulateAsync(new DirectoryInfo("Simulation"), false, false, TimeSpan.FromSeconds(100.0));

var spcSV = spc.GetEphemeris(end, earth, Frames.Frame.ICRF, Aberration.None).ToStateVector();
var moonSV = moon.GetEphemeris(end, earth, Frames.Frame.ICRF, Aberration.None).ToStateVector();

var delta = spcSV - moonSV;
var deltaP = delta.Position.Magnitude();
var deltaV = delta.Velocity.Magnitude();

Assert.True(deltaP < 35);
Assert.True(deltaV < 2.7E-04);

If it's ok for you I will close the issue ? Thank you again for your contribution. Sylvain

[szabolcsvelkei]

Thanks @sylvain-guillet ,I think that this is a perfect solution of the problem! Sorry for the late reply. The issue can be closed.

Note: And yes I modified the previous version a little to bypass barycenters.