The ECI reference frame and tidal forces

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
Code:
Introduction
In the thread in the Maths & Physics section J4 - J5 perturbations, it was pointed out that that the gravitational perturbation calculations needed to be carried out in the Earth-centred Earth-fixed (ECEF) reference - i.e., one that is centred on the Earth but rotates along with the Earth.

However, to calculate the trajectory of a satellite in orbit reasonably close to the Earth, it is conventional to integrate the equations of motion in the Earth-centred inertial (ECI) reference frame. Although this reference frame includes a reference to 'inertial' in its name, it isn't really an inertial reference frame and if we want to do things properly we should really add tidal accelerations (due to the Sun and the Moon) in order to calculate the right trajectory.

In addition to dealing with tidal forces in the ECI reference frame, somewhere in the calculation of gravitational forces, one also needs to convert from the ECI reference frame to the ECEF reference frame (and back again) in order to add in the 'non-spherical' gravitational perturbations. as discussed in the aforementioned thread.

However, if we are going to do stuff in the ECI and ECEF reference frames, we first need to understand the ECI reference frame first.


What is the ECI refernce frame?
Orbiter bases all of its internal calculations on a 'global' reference frame notionally centred on the barycentre of the Solar System [see comment ^1 below]. The position and velocity of all bodies and vessels in Orbiter are tracked using this global reference frame. This reference to a base 'global' reference frame is conventional and JPL and all other relevant bodies use the same basic idea that all events in space-time are reference and related by some form of universal reference frame.

In Orbiter, this global reference frame is treated as a Newtonian inertial reference frame so that if one were able to remove all of the gravitational force from all of the planets and the Sun in the program, a spacecraft at rest in this reference would remain at rest [see comment ^2 below] thereby satisfying one of the standard Newtonian tests of an inertial reference frame.

However, what it also means is that if one wishes to describe the motion of a spacecraft near Earth, one is using a coordinate system referenced to the Solar System barycentre - a point in space roughly 1.5 x 10^11 m away. Although there is no reason why one couldn't carry out all of one's calculations in this global and inertial reference frame, it clearly isn't a natural representation of a satellite in orbit around the Earth where it would be far more natural to represent the equations of motion of the satellite in coordinates centred on the Earth. And if the Earth were an isolated (large) mass, then because of the huge difference in mass between the Earth and the satellite, a reference frame centred on the Earth would also be an inertial reference frame. In which case we would simply be left with the standard Keplerian equations of motion for gravity without further ado.

The problem is, though, that the Earth is not an isolated mass. The Earth and Moon are in orbit around the Earth-Moon barycentre and, in turn, the Earth-Moon barycentre is in orbit around the Sun. (To a very good approximation, we could ignore the gravitational contributions from the other planets in the Solar System but if we wanted to be pedantic about the thing, we may also want to take into account the gravitational contribution of all of the other planets as well.)

These more complex motion of the Earth created by the gravitational influence of other planetary bodies means that a non-rotating reference frame centred on the Earth is an accelerated reference frame because the Earth is in 'free fall' in its local gravitational field induced by the Sun, Moon and other bodies. Now if you could free fall in this accelerating frame of reference exaclty alongside the centre of the Earth, you would not notice this acceleration because you and the Earth would be in free fall together. In practice, though, you must be at some non-zero distance from the Earth (at least 6371 km, presumably) and so you cannot be in free fall exactly alongside the Earth. And this means that you experience a slightly different free fall acceleration from the Earth, which means that the acceleration of the Earth appears to be slightly different from your own. The closer you are to the centre of the Earth, the smaller this effect. This anomalous acceleration differential between two closeby bodies in free fall is usually termed a 'tidal force' because from an observer's standpoint on the surface of the Earth, there appears to be an anomalous tidal force tugging on things (witness the behaviour of the oceans).

The ECI reference frame is, then, a non-rotating reference frame centred on the Earth that takes account of the accelerated motion of the Earth by adding tidal contributions due to the differential gravitational influence due to other massive bodies (e.g., the Sun and the Moon) into the equations of motion in the ECI reference frame. These tidal contributions are 'fictitious' forces (in much the same way that the centrifugal force in a rotating reference frame is a 'fictitious' force) but everything works out so long as we make sure that we include these tidal forces in our ECI reference frame calculations.

OK, so how do we calculate the tidal forces so that we can include them in our calculations?


Derivation of the ECI equations of motion
To begin with, let's simplify things and ignore the gravitational contribution from the Sun and just focus on the Earth-Moon system. In this case, if we fix the origin of our coordinate system on the Earth-Moon barycentre, we should have an inertial reference frame. In this case, we can write down our equations of motion as:

[MATH]x''(t) = -\frac{\mu _e}{\Delta r_e^3}\left(x(t)-x_e(t)\right)-\frac{\mu _m}{\Delta r_m^3}\left(x(t)-x_m(t)\right)[/MATH]
[MATH]y''(t) = -\frac{\mu _e}{\Delta r_e^3}\left(y(t)-y_e(t)\right)-\frac{\mu _m}{\Delta r_m^3}\left(y(t)-y_m(t)\right)[/MATH]
[MATH]z''(t) = -\frac{\mu _e}{\Delta r_e^3}\left(z(t)-z_e(t)\right)-\frac{\mu _m}{\Delta r_m^3}\left(z(t)-z_m(t)\right)[/MATH]
where:

* [MATH]\left(x(t), y(t), z(t)\right)[/MATH] is the location of the satellite in the EMB inertial reference frame;

* [MATH]\left(x_e(t), y_e(t), z_e(t)\right)[/MATH] is the location of the Earth in the EMB inertial reference frame;

* [MATH]\left(x_m(t), y_m(t), z_m(t)\right)[/MATH] is the location of the Moon in the EMB inertial reference frame;

