Ball script circular orbit setup

[seanmcleod]

@jonsberndt while debugging the issue with regards to the argument of perigee I noticed that the eccentricity for the ball script wasn't 0, it was 0.0015. The code use 1e-8 as a cut-off in determining whether the eccentricity is close enough to 0 or not.

I sort of wondered/assumed you had tried to set up a perfect circular equatorial orbit, but wasn't 100% certain if that was your intention and/or whether there was some floating point lack of precision playing a part. I didn't pursue it too much at the time.

However taking another look now I noticed the following when looking more closely at the ball script etc.

In ball_orbit.xml you specifically set the gravity model to 1, which is the WGS84 ellipse model.

jsbsim/scripts/ball_orbit.xml

Line 23 in 5c31c25

<property value="1"> simulation/gravity-model </property>

Which is the default gravity model.

jsbsim/src/models/FGInertial.cpp

Lines 56 to 62 in 5c31c25

	// Earth defaults
	double RotationRate = 0.00007292115;
	GM = 14.0764417572E15; // WGS84 value
	J2 = 1.08262982E-03; // WGS84 value for J2
	a = 20925646.32546; // WGS84 semimajor axis length in feet
	b = 20855486.5951; // WGS84 semiminor axis length in feet
	gravType = gtWGS84;

Then in reset00_v2.xml there are the following values and comments with regards to setting the initial velocity.

jsbsim/aircraft/ball/reset00_v2.xml

Lines 16 to 19 in 5c31c25

	<velocity unit="FT/SEC" frame="BODY">
	<!--x> 23869.978340 </x--> <!-- For use with standard grav model -->
	<x> 23889.145167 </x> <!-- For use with WGS84 grav model -->
	</velocity>

If I simply switch the x velocity to use the value calculated supposedly for the spherical earth gravity model then the eccentricity drops down to effectively 0.

Eccentricity - when using 23889.145167 supposedly WGS84 gravity model
0.0015065501444178064

Eccentricity - when using 23869.978340 supposedly standard (spherical earth) gravity model
5.2821191864893535e-11

So are the two x velocity values and their comments simply switched around, mixed up in reset00_v2.xml?

@bcoconni I noticed that you tweaked the values in the last commit for this file in 2016, although you were simply tweaking the values after the decimal point - 36de9ac

[jonsberndt]

Hi, Sean, Yes, it’s very possible that things got switched around. I’ll try and take a look at it soon, but can’t get to it today. Jon

[seanmcleod]

So I had a quick look at calculating the orbital velocity.

https://en.wikipedia.org/wiki/Vis-viva_equation

$$ v^2 = GM \left(\frac{2}{r} - \frac{1}{a}\right) $$

Given a circular orbit a = r.

jsbsim/src/models/FGInertial.cpp

Lines 58 to 62 in 5c31c25

	GM = 14.0764417572E15; // WGS84 value
	J2 = 1.08262982E-03; // WGS84 value for J2
	a = 20925646.32546; // WGS84 semimajor axis length in feet
	b = 20855486.5951; // WGS84 semiminor axis length in feet
	gravType = gtWGS84;

Using the Vis-viva equation:

>>> import math
>>>
>>> alt = 800000
>>> GM = 14.0764417572e15
>>> earth_a = 20925646.32546
>>>
>>> r = earth_a + alt
>>> r
21725646.32546
>>>
>>> v = math.sqrt(GM*(2/r - 1/r))
>>>
>>> v
25454.237455218088
>>>

Which is a much higher velocity. So I must be missing something in terms of calculating the required orbital velocity.

I double-checked my calculation for r, by checking vInertialPosition in the debugger.

jsbsim/src/models/FGPropagate.cpp

Lines 499 to 502 in 5c31c25

	void FGPropagate::ComputeOrbitalParameters(void)
	{
	const FGColumnVector3 Z{0., 0., 1.};
	FGColumnVector3 R = VState.vInertialPosition;

Which matches my calculation of r above.

VState.vInertialPosition
{data=0x00000000006f9c20 {21725646.325459998, 0.0000000000000000, 0.0000000000000000} }

[bcoconni]

The thing is that the velocity provided in ball/reset00_v2.xml is in the ECEF frame while the vis-viva computes the formula in the ECI frame. So you must transform the velocity from ECI to ECEF:
$$v_e=v_i-r\omega_{Earth}$$

jsbsim/src/models/FGInertial.cpp

Lines 56 to 57 in 5c31c25

	// Earth defaults
	double RotationRate = 0.00007292115;

Once this is done, you get the value from the ball script initial conditions

>>> omega = 0.00007292115
>>> v - r*omega
23869.97834067227

[seanmcleod]

Once this is done, you get the value from the ball script initial conditions

Ah, thanks @bcoconni. Although there is still the mix-up in the initial conditions file between which one is the correct one to use for WGS84.

Specifying a couple more decimal places for the velocity decreases the eccentricity even closer to 0.

v = 23869.978340          Eccentricity = 5.2821191864893535e-11
v = 23869.97834067227     Eccentricity = 6.6613381477509392e-16

I'll update the initial conditions file.