TestInitialConditions.py test fails on Windows about jsbsim HOT 13 CLOSED

Refreshing my basic/limited quaternions knowledge I guess I can see that the results I saw with the C++ debugger are a valid set of Euler angles (90, 90, 90) for the quaternion, but I don't see how the results from the Python test case (0, 90, 90) are.

from jsbsim.

According to Euler angle theory, a value of theta = 90 deg is not allowed in the 3-2-1 convention. It's a singularity case, when the fuselage is aligned with the local vertical and nose pointing upwards.

As you see, for such a given attitude, with a given position of the wings w.r.t. Earth axes, if you assign 90 deg to theta, there are infinite possible combinations of (phi, psi). So, there are infinite triplets of Euler angles with an elevation angle theta relaxed to a value of 90 degrees.
On the contrary, this attitude is well represented by a unique (Euler-Rodriguez normalized) quaternion value.

from jsbsim.

@seanmcleod Let's keep the Python thing apart first. Do you get the same results if you execute the script directly with the JSBSim executable ? On my Linux box (fedora 27), when I run the script with the following command

$ JSBSim scripts/mk82_script.xml  --logdirectivefile=data_output/position.xml

then I get the following data in position.csv

Do you get the same result from the C++ executable ?

In my experience, the results of some computations is sensitive to math optimizations made by the compiler. Do you get the same result if you disable all the optimizations from the compiler ?

from jsbsim.

@agodemar nice diagram! I see Jon mentions briefly some of the issues with a pitch attitude of 90 degrees in the JSBSim manual in terms of the rocket launch example in terms of rocket guidance, and avoiding some of the issues etc.

@bcoconni I was wondering if some of the differences between the Windows build and the Linux build could be due to some slight differences in terms of some of the trigonometric functions like atan etc. But as you say it may also depend on compiler optimisations.

In my Python testing I was using a release version of JSBSim, since I previously had issues building a debug version of JSBSim and getting the debug version to build/link with the Python libs.

For my reports in the original post above regarding the C++ output I was seeing I was using a debug build of JSBSim running the mk82_script directly (no Python test case code involved) and reporting based on inspecting/executing code in the C++ debugger.

So I was using different compiler optimisations between the two. Also I wasn't looking at output logged, rather the variables in the debugger while executing Load_v1().

What I'll do is log the output as per your example and try both a release and debug version and report back with the results.

from jsbsim.

Some further thoughts: I just reviewed the code that converts quaternions to Euler angles and I think that the code at stake is

src/math/FGMatrix33.cpp lines 166-169

  if (data[8] == 0.0)
    mEulerAngles(1) = 0.5*M_PI;
  else
    mEulerAngles(1) = atan2(data[7], data[8]);

I guess that the initial conditions in the script mk82_script.xml corresponds to a case where both data[7] and data[8] are very close to 0.0

1. If data[8]==0.0 then the Euler angle phi is set to π/2 (this is what you get on Windows ?)
2. If the numerical round off is such that data[8] is slightly different from 0.0 then the second case of the condition applies and set phi to 0.0 because data[7]==0.0 (this is what I get on Linux ?)

from jsbsim.

So I ran both the Windows x64 release and debug builds.

.\Release\JSBSim.exe .\scripts\mk82_script.xml --logdirectivefile=data_output\position.xml
.\Debug\JSBSim.exe .\scripts\mk82_script.xml --logdirectivefile=data_output\position.xml

The data in position.csv is identical between the release and debug builds, so no differences in terms of compiler optimisations between the two Windows builds.

So what is interesting is the difference I see in the data logged to the .csv file and the data I see when I set a breakpoint in FGInitialCondition::Load_v1() and then execute orientation.GetEulerDeg(ePhi).

from jsbsim.

As @agodemar mentioned we are here in a "gimbal lock" where there is not a unique set of Euler angles to describe the position.

I am pretty sure that we are hitting a round off problem. The result of the formula I have shown above is summed up in the matrix below for a number of scenarios where data[7] and data[8] are either equal or very close to zero.

It illustrates that for very small variations of data[8] the angle phi has very different values. If it is confirmed that it is indeed the cause of the problem then I am unsure about the way to proceed:

1. Rewrite the Euler angles computations so that they are not so sensitive to round off (this is theoretically impossible but we might be able to alleviate the problem in most cases)
2. Force the Euler angles to special/arbitrary values for the gimbal lock so that we can predict the output for the script mk82_script
3. Modify TestInitialConditions.py to not check phi (or psi) when theta is equal to 90 deg.

from jsbsim.

At the time when I initially came across this issue I tested a couple of values for theta in reset00.xml to see how sensitive the test was to values around 90 degrees. For example a value of 89.9 for theta will result in the test passing.

I'd vote for option 3 if it is confirmed this issue is related to the round off example you've presented.

from jsbsim.

I ran the debugger and took a look at the values for data[7] and data[8] at the time that the output code requested the Euler angles to output them to the .csv file. And sure enough as you suspected:

Which ties up with the value of 180 degrees for phi that I see in my .csv file.

from jsbsim.

@bcoconni, @agodemar should I go ahead and submit a commit with option 3 above implemented?

from jsbsim.

@seanmcleod Sure, go ahead.

from jsbsim.

@bcoconni I submitted a pull request for this - #43

from jsbsim.

PR #43 has been merged (commit 470fdce)

from jsbsim.