How to get position and attitude of plane during simulation

[yuki-2000]

Hello.

I am currently creating a program on Google Colab to visualize the movement of airplanes as videos, because it is easier to intuitively understand the movement of the aircraft than to view the changes in each parameter over time on a graph.
I need position(x,y,z) variables and attitude variables(psi, theta, phi) in orthogonal coordinate system.
I try to search in the document and programs, and I found 'position/distance-from-start-lon-mt', 'position/distance-from-start-lat-mt', 'position/h-sl-meters' and 'attitude/phi-rad', 'attitude/theta-rad', 'attitude/psi-rad'
However, I have never learned how to express position and orientation taking into account the shape of the Earth, and I am confused.
I read the official documentation's Section 3: Formulation Manual, https://jsbsim.sourceforge.net/JSBSimReferenceManual.pdf.

My question is,

1. By default, in initial condition, where and in what direction is the aircraft facing in what coordinates ?
2. Can I use 'position/distance-from-start-lon-mt', 'position/distance-from-start-lat-mt', 'position/h-sl-meters' as x,y,z ?
3. How can I use rotation matrix from python? and which rotation matrix should I use to calculate the attitude in a fixed Cartesian coordinate system (probably North East Down, not considering the spherical shape of the Earth)?

To be more specific, this is the kind of program I want to create. (sorry for not English.)
https://colab.research.google.com/drive/1qv5Hll_v4UWcC6YVJDv2KYifhyz7vU5Q?usp=sharing
Videos are made using only python and you can play them in the google colab.

[seanmcleod70]

In terms of initializing the aircraft in terms of position and attitude here are the relevant initial condition properties you can use.

jsbsim/src/initialization/FGInitialCondition.h

Lines 181 to 187 in b005a91

	@property ic/theta-deg (read/write) Pitch angle initial condition in degrees
	@property ic/phi-deg (read/write) Roll angle initial condition in degrees
	@property ic/psi-true-deg (read/write) Heading angle initial condition in degrees
	@property ic/lat-gc-deg (read/write) Latitude initial condition in degrees
	@property ic/long-gc-deg (read/write) Longitude initial condition in degrees
	@property ic/h-sl-ft (read/write) Height above sea level initial condition in feet
	@property ic/h-agl-ft (read/write) Height above ground level initial condition in feet

jsbsim/src/initialization/FGInitialCondition.h

Lines 211 to 214 in b005a91

	@property ic/phi-rad (read/write) Roll angle initial condition in radians
	@property ic/psi-true-rad (read/write) Heading angle initial condition in radians
	@property ic/lat-gc-rad (read/write) Geocentric latitude initial condition in radians
	@property ic/long-gc-rad (read/write) Longitude initial condition in radians

jsbsim/src/initialization/FGInitialCondition.cpp

Lines 1586 to 1593 in b005a91

	PropertyManager->Tie("ic/lat-geod-rad", this,
	&FGInitialCondition::GetGeodLatitudeRadIC,
	&FGInitialCondition::SetGeodLatitudeRadIC);
	PropertyManager->Tie("ic/lat-geod-deg", this,
	&FGInitialCondition::GetGeodLatitudeDegIC,
	&FGInitialCondition::SetGeodLatitudeDegIC);
	PropertyManager->Tie("ic/geod-alt-ft", &position,
	&FGLocation::GetGeodAltitude);

2. Can I use 'position/distance-from-start-lon-mt', 'position/distance-from-start-lat-mt', 'position/h-sl-meters' as x,y,z ?

Yes, that should work.

3. How can I use rotation matrix from python? and which rotation matrix should I use to get the ?

You didn't finish your question. Do you need a rotation matrix, as opposed to simply making use of the attitude angles theta, phi, psi that JSBSim supplies on each time step?

Did you write the Colab Jupyter notebook you linked to? Looks useful.

[yuki-2000]

Thank you for the answer.

You didn't finish your question.

Sorry. I fixed the comment and the true sentence is "How can I use rotation matrix from python? and which rotation matrix should I use to calculate the attitude in a fixed Cartesian coordinate system (probably North East Down, not considering the spherical shape of the Earth)?"
This means that it is unclear which coordinate system these three Euler angles are based on. If it is the simplest NED (North, East, Down) frame, I can understand, but if it is the ECEF frame or ECI frame, I don't know how to interpret those Euler angles.
In 3D plot, I want to represent the aircraft's trajectory and attitude with arrows, and coordinate conversion is necessary to create those arrows.

Did you write the Colab Jupyter notebook you linked to? Looks useful.

Yes, I made this last year to study the airplane dynamics. I wonder if there is a video in the examples such and ladder kick .ipynb, or this types of lecture file of the coordinate system on which it is based.

[agodemar]

@yuki-2000 here is a project with Jupyter notebooks where a lot of work has been done to visualize trajectories in various simulation scenarios: WiReS/calc

This work was then refined and the results were used for this open access publication: Varriale C., De Marco A., et al., "Flight Load Assessment for Light Aircraft Landing Trajectories in Windy Atmosphere and Near Wind Farms", DOI: 10.3390/aerospace5020042.

Here is a concise tutorial on how to handle position: TUTORIAL 02

Here is an example that builds on previous tutorials: Waypoint based, automatic guidance 1

In this folder WiReS/calc/python you'll find the Python sources that implement the low-level functions.

Example of function: WiReS/calc/python/geography.py

[yuki-2000]

Thank you for sharing the codes. They loos very useful for me. It seems these are based on output xml files, so I will change them to more pythonic code like using fdm[] during simulation.
especially https://github.com/agodemar/WiReS/blob/master/calc/06_ManeuverTest_WaypointGuidance.ipynb is very interesting. I want to understand this in the future.

