Orbital Aero Model: Orientations & Rotations

Orientations in simulations are commonly described via Euler angles or Quaternions.

Euler Angles

Euler angles are a succession of 3 rotations about the axes of a coordinate system. You can read the three Euler angles and visualize in your mind how the object is orientated in space. The order of the rotations matters, but there is no universal order for the rotations. The downside with Euler Angles is they suffer from gimbal lock - at pitch angles of +90 or -90 the rotations fall apart and special conditions have to be made.

https://en.wikipedia.org/wiki/Euler_angles

Quaternions

Quaternions are far less intuitive than Euler angles. A Quaternion representation of an object's orientation consists of 4 numbers: w, x, y, and z (or a, b, c, d, or s, x, y, z or a scalar & a vector or ...). However, even if you know those numbers it's (near?) impossible to visualize in your mind how the object is orientated in space. Rotations of Quaternions are more efficient than Euler angles for a computer to calculate, and they don't have the gimbal lock problem.

https://en.wikipedia.org/wiki/Quaternion

Angles in Orbital Aero Model

For any project doing significant work with angles, you definitely want to use Quaternions internally. However when communicating angles to the user, it's better to use Euler angles.

Using Quaternions is actually fairly simple, but getting a solution off the Internet that "just works" is more difficult than it should be. The options were either gigantic libraries that I only needed a tiny piece of, smaller convoluted source code snippets, or straightforward code snippets that actually failed under certain circumstances!

My solution was to take the best pieces of several sources and simplify them down to just what I needed. A way to represent an orientation, a way to rotate the object in space, and a way to rotate a vector (such as thrust) by an orientation.

Body Axis

A body's axis is defined as follows:
x forward
y right
z down

For my Euler angles, I'm using the convention of rotation of z followed by y followed by x. This matches the IEEE 1278.1 DIS (Distributed Interactive Simulation) standard used by military simulation networks (easy for me to communicate with those simulations).

A body is represented in space by its location vector from the world axis and then the rotation from the world to the body axis.

Using the Rotation Code

The Body class was modified to contain an orientation with the Quaternions class and an angular velocity with the Euler angles class:


class BodyModel
{
   .
   .
   Quaternion orientation_quaternions;
   Euler angularVelocity_radiansPerSecond;
};

In the body's constructor, I clear out the orientation and the angular velocity:


this->orientation_quaternions.clear();

this->angularVelocity_radiansPerSecond.clear();

Quaternion.clear() sets the Quaternions values to: w=1, x=0, y=0, and z=0 which is the equivalent to an Euler angle of 0,0,0 (no rotation).

I'm using SFML (Simple Fast Media Library) to get joystick inputs to control one of the bodies. Each frame I get the inputs and store them in the orbital aero object. SFML returns axis positions +/- 100, and I'm using a 10% dead zone so I don't have excess twitching when things should be static.


sf::Joystick::update();
if (sf::Joystick::isConnected(0))
{
   float joystickDeadZone = 10;
  float joystickX = sf::Joystick::getAxisPosition(0, sf::Joystick::X);
  float joystickY = sf::Joystick::getAxisPosition(0, sf::Joystick::Y);
  float joystickR = sf::Joystick::getAxisPosition(0, sf::Joystick::R);
  float joystickZ = sf::Joystick::getAxisPosition(0, sf::Joystick::Z);
  if (abs(joystickZ) > joystickDeadZone)
  {
  orbitalAero.joystickInputs.roll = (joystickZ - (joystickDeadZone * sgn(joystickZ))) / (100 - joystickDeadZone);
  }
  else
  {
orbitalAero.joystickInputs.roll = 0.0;
  }
  if (abs(joystickR) > joystickDeadZone)
  {
  orbitalAero.joystickInputs.pitch = (joystickR - (joystickDeadZone * sgn(joystickR))) / (100 - joystickDeadZone);
  }
  else
  {
  orbitalAero.joystickInputs.pitch = 0.0;
  }
  if (abs(joystickX) > joystickDeadZone)
  {
      orbitalAero.joystickInputs.yaw = (joystickX - (joystickDeadZone * sgn(joystickX))) / (100 - joystickDeadZone);
  }
  else
  {
  orbitalAero.joystickInputs.yaw = 0.0;
  }

  if (abs(joystickY) > joystickDeadZone)
  {
  orbitalAero.joystickInputs.thrust = -(joystickY - (joystickDeadZone * sgn(joystickY))) / (100 - joystickDeadZone);
  }
  else
  {
  orbitalAero.joystickInputs.thrust = 0.0;
  }
} // Joystick::isConnected(0)