* [MATH]\Delta r_e^3[/MATH] is the standard Euclidian distance of the satellite from the centre of the Earth;

* [MATH]\Delta r_m^3[/MATH] is the standard Euclidian distance of the satellite from the centre of the Moon;

* [MATH]\mu_e[/MATH] is the gravitational parameter ('GM') for the Earth; and

* [MATH]\mu_m[/MATH] is the gravitational parameter ('GM') for the Moon.

What you see here in these equations is on the right-hand side of each a force due to the gravitational attraction of the Earth; and force due to the gravitation attraction of the Moon. Although it looks a little complex, each is just the component of the normal Keplerian contribution from each body as one might expect.

So, far we've been working in the inertial reference frame centred on the EMB. Let's shift the origin of our coordinate system to the centre of the Earth by introducing the new variables:

[MATH]X(t) = x(t) - x_e(t)[/MATH][MATH]Y(t) = y(t) - y_e(t)[/MATH][MATH]Z(t) = z(t) - y_e(t)[/MATH]
[MATH]X_m(t) = x_m(t) - x_e(t)[/MATH][MATH]Y_m(t) = y_m(t) - y_e(t)[/MATH][MATH]Z_m(t) = z_m(t) - y_e(t)[/MATH]
These new spatial coordinates are just a re-expression of the position of the satellite and the Moon relative to the centre of the Earth rather than the EMB. If we make these substitutions in the above equations of motion, we get:

[MATH]X''(t) + x_e''(t) = -\frac{\mu _e}{\Delta r_e^3}X(t)-\frac{\mu _m}{\Delta r_m^3}\left(X(t)-X_m(t)\right)[/MATH]
[MATH]Y''(t) + y_e''(t) = -\frac{\mu _e}{\Delta r_e^3}Y(t)-\frac{\mu _m}{\Delta r_m^3}\left(Y(t)-Y_m(t)\right)[/MATH]
[MATH]Z''(t) + z_e''(t) = -\frac{\mu _e}{\Delta r_e^3}Z(t)-\frac{\mu _m}{\Delta r_m^3}\left(Z(t)-Z_m(t)\right)[/MATH]
In this expression, the values [MATH]\Delta r_e[/MATH] and [MATH]\Delta r_m[/MATH] don't change although we now calculate them as:

[MATH]\Delta r_e = \sqrt{X(t)^2 + Y(t)^2 + Y(t)^2}[/MATH]
[MATH]\Delta r_m = \sqrt{\left(X(t)-X_m(t)\right)^2 + \left(Y(t)-Y_m(t)\right)^2 + \left(Z(t)-Z_m(t)\right)^2}[/MATH]
This just leaves the terms like [MATH]x_e''(t)[/MATH[. What is this? If we remember that the Earth is indeed accelerating due to gravitational from the Moon only in this simple Earth-Moon system, we can immediately write: [MATH]x_e''(t) = -\frac{\mu_m}{\Delta r_{e,m}^3}\left(x_e(t) - x_m(t)\right)[/MATH]
[MATH]y_e''(t) = -\frac{\mu_m}{\Delta r_{e,m}^3}\left(y_e(t) - y_m(t)\right)[/MATH]
[MATH]z_e''(t) = -\frac{\mu_m}{\Delta r_{e,m}^3}\left(z_e(t) - z_m(t)\right)[/MATH]
where [MATH]\Delta r_{e,m}[/MATH] is the normal Euclidean distance between the Earth and the Moon. The right-hand side of these equations can be thought of as having been expressed in terms of the initial EMB-centric coordinates. If we re-express this in terms of the Earth-centric coordinates, we get:

[MATH]x_e''(t) = +\frac{\mu_m}{\Delta r_{e,m}^3}X_m(t)[/MATH]
[MATH]y_e''(t) = +\frac{\mu_m}{\Delta r_{e,m}^3}Y_m(t)[/MATH]
[MATH]z_e''(t) = +\frac{\mu_m}{\Delta r_{e,m}^3}Z_m(t)[/MATH]
with

[MATH]\Delta r_{e,m} = \sqrt{X_m^2 + Y_m^2 + Z_m^2}[/MATH]
so that in the ECI reference frame, our equations of motion become:

[MATH]X''(t) = -\frac{\mu _e}{\Delta r_e^3}X(t) + \left(\frac{\mu _m}{\Delta r_m^3}\left(X_m(t)-X(t)\right) - \frac{\mu_m}{\Delta r_{e,m}^3}X_m(t)\right)[/MATH]
[MATH]Y''(t) = -\frac{\mu _e}{\Delta r_e^3}Y(t) + \left(\frac{\mu _m}{\Delta r_m^3}\left(Y_m(t)-Y(t)\right) - \frac{\mu_m}{\Delta r_{e,m}^3}Y_m(t)\right)[/MATH]
[MATH]Z''(t) = -\frac{\mu _e}{\Delta r_e^3}Z(t) + \left(\frac{\mu _m}{\Delta r_m^3}\left(Z_m(t)-Z(t)\right) - \frac{\mu_m}{\Delta r_{e,m}^3}Z_m(t)\right)[/MATH]
So, let's quickly recap. The above are the equations of motion for the satellite in the ECI reference frame. On the right-hand side in each equation, we have two terms. The first is the usual Keplerian terms due to the gravitational attraction of the Earth; and the second is the tidal force due to the gravitational influence of the Moon in this case. This tidal force is clearly the difference between the acceleration of the satellite due to the gravitational influence of the Moon and the acceleration of the Earth due to the gravitational influence of the Moon. It is a small differential that gets smaller the closer to the centre of the Earth we get.


A generalisation of the ECI equations of motion
OK, so far we have worked out the ECI reference frame equations of motion for the simple Earth-Moon system. What happens if we include the gravitational attraction of other bodies - most notably the Sun? It turns out that if we add other gravitating bodies into the mix, the equations of motion generalise to:

[MATH]X''(t) = -\frac{\mu _e}{\Delta r_e^3}X(t) + \sum_i\left(\frac{\mu _i}{\Delta r_i^3}\left(X_i(t)-X(t)\right) - \frac{\mu_i}{\Delta r_{e,i}^3}X_i(t)\right)[/MATH]
[MATH]Y''(t) = -\frac{\mu _e}{\Delta r_e^3}Y(t) + \sum_i\left(\frac{\mu _i}{\Delta r_i^3}\left(Y_i(t)-Y(t)\right) - \frac{\mu_i}{\Delta r_{e,i}^3}Y_i(t)\right)[/MATH]
[MATH]Z''(t) = -\frac{\mu _e}{\Delta r_e^3}Z(t) + \sum_i\left(\frac{\mu _i}{\Delta r_i^3}\left(Z_i(t)-Z(t)\right) - \frac{\mu_i}{\Delta r_{e,i}^3}Z_i(t)\right)[/MATH]
Here, the index [MATH]i[/MATH] runs over all of the bodies other than the Earth that exert a gravitational influence over the satellite (and the Earth). For the Earth, it is usually sufficient to include just the Sun and the Moon so that we can set [MATH]i=1[/MATH] for the Sun and [MATH]i=2[/MATH] for the Moon.


How do we calculate the tidal force terms?
So far, we have worked out the equations of motion of a satellite in the ECI reference frame. To calculate the acceleration of the satellite in this frame, we need to calculate the acceleration due to the Earth (the first term in the equations on the right-hand side); and the tidal force terms for each of the other gravitating bodies that we want to include in our gravitational model.

Note, though, that these are not equations of motion for the other gravitating bodies - e.g., the Moon and the Sun. They are only equations of motion for the satellite in which case we presume that we can calculate the position of these gravitating bodies relative to the Earth for any time, [MATH]t[/MATH]. In practice, what we do is use an ephemeris system (such as Orbiter's internal VSOP87 calculations) that allow us to calculate the position of the relevant gravitating bodies (e.g., the Sun and the Moon) in a global, inertial reference frame; and then subtract of the position of the position of the Earth to calculate the same positions relative to the Earth.


An example calculation
So that we can see how this works - and the magnitude of the various components - let's calculate the forces - Keplerian and tidal - acting on a satellite in Low Earth Orbit using the above framework. To be specific, let's use Orbiter's default "Docked at ISS" scenario for the Delta-Glider and use a simple Lua script to extract the relevant information about the position of the ISS, the Moon and the Sun relative to the Earth. From this, we'll calculate the Keplerian and tidal contributions to the ISS's acceleration in the ECI reference frame.

The Lua script is:
Code:
oapi.set_pause(true)

mjd     = oapi.get_simmjd()

iss     = vessel.get_handle("ISS")
earth   = oapi.get_objhandle("Earth")
moon    = oapi.get_objhandle("Moon")
sun     = oapi.get_objhandle("Sun")

-- get the current location of ISS
q_iss   = oapi.get_globalpos(iss)

-- get the current location of Earth
q_ear   = oapi.get_globalpos(earth)

-- get the current location of the Moon
q_mon   = oapi.get_globalpos(moon)

-- get the current location of Earth
q_sun   = oapi.get_globalpos(sun)

-- print out the results to a text file
file = io.open("testmb2012.txt", "w")
io.output(file)
file:write("q_v  ", q_iss.x, "  ", q_iss.y, "  ", q_iss.z, "\n")
print()
file:write("q_e  ", q_ear.x, "  ", q_ear.y, "  ", q_ear.z, "\n")
print()
file:write("q_m  ", q_mon.x, "  ", q_mon.y, "  ", q_mon.z, "\n")
print()
file:write("q_s  ", q_sun.x, "  ", q_sun.y, "  ", q_sun.z, "\n")
print()
file:close()

Since we haven't fixed the date at which these positions are sampled, if you were to run this script you will probably get slightly different values - but in my case the result of this sampling yields the following output:


Code:
q_v  -148373734893.55  26948501.361495  16086774535.065
q_e  -148375953122.08  20958013.704895  16088895065.916
q_m  -148619322792.77  47628076.771252  15794821861.364
q_s  -605098472.74882  21242703.919854  -772982266.62465


This gives the x-y-z coordinates of the ISS, the Earth, the Moon and the Sun in Orbiter's global inertial coordinate system. Values are in SI units metres. From this, we can calculate the coordinates of the ISS, the Moon and the Sun in the ECI reference frame simply by subtracting the Earth coordinates from the three others:

[MATH](X, Y, Z) = (2218228.53, 5990487.66, -2120530.85 )\,\text{m}[/MATH]
[MATH](X_m, Y_m, Z_m) = (2296995275914.85, 26670063.07, -294073204.55 )\,\text{m}[/MATH]
[MATH](X_s, Y_s, Z_s) = (147770854649.33, 284690.21, -16861877332.54 )\,\text{m}[/MATH]
To be sure, some of these are big numbers - but that's just a reflection of using SI units for these calculation. Now, if we use the following values for the gravitational parameters for the Earth, Moon and Sun as listed in Wikipedia.

[MATH]\mu_e = 3.986004418\times 10^{14}[/MATH][MATH]\mu_m = 4.9048695\times 10^{12}[/MATH][MATH]\mu_s = 1.32712440018 \times 10^{20}[/MATH]
then we can calculate the Keplerian gravitational acceleration due to the Earth (the first term on the right-hand side of the equations of motion in the ECI reference frame):

[MATH]\left(g_{e,x},g_{e,y},g_{e,z}\right) = \left(-2.89969,-7.83083,2.77198\right)[/MATH] ms^-2

The magnitude of the gravitational acceleration is around 8.80 ms^-2. Which is what one would expect for an object at an orbital altitude slightly less than 400 km.

What about the tidal contribution from the Moon? If we work our way through the calculations, we calculate that the tidal acceleration due to the Moon is:

[MATH]\left(t_{m,x},t_{m,y},t_{m,z}\right) = 10^{-6}\,\left(-0.29185, -0.516213, 0.069545\right) [/MATH] ms^-2

Note the factor of [MATH]10^{-6}[/MATH] in front of the tidal acceleration vector. This says that this is a pretty small acceleration. It isn't zero, but it certainly isn't going to have a huge impact on the satellite's trajectory. Unless we want to do high fidelity trajectory planning for objects in Low Earth Orbit, we can do what Orbiter tools such as TransX normally do in this case which is to ignore it.

And what about the tidal contribution from the Sun? Again, if we work our way though the calculations, we find that the tidal contribution due to the Sun is:

[MATH]\left(t_{s,x},t_{s,y},t_{s,z}\right) = 10^{-6}\,\left(0.204405, -0.241657, 0.052007\right) [/MATH] ms^-2

Again note the leading factor of [MATH]10^{-6}[/MATH]. In other words, the tidal acceleration from the Sun is also very small (and somewhat smaller that the Moon's - about 50% of the Moon's tidal contribution in fact [see comment ^3 below]) so it too can generally be ignored for objects in Low Earth Orbit unless you are wanting to do high fidelity trajectory calculations.


In summary
And that's about all there is to it. So long as one includes the tidal forces attributable to the effects of other gravitating bodies, the ECI reference frame is equivalent to working in a 'pure' inertial coordinate system. And generally speaking the numbers used in the calculations are a lot more 'human friendly than those that appear when using the global inertial reference frame.


[^1]: Orbiter inherits this reference frame the underlying VSOP87 ephemeris system used to calculate the future past, present and future position of the planets (or, more correctly for the Earth Moon system, say, the position of the Earth-Moon barycentre (EMB)).

[^2]: One can approximate this scenario by adjusting the mass of all of the objects in the configuration files to just 1 kg each. Now, the motion of the Sun, planets ad moons will not be correct because those objects will continue to move along rails as if they were massive bodies - but if you use the scenario editor to place a vessel at an arbitrary point in space with zero velocity, then it should remain at rest.

[^3]: If this were not the case, there would be no such thing as spring and neap tides.
 
Last edited:

cristiapi

New member
Joined
May 26, 2014
Messages
222
Reaction score
0
Points
0
Location
Ancona
What about the tidal contribution from the Moon? If we work our way through the calculations, we calculate that the tidal acceleration due to the Moon is:

[MATH]\left(t_{m,x},t_{m,y},t_{m,z}\right) = 10^{-6}\,\left(-0.29185, -0.516213, 0.069545\right) [/MATH] ms^-2

Note the factor of [MATH]10^{-6}[/MATH] in front of the tidal acceleration vector. This says that this is a pretty small acceleration. It isn't zero, but it certainly isn't going to have a huge impact on the satellite's trajectory. Unless we want to do high fidelity trajectory planning for objects in Low Earth Orbit, we can do what Orbiter tools such as TransX normally do in this case which is to ignore it.

And what about the tidal contribution from the Sun? Again, if we work our way though the calculations, we find that the tidal contribution due to the Sun is:

[MATH]\left(t_{s,x},t_{s,y},t_{s,z}\right) = 10^{-6}\,\left(0.204405, -0.241657, 0.052007\right) [/MATH] ms^-2

I probably missed some calculations (I also see two lines with the error: "[LaTeX Error: Syntax error]"), but while for the Earth I get 8.8 m/s2 (as you wrote), I get about 5.8 mm/s2 for the Sun and about 33 [MATH]\mu[/MATH]m/s2 for the Moon.
 

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
(I also see two lines with the error: "[LaTeX Error: Syntax error]")

Yep. Thanks. Now fixed.

(I generally write these notes in a LaTeX-friendly Markdown script before transcribing them to OF. Because of OF limitations, this entails a manual editing of each and every LaTeX line and sometimes my editing isn't quite perfect).


I get about 5.8 mm/s2 for the Sun and about 33 \mum/s2 for the Moon.

Here's my C code snippet that I used to check the posted results:

Code:
#include <stdio.h>
#include <math.h>

int main(int argc, const char * argv[]) {
    
    // define arrays to house the position vectors of the ISS, the Earth, the Moon and
    // the Sun in the global coordinate system
    double iss[3];
    double ear[3];
    double mon[3];
    double sun[3];

    // inirialise the arrays with values extracted from Orbiter using the Lua script
    iss[0] = -148373734893.55;    iss[1] = 26948501.361495; iss[2] = 16086774535.065;
    ear[0] = -148375953122.08;    ear[1] = 20958013.704895; ear[2] = 16088895065.916;
    mon[0] = -148619322792.77;    mon[1] = 47628076.771252; mon[2] = 15794821861.364;
    sun[0] =    -605098472.74882; sun[1] = 21242703.919854; sun[2] =  -772982266.62465;
    
    // define arays to house poition vectors of the ISS, the Moon and the Sun in the ECI reference frame
    double Q [3];
    double Qm[3];
    double Qs[3];
    
    // calcuate the position vectors in the ECI reference frame
    Q [0] = iss[0] - ear[0]; Q [1] = iss[1] - ear[1]; Q [2] = iss[2] - ear[2];
    Qm[0] = mon[0] - ear[0]; Qm[1] = mon[1] - ear[1]; Qm[2] = mon[2] - ear[2];
    Qs[0] = sun[0] - ear[0]; Qs[1] = sun[1] - ear[1]; Qs[2] = sun[2] - ear[2];
    
    // set the standard gravitational parameters for the Earth, the Moon and the Sun
    double mu_e = 3.986004418e14;
    double mu_m = 4.9048695e12;
    double mu_s = 1.32712440018e20;
    
    // calculate the distance of the iss from the Earth, the Moon and the Sun
    double r_e = sqrt(Q[0]*Q[0] + Q[1]*Q[1] + Q[2]*Q[2]);
    double r_m = sqrt((Q[0]-Qm[0])*(Q[0]-Qm[0]) + (Q[1]-Qm[1])*(Q[1]-Qm[1]) + (Q[2]-Qm[2])*(Q[2]-Qm[2]));
    double r_s = sqrt((Q[0]-Qs[0])*(Q[0]-Qs[0]) + (Q[1]-Qs[1])*(Q[1]-Qs[1]) + (Q[2]-Qs[2])*(Q[2]-Qs[2]));
    
    // calculate the distance of the Earth from the Moon and the Sun
    double r_em = sqrt(Qm[0]*Qm[0] + Qm[1]*Qm[1] + Qm[2]*Qm[2]);
    double r_es = sqrt(Qs[0]*Qs[0] + Qs[1]*Qs[1] + Qs[2]*Qs[2]);
    
    
    // calculate normal Keplerian gravitational acceleration due to the Earth
    double g[3];
    double k;
    
    k = - mu_e / r_e / r_e / r_e;
    
    g[0] = k * Q[0];
    g[1] = k * Q[1];
    g[2] = k * Q[2];
    
    
    // calculate the tidal acccleration due to the Moon
    double tm[3];
    double k1, k2;
    
    k1 = mu_m / r_m  / r_m  / r_m ;
    k2 = mu_m / r_em / r_em / r_em;

    tm[0] = 1.e6 * (k1 * (Qm[0] - Q[0]) - k2 * Qm[0]);
    tm[1] = 1.e6 * (k1 * (Qm[1] - Q[1]) - k2 * Qm[1]);
    tm[2] = 1.e6 * (k1 * (Qm[2] - Q[2]) - k2 * Qm[2]);
    
    
    // calculate the tidal acccleration due to the Sun
    double ts[3];
    
    k1 = mu_s / r_s  / r_s  / r_s ;
    k2 = mu_s / r_es / r_es / r_es;
    
    ts[0] = 1.e6 * (k1 * (Qs[0] - Q[0]) - k2 * Qs[0]);
    ts[1] = 1.e6 * (k1 * (Qs[1] - Q[1]) - k2 * Qs[1]);
    ts[2] = 1.e6 * (k1 * (Qs[2] - Q[2]) - k2 * Qs[2]);
    
    
    // print out the results
    printf("Earth gravity:       (%2.6f, %2.6f, %2.6f) ms^-2 \n", g[0], g[1], g[2]);
    printf("\n");
    printf("Moon tidal acc.:     (%2.6f, %2.6f, %2.6f) micro-ms^-2 \n", tm[0], tm[1], tm[2]);
    printf("\n");
    printf("Sun tidal acc.:      (%2.6f, %2.6f, %2.6f) micro-ms^-2 \n", ts[0], ts[1], ts[2]);
    printf("\n");
    
    return 0;
}

Running this code produces:

Code:
Earth gravity:       (-2.899691, -7.830827, 2.771980) ms^-2 

Moon tidal acc.:     (-0.291854, -0.516213, 0.069545) micro-ms^-2 

Sun tidal acc.:      ( 0.204405, -0.241658, 0.052008) micro-ms^-2 

Program ended with exit code: 0

I don't think I've made a mistake in the maths or the coding - but if I have, please let me know.
 

cristiapi

New member
Joined
May 26, 2014
Messages
222
Reaction score
0
Points
0
Location
Ancona
My mistake, I was simply calculating the acceleration exerted by the celestial body on the ISS; now I got the point.

Do I just need to add those accelerations in the simulation?
 

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
Do I just need to add those accelerations in the simulation?

Yes, basically, the sum of the three accelerations is just the right-hand side of the following equations of motion in the ECI reference frame:

[MATH]X''(t) = -\frac{\mu _e}{\Delta r_e^3}X(t) + \sum_i\left(\frac{\mu _i}{\Delta r_i^3}\left(X_i(t)-X(t)\right) - \frac{\mu_i}{\Delta r_{e,i}^3}X_i(t)\right)[/MATH]
[MATH]Y''(t) = -\frac{\mu _e}{\Delta r_e^3}Y(t) + \sum_i\left(\frac{\mu _i}{\Delta r_i^3}\left(Y_i(t)-Y(t)\right) - \frac{\mu_i}{\Delta r_{e,i}^3}Y_i(t)\right)[/MATH]
[MATH]Z''(t) = -\frac{\mu _e}{\Delta r_e^3}Z(t) + \sum_i\left(\frac{\mu _i}{\Delta r_i^3}\left(Z_i(t)-Z(t)\right) - \frac{\mu_i}{\Delta r_{e,i}^3}Z_i(t)\right)[/MATH]
The result is the left-hand side which is the acceleration of the satellite in the ECI reference frame.
 

cristiapi

New member
Joined
May 26, 2014
Messages
222
Reaction score
0
Points
0
Location
Ancona
While for the very low orbit satellites nothing changes (as probably expected), for GPS satellites the result is worse.

Please, could you confirm the following?
All vectors and y[] is the position in the global reference frame:

r_e = y[sat] - y[earth] (and not y[earth] - y[sat])
r_m = y[sat] - y[moon]
r_s = y[sat] - y[sun]
r_em= y[moon] - y[earth]
r_es= y[sun ] - y[earth]
 

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
Please, could you confirm the following?
All vectors and y[] is the position in the global reference frame:

Code:
r_e = y[sat] - y[earth] (and not y[earth] - y[sat])
r_m = y[sat] - y[moon]
r_s = y[sat] - y[sun]
r_em= y[moon] - y[earth]
r_es= y[sun ] - y[earth]

Starting from a clean sheet, I re-derived the ECI equations of motion in Mathematica using a slightly different approach (Lagrangian mechanics) and got the same result. So, unless I'm having a brain-freeze, I think the equations of motion are OK.

Also note that:

Code:
temp = y[sat] - y[earth];
temp = dot product (temp, temp)
re     = sqrt(temp)

etc., for all of the distances.

As a further check, I'll set up a numerical experiment in Orbiter over the weekend. Starting from a set of state vectors for a satellite in something resembling a geostationary orbit, I'll integrate the ECI equations of motion in Mathematic; and then I'll let Orbiter do the same. And we'll see if the theory matches the 'experimental' result.
 

cristiapi

New member
Joined
May 26, 2014
Messages
222
Reaction score
0
Points
0
Location
Ancona
I have no doubt that your equations of motion are OK, I just wanted to ask for the correctness of my formulas for r_xx, just to be sure that, for example, r_em= y[moon] - y[earth] is correct and it's not r_em= y[earth] - y[moon].
Well, I'll try to do some other checks.
 

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
The signs are correct, I think.

But note that r_xx is always greater than zero because one is taking the dot product of the vector and then the positive square root. So, the sign shouldn't matter.

I'll double check via numerical experiment over the weekend that I haven't slipped up on a sign anywhere.
 

cristiapi

New member
Joined
May 26, 2014
Messages
222
Reaction score
0
Points
0
Location
Ancona
As I wrote, they are all vectors and if I change the sign, the result is totally different (and a bit strange).
 

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
As I wrote, they are all vectors and if I change the sign, the result is totally different (and a bit strange).

Indeed, your expressions are vectors.

In vector form, the ECI equations of motion are:

[MATH]\mathbf{Q}''(t) = -\frac{\mu _e}{\Delta r_e^3}\mathbf{Q}(t) + \sum_i\left(\frac{\mu _i}{\Delta r_i^3}\left(\mathbf{Q}_i(t)-\mathbf{Q}(t)\right) - \frac{\mu_i}{\Delta r_{e,i}^3}\mathbf{Q}_i(t)\right)[/MATH]
where [MATH]\mathbf{Q}(t)[/MATH] is the position vector of the satellite in ECI coordinates such that:

[MATH]\mathbf{Q}(t) = \left(X(t), Y(t), Z(t)\right)[/MATH]
and [MATH]\mathbf{Q}_i(t)[/MATH is the position vector of body [MATH]i[/MATH] in ECI coordinates such that:

[MATH]\mathbf{Q}_i(t) = \left(X_i(t), Y_i(t), Z_i(t)\right)[/MATH]
For the satellite, the ECI coordinates are related to the global coordinates by:

[MATH]\mathbf{Q}(t) = \mathbf{q}(t) - \mathbf{q}_e(t)[/MATH]
and for body [MATH]i[/MATH] by:

[MATH]\mathbf{Q}_i(t) = \mathbf{q}_i(t) - \mathbf{q}_e(t)[/MATH]
where

[MATH]\mathbf{q}(t) = \left(x(t), y(t), z(t)\right)[/MATH]
and

[MATH]\mathbf{q}_i(t) = \left(x_i(t), y_i(t), z_i(t)\right)[/MATH]
Note, though, that [MATH]\Delta r_e[/MATH], [MATH]\Delta r_i[/MATH] and [MATH]\Delta r_{e,i}[/MATH] are all scalars.
 

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
I've had a chance to run a test of the tidal contributions in the equations of motion in the ECI reference frame.

A brief description of the test
Before presenting the results, I should describe the specific test performed.

* I placed a Delta-Glider in a circular orbit around the Earth with an orbital radius of 40,000 km and coplanar with the Earth's equator

* A Lua script was run to extract the initial ECI coordinates of the Delta-Glider. Orbiter 2016 was then run for around 6,000 seconds and the Lua script was run for a second time to extract the final ECI coordinates of the Delta-Glider

* Using the initial ECI coordinates, the equations of motion of the Delta-Glider was set up and integrated forwards in time using Mathematica for the same length of time as Orbiter 2016 integrated forwards.

* A comparison was then made of the separation of the two final ECI coordinates - i.e., the difference between the final position as calculated by Orbiter 2016; and the final position as calculated by Mathematica using the ECI reference frame equations of motion.

This test was run twice - once with the tidal terms of the gravitational force in the Mathematica model; and once with those terms removed.

(N.B. In Orbiter 2016, non-spherical gravitational forces and radiation pressure were switched off to permit a like-for-like comparison)


The results
The results are as follows:

* With the tidal terms included in the Mathematica model, after 6,000 seconds and integrating forward from the same state vector in the ECI reference frame, the separation of the position of the Delta-Glider as calculated by Mathematica and by Orbiter 2016 was 6.04 m. (During this time, the Delta-Glider moved roughly 20 million metres relative to the Earth.)

* With the tidal terms excluded in the Mathematica model, after 6,000 seconds and integrating forward from the same state vector in the ECI reference frame, the separation of the position of the Delta-Glider as calculated by Mathematica and by Orbiter 2016 was 62.33 m.

In other words, excluding the tidal terms from the ECI equations of motion, leads to an inferior estimate of the final position of the Delta-Glider. Hence, the tidal terms are material for orbits in near GSO.

Moreover, if one changes the sign of the tidal terms in the ECI equations of motion, the final error in position rises to 122 m - or roughly twice the error arising from excluding the tidal terms from the equations of motion. This suggests that the sign of the tidal terms in the ECI equations of motion as described above is correct.

---------- Post added at 06:15 AM ---------- Previous post was at 04:11 AM ----------

As a bit of an addendum to the above (and as an aide-memoire for me), the equations integrated in Mathematica were (in vector form - see above):

[MATH] \mathbf{Q}''(t) = -\mu _e\,\frac{\mathbf{Q}(t)}{\left|\mathbf{Q}(t)\right|^3} + \sum_i\mu _i\,\left(\frac{\mathbf{Q}_i(t)-\mathbf{Q}(t)}{\left|\mathbf{Q}_i(t)-\mathbf{Q}(t)\right|^3} - \frac{\mathbf{Q}_i(t)}{\left|\mathbf{Q}_i(t)\right|^3}\right) [/MATH]
Now, to integrate these equations of motion, one needs to be able to calculate [MATH]\mathbf{Q}_i(t)[/MATH] - i.e., the position vector of body [MATH]i[/MATH] at time [MATH]t[/MATH] n the ECI reference frame. In this case, the two relevant gravitating bodies are the Moon and the Sun. Rather than use Mathematica's somewhat dubious and slow ephemeris calculator, it is possible to augment the equations of motion of the satellite in the ECI reference frame with equations of motion for the Moon and the Sun in the ECI reference frame.

With a little effort, one can show that the equations of motion for the Moon and the Sun in the ECI reference frame are:

[MATH] \mathbf{Q}_m''(t) = -(\mu _e+\mu_m)\,\frac{\mathbf{Q}_m(t)}{\left|\mathbf{Q}_m(t)\right|^3} + \mu _s\,\left(\frac{\mathbf{Q}_s(t)-\mathbf{Q}_m(t)}{\left|\mathbf{Q}_s(t)-\mathbf{Q}_m(t)\right|^3} - \frac{\mathbf{Q}_s(t)}{\left|\mathbf{Q}_s(t)\right|^3}\right) [/MATH]
and

[MATH] \mathbf{Q}_s''(t) = -(\mu _e+\mu_s)\,\frac{\mathbf{Q}_s(t)}{\left|\mathbf{Q}_s(t)\right|^3} + \mu _m\,\left(\frac{\mathbf{Q}_m(t)-\mathbf{Q}_s(t)}{\left|\mathbf{Q}_m(t)-\mathbf{Q}_s(t)\right|^3} - \frac{\mathbf{Q}_m(t)}{\left|\mathbf{Q}_m(t)\right|^3}\right) [/MATH]
In these last two equations, the first term is the standard Keplerian term where the mass of the Earth is augmented by the mass of the Moon in the first case; and the mass of the Sun in the second.

The second term is a tidal term - in the first equation, the tidal force exerted by the Sun on the Moon; and in the second, the tidal force exerted by the Moon on the Sun. The latter is negligible, but for completeness the term was included in the equations of motion.

The Lua Script used to extract relevant information from Orbiter 2016 for the initial conditions is as follows:

Code:
oapi.set_pause(true)

mjd     = oapi.get_simmjd()

iss       = vessel.get_handle("GL-02")
earth		  = oapi.get_objhandle("Earth")
moon		  = oapi.get_objhandle("Moon")
sun       = oapi.get_objhandle("Sun")

-- get the current location of ISS
q_iss 	= oapi.get_globalpos(iss)
p_iss 	= oapi.get_globalvel(iss)
	
-- get the current location of Earth
q_ear 	= oapi.get_globalpos(earth)
p_ear   = oapi.get_globalvel(earth)
	
-- get the current location of the Moon
q_mon 	= oapi.get_globalpos(moon)
p_mon   = oapi.get_globalvel(moon)

-- get the current location of Earth
q_sun 	= oapi.get_globalpos(sun)
p_sun   = oapi.get_globalvel(sun)


file = io.open("testmb2012.txt", "w")
io.output(file)
file:write(mjd,"\n")
file:write("qv = \{", q_iss.x, ", ", q_iss.y, ", ", q_iss.z, "\}\n")
file:write("pv = \{", p_iss.x, ", ", p_iss.y, ", ", p_iss.z, "\}\n")
print()
file:write("qe = \{", q_ear.x, ", ", q_ear.y, ", ", q_ear.z, "\}\n")
file:write("pe = \{", p_ear.x, ", ", p_ear.y, ", ", p_ear.z, "\}\n")
print()
file:write("qm = \{", q_mon.x, ", ", q_mon.y, ", ", q_mon.z, "\}\n")
file:write("pm = \{", p_mon.x, ", ", p_mon.y, ", ", p_mon.z, "\}\n")
print()
file:write("qs = \{", q_sun.x, ", ", q_sun.y, ", ", q_sun.z, "\}\n")
file:write("ps = \{", p_sun.x, ", ", p_sun.y, ", ", p_sun.z, "\}\n")
print()
file:close()

oapi.set_pause(false)
 

cristiapi

New member
Joined
May 26, 2014
Messages
222
Reaction score
0
Points
0
Location
Ancona
Here's my experiment:
1) calculate the state of the GPS PRN 07 satellite with the SDP4 propagator for the TLE 18001.44476542
2) propagate the state with the DOPRI853 integrator
3) for all the subsequent TLEs, calculate the difference between the DOPRI position and the TLE position for the TLE epoch.

After 129.6 days of propagation (the last used TLE is 18131.07605472) I get the following errors:

no tidal accelerations: 4.4 km in the radius vector and 0 in the ground track (the satellite location on the Earth surface is exactly the same obtained with the TLE);

with tidal acceleration: 216.6 km in the radius vector and 52 km in the ground track; most of the error is along the ECI Z axis (212.2 km).
 

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
Code:
Here's my experiment:
1) calculate the state of the GPS PRN 07 satellite with the SDP4 propagator for the TLE 18001.44476542
2) propagate the state with the DOPRI853 integrator
3) for all the subsequent TLEs, calculate the difference between the DOPRI position and the TLE position for the TLE epoch.