[seanmcleod70]

I've added some NEU (North East Up) from initial location properties to JSBSim which you could then use instead of position/distance-from-start-lon-mt or having to implement the equivalent function from @agodemar's TUTORIAL 02.

Will post a link to the proposed changes once I've done a bit more testing.

[yuki-2000]

Thank you very much. Then, I wait for you, studying https://github.com/agodemar/WiReS/blob/master/calc/02_Tutorial02_ComplexManeuver-UnitConversion-Trajectory.ipynb

[seanmcleod70]

This means that it is unclear which coordinate system these three Euler angles are based on. If it is the simplest NED (North, East, Down) frame, I can understand, but if it is the ECEF frame or ECI frame, I don't know how to interpret those Euler angles.

  PropertyManager->Tie("attitude/phi-rad", this, ePhi, &FGPropagate::GetEuler);
  PropertyManager->Tie("attitude/theta-rad", this, eTht, &FGPropagate::GetEuler);
  PropertyManager->Tie("attitude/psi-rad", this, ePsi, &FGPropagate::GetEuler);

  PropertyManager->Tie("attitude/phi-deg", this, ePhi, &FGPropagate::GetEulerDeg);
  PropertyManager->Tie("attitude/theta-deg", this, eTht, &FGPropagate::GetEulerDeg);
  PropertyManager->Tie("attitude/psi-deg", this, ePsi, &FGPropagate::GetEulerDeg);

And the documentation.

  /** Retrieves a vehicle Euler angle component.
      Retrieves an Euler angle (Phi, Theta, or Psi) from the quaternion that
      stores the vehicle orientation relative to the Local frame. The order of
      rotations used is Yaw-Pitch-Roll. The Euler angle with subscript (1) is
      Phi. Various convenience enumerators are defined in FGJSBBase. The
      relevant enumerators for the Euler angle returned by this call are,
      ePhi=1, eTht=2, ePsi=3 (e.g. GetEuler(eTht) returns Theta).
      units radians
      @return An Euler angle.
  */
  double GetEuler(int axis) const { return VState.qAttitudeLocal.GetEuler(axis); }

[yuki-2000]

Retrieves an Euler angle (Phi, Theta, or Psi) from the quaternion that stores the vehicle orientation relative to the Local frame. The order of rotations used is Yaw-Pitch-Roll.

This is the exactly the information what I look for !!
Thank you very much @seanmcleod70

[seanmcleod70]

This is the exactly the information what I look for !!

In general, you won't find everything detailed in and documented in the reference manual which is very old by now. So you often need to search through the source code to look at the comments and to look at the source code itself.

[seanmcleod70]

I've pushed a commit to a branch in my JSBSim fork showing what I'm thinking of submitting as a PR to the JSBSim repo.

seanmcleod70@b78ff97

Want to do a bit more testing first. Let me know if you have any comments.

In the meantime I'd suggest you continue using position/distance-from-start-lon-mt, position/distance-from-start-lat-mt, position/h-sl-meters. For small distances from the start there shouldn't be much difference between using those properties compared to the new NEU properties I'm proposing to expose.

[yuki-2000]

For small distances from the start there shouldn't be much difference between using those properties compared to the new NEU properties I'm proposing to expose.

Thank you.
Is this because the radius of the earth is very big and there is only a very small difference between arc length and chord length (or the strait line length along the NED frame direction), considering only the thousand of feet ? Is my understanding correct?

[seanmcleod70]

So in the NED/NEU case there is a tangential plane at the point chosen for the origin of the local frame.

In my proposed set of new properties the origin of the local plane is the starting/initial conditions (lat, lon) at 0 alt. The distance/vector then between the aircraft's current position and the origin is based on the subtracting the ECEF coordinates of the origin with the aircraft's current ECEF coordinates.

In terms of position/distance-from-start-lon-mt etc.

  /** Get the geodetic distance between the current location and a given
      location. This corresponds to the shortest distance between the two
      locations. Earth curvature is taken into account.
      @param target_longitude the target longitude in radians
      @param target_latitude the target geodetic latitude in radians
      @return The geodetic distance in feet between the two locations */
  double GetDistanceTo(double target_longitude, double target_latitude) const;

So yes, taking the earth's ellipsoid into account, and being the shortest distance it also means things like great circle distance which will show up for points quite far apart.

I'll reply to you later in terms of the video options.

[yuki-2000]

This is an example of video on Google Colab that I had in mind. (Output of the last cell)
The position and heading direction is animated.
https://drive.google.com/file/d/1ZAPO9KAV50Sssa9DDurwwt5qaoRbPTj4/view?usp=sharing

Aside from the dirtiness of my code, is there a demand for videos like this in example?

@agodemar I used your codes and that was very helpful. Thank you.

[yuki-2000]

The roll is not animated with this arrow way though.

[agodemar]

@yuki-2000 the video is not visible

[seanmcleod70]

@agodemar for some reason the link he shared was for the .ipynb file on Google Drive and not a direct Colab link. So you see the plain text contents of the .ipynb file, you have to then click on the "Open with Google Collaborat..." link at the top.

Or alternatively here is a colab link.

https://colab.research.google.com/drive/1ZAPO9KAV50Sssa9DDurwwt5qaoRbPTj4

Once all the cells have then run the animated video is in the last cell.

[yuki-2000]

@seanmcleod70 Thank you.
Sometimes, this happens.

[yuki-2000]

https://colab.research.google.com/drive/1mhslHUIR-Vtg1aPb0uQyfRbLl5CkkYk6?usp=sharing

@agodemar I made the notebook as a colab notebook. Please check it again.