Then in the body's control system called each frame (not described yet), I assign angular velocities and thrust based on the joystick inputs. The following / 10.0 is just a scaling so the body doesn't rotate too quickly for a given input.

hostBody.angularVelocity_radiansPerSecond.phi = hostOrbitalAero.joystickInputs.roll / 10.0;

hostBody.angularVelocity_radiansPerSecond.theta = hostOrbitalAero.joystickInputs.pitch / 10.0;

hostBody.angularVelocity_radiansPerSecond.psi = hostOrbitalAero.joystickInputs.yaw / 10.0;

Currently I'm changing the angular velocities directly based on the joystick, but in the near future I'll model some type of thrusters / other mechanisms to give an off-axis thrust that will cause the body to rotate. I need to get moment of inertia in for that!

Joystick thrust is considered to be along the body x (forward) axis; however for the processForces calculations later I need that thrust along the world axis.

The following code converts the thrust from the body axis to the world axis via the .rotatedBy call. (If the two axis are aligned, then those vectors will be the same. But if the body is turned at some angle relative to the world axis, then the body's thrust vector will also be rotated by the same angle.)

double thrust_newtons = hostOrbitalAero.joystickInputs.thrust * 1000.0;

// All thrust is along the x (forward) axis of the entity.

// Convert thrust from local axis to world axis.

Vector thrustBody_newtons(thrust_newtons, 0, 0);

hostBody.forceThrust_newtons = thrustBody_newtons.rotatedBy(hostBody.orientation_quaternions.inverse());

// Alternate way to rotate the vector using Quaternions converted to Euler angles.

//hostBody.forceThrust_newtons = thrustBody_newtons.rotatedBy(hostBody.orientation_quaternions.inverse().toEulerAngles());

When summing the forces on a body, my summation starts with the thrust force (relative to world axis) and then adds the effects of gravity and atmosphere from all the other bodies (also relative to world axis).

// Function to loop through all bodies and sum the forces acting on them, including forces from
// thrust, gravity, atmospheric drag, and friction.
bool OrbitalAeroModel::processForces()
{
// Loop through all bodies and calculate their forces upon each other.

// The outer loop is for the primary body that we're calculating forces for.
map <uint64, BodyModel>::iterator iteratorPrimaryBody = tableBodies.begin();
while (iteratorPrimaryBody != tableBodies.end())
{
BodyModel &primaryBody = iteratorPrimaryBody->second;

// Only process the primary body if it's active.
if (primaryBody.active)
{
// Create a vector to contain the sum of all forces.
Vector sumForces_newtons;

// Start with the thrust force of any onboard engines.
sumForces_newtons = primaryBody.forceThrust_newtons;
.
.

And finally the angular velocities change the orientation in the processKinematics loop.

//  Translate and rotate bodies based on their mass, moments of inertia, and forces.

bool OrbitalAeroModel::processKinematics(double timeStep_seconds)

{

// Loop through all bodies and reposition them based on the forces.

map <uint64, BodyModel>::iterator iteratorPrimaryBody = tableBodies.begin();

while (iteratorPrimaryBody != tableBodies.end())

{

BodyModel &primaryBody = iteratorPrimaryBody->second;

if (primaryBody.active)

{

.

.

// Rotation

if (primaryBody.angularVelocity_radiansPerSecond.isNonZero())

{

// Rotation around local axis.

primaryBody.orientation_quaternions.rotateAboutLocalAxis(primaryBody.angularVelocity_radiansPerSecond * timeStep_seconds);

}

} // (primaryBody.active)

iteratorPrimaryBody++;

} // (iteratorPrimaryBody != tableBodies.end())

return true;

}

The angular velocities are relative to the local axis of the entity. So when I do my Quaternion rotation, I call the rotateAboutLocalAxis function (as opposed to rotateAboutWorldAxis).

Testing

I've gone away from my .NET app to visualize the bodies in motion. It worked by having a .NET program consume a CLI wrapper library of my native C++ aero library. Too complicated for something that wasn't ideal long term anyways. Soon I'll be putting together a much nicer 3-D viewer. In the mean time, I've been outputting my bodies via network packets to 3-D visualization software - I'll discuss that in another post.

It has been fun modeling a spacecraft over the Earth and flying it around with the joystick. I'll try to get a video up.

Rotation Classes

Rather than overload this post with both the approach and class code, I'll post my Vector, Euler, Quaternion, and Matrix code in the next few posts.

Copyright (c) 2017 Clinton Kam
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.