I don't really know much about SDP4 but my brief reading of web-based literature suggests that they are built for speed and not necessarily high long-term accuracy. In other words, SDP4 may not be the best performance benchmark to use.
 
Last edited:

cristiapi

New member
Joined
May 26, 2014
Messages
222
Reaction score
0
Points
0
Location
Ancona
I certainly agree. The SGP4/SDP4 propagator should be used only for short term propagation (few days). For that reason, I use SDP4 just to calculate the state of the satellite for the TLE epoch (I don't use SDP4 for the propagation), while I use the DOPRI integrator to propagate the first TLE:
- obtain the initial state from the 1st TLE with SDP4 for the TLE epoch
- propagate the state to the last TLE with DOPRI
- calculate the distance from the position obtained with DOPRI and the position obtained with SDP4 from the last TLE for the TLE epoch.
 

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
I also don't know much about the DOPRI integrator either.

Most integrators are subject to errors of one kind or another. It's always a good idea to run an integrator against a known analytical solution (e.g., pure Keplerian elliptical motion, to see how accurate the integrator is under various conditions. If the integrator is accurate under controlled modelling assumptions where comparison with exact analytical results is possible, then one can add in perturbations of various kinds with confidence that the underlying integration algorithm is accurate and adding small perturbations shouldn't change the overall accuracy of the integration very much.

For example, you can set up initial conditions corresponding to a pure circular motion around a single gravitating body of about the size and mass of the Earth at a GSO-type radius and then integrate forward for a few hundred days. In this case, you can ignore J_n perturbations, atmospheric drag, and lunar/solar perturbations - i.e., pure Keplerian (circular) motion. How close to the analytical result is the integrated solution? Are we talking a cumulative position error of1m, or 1,000 km over the integration period using DOPRI?

(Or, to put it another way, DOPRI may itself be a source of integration error.)
 

cristiapi

New member
Joined
May 26, 2014
Messages
222
Reaction score
0
Points
0
Location
Ancona
Initial state (m, m/s):
POS= 26400e3, 0, 1000
VEL= 0, 3700, 50

GM= 398600.5e9 m3/s2 (WGS-84)

after 1000 days I get (the variable integration step is about 183 s):
POS= 19945923.2217823 16361506.2589119 221856.963325989
VEL= -2587.93081593407 2774.37916220418 37.3935826176335

the expected state is:
POS= 19945923.311264 16361506.1630025 221856.962033213
VEL= -2587.93079992651 2774.37917533393 37.393582795654

the errors are:
POS: 0.131176 m
VEL: 2.07042e-5 m/s

If I use the SPICE library conics_c() function to obtain the state after 1000 days, the result is slightly different and the position error becomes 0.132202 m.
 
Last edited:

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
OK, I guess that this means that we can rule out DOPRI as a significant source of error.

Let's focus on the TLEs.. My understanding is that you are using the TLEs from one epoch to inform the initial state vectors of the DOPRI integrator and then, after computing the final state vectors using the DOPRI integrate, comparing those with state vectors from the TLEs of a corresponding later epoch.

Again, I don't know much about the TLE tables, and I could easily be wrong, but the I suspect that the TLE's are largely empirical: using a series of ground-based observations, parameters of an integration model such as SGP4/SDP4 are adjusted to get the best possible match with the observations with respect to the trajectory propagation of the same integration model. This means a few things:

* the empirical orbital parameters calculated using this approach are not the actual orbital parameters of the satellite - there is both a measurement uncertainty and a model fitting uncertainty. Back of the envelope calculations suggest that you would only need a 0.01 m/s uncertainty in initial velocity components to give a future position uncertainty of a few hundred kilometres after integrating forward a few hundred days.

* the orbital parameters are a best fit to a specific integration model - probably SGP4/SDP4. If you then run the same 'best fit' orbital parameters through a different model such as DOPRI, it may yield to inferior forward predictions of the state vector because the orbital parameters were not fitted to observation using that model.

One of the advantages of Orbiter is that we can access the true state vectors of a satellite with the full precision of floating point arithmetic. Empirical error is all but irrelevant. And because we know all of the model parameters - body masses, positions and J_n components, etc., we can use Orbiter to check that any given integration model that we have built is well formulated.

Orbiter functions as a proxy for a real physical system but one from which observation error is essentially eliminated. We can use Orbiter to check that we can build or have built a well-formulated integration model. Finally, having done this, we can then gauge the significance of empirical measurement issues. But my best guess at the moment is that you are hitting up against orbital parameter measurement uncertainty problems inherent in the SGP4/SDP4 fitting of orbital parameters.
 

cristiapi

New member
Joined
May 26, 2014
Messages
222
Reaction score
0
Points
0
Location
Ancona
Are you saying that checking a simulation with another simulation is better than checking a simulation against reality? Surely TLEs are not perfect, but their uncertainty for the TLE epoch is far smaller than that of any other possible general purpose simulation and I wouldn’t define TLEs as “largely empirical” because they use exactly the same method used by anyone who want to integrate the state of a celestial body (JPL, IAA-RAS, IMCCE, me :)).
Before running my simulation I do a best fit with TLEs for the TLE epoch (I also use some additional constraints); if the fit is good (“good” depends on the satellite), I’m sure that the simulation is reasonably accurate.
If I use your tidal acceleration equations, the fit is always too bad; I still don’t know why, but I’m doing some additional check…
 
Top