HIGH FIDELITY, REAL-TIME OBJECT BODY STATE GENERATOR WITH DATA- DRIVEN MODELING BACKGROUND Systems for simulating movement of objects have become widely used both for consumer entertainment and for professional uses. For example, some video games for consumers simulate a flight of an aircraft by rendering a three-dimensional space. Some other flight simulator systems are available for training pilots in professional settings with enhanced features. While these systems provide at least near real-time motion of the object (e.g., an aircraft) in motion by periodically updating a display of the object, differences arise from accuracy in rendering positions and directions over time. There are tradeoffs at least between accuracy of determining motion of the object and computing resources needed to achieve at least the near real-time processing. As will be discussed further below, determining and updating movement (e.g., positioning and directions) of an object is a critical task in simulating motion of the object. Accordingly there has been a need to improve accuracy while also enabling real-time processing of movement of objects. It is with respect to these and other general considerations that the aspects disclosed herein have been made. Also, although relatively specific problems may be discussed, it should be understood that the examples should not be limited to solving the specific problems identified in the background or elsewhere in this disclosure. SUMMARY As described further herein, there are several challenges for simulating movement of objects. Definition of rigid body, local and global positioning, incorporating impacts of external force into the movement, and errors caused by estimating a trajectory based on interpolation. The present disclosure relates to simulating movement of an object (e.g., a rigid body) along a trajectory as specified by input. For example, the object includes an aircraft, a car, and the like. In particular, the disclosed technology receives customer data that is associated with an object. For example, the customer data includes a mass, an inertia tensor, and force/moment co-efficient associated with a particular model of an aircraft. The disclosed technology generates a rigid body model based on the customer data. The rigid body model includes a combination of three rotational data based on Euler angles (e.g., 312 Intrinsic Euler) and three translational data (e.g., x, y, and z) which together compose a Six Degrees of Freedom dynamics system with inertial properties. The rigid body includes customizable external force (e.g., Newton-Euler Body Equations of Motion). The rigid body further indicates a constant inertial tensor for constant inversion calculation. The present disclosure receives input command (e.g., a velocity of an object, an altitude, heading, and time tick at a particular time (e.g., i-th time tick). Based on input command, the disclosed technology determines position, velocity, and acceleration of the object (e.g., the aircraft) at the time-tick. In particular, a tracking controller, based on nonlinear inversion dynamics (NDI) formulation, generates an angle of attack (AoA), sideslip, roll data, and the body forces (thrust, lift, drag, and local weight) of a rigid body based on a set of input velocity, altitude, and heading of the rigid body. A “slow rotational dynamics” NDI controller generates angular rate commands (e.g., Angular velocities P, R, Q) based on the AoA, the sideslip, and roll data. A “fast rotation dynamics” NDI controller generates moment commands (e.g., Roll, Pitch, and Yaw) based on the angular rate commands. In combination of the moment commands with the thrust, the disclosed technology determines a body state at the time-tick. The body state includes a position, a velocity, and acceleration of the object (e.g., an aircraft) at the specified time. In particular, the disclosed technology includes a real-time strategy (RTS) for formulating rotation using the 312 Intrinsic Euler as the present disclosure converts traditional reference systems and angles into the 312 Intrinsic Euler. The disclosed technology uses an Equations of Motion (EOM) expressed for a body frame, defining forces and inertia and quaternion orientation. In aspects, coordinates associated with RTS may be based on, but not limited to, a Y forward, Z up coordinate system. In describing the present disclosure, the term “physics” of an object refers to physical aspects of the object including motion dynamics of the object and effects of external forces (e.g., wind) upon the object. Physical aspects of an object may include but not limited to position, velocity, acceleration, altitude information, angular rate, moment, and object body forces. This Summary is provided to introduce a selection of concepts in a simplified form, which is further described below in the Detailed Description. This Summary is not intended to identify key features or essential features of the claimed subject matter, nor is it intended to be used to limit the scope of the claimed subject matter. Additional aspects, features, and/or advantages of examples will be set forth in part in the following description and, in part, will be apparent from the description, or may be learned by practice of the disclosure. BRIEF DESCRIPTIONS OF THE DRAWINGS Non-limiting and non-exhaustive examples are described with reference to the following figures. FIG.1 illustrates an overview of an example system implementing a physics simulator accordance with aspects of the present disclosure. FIG.2 illustrates an overview of an example rigid body data model in accordance with aspects of the present disclosure. FIG.3 illustrates an overview of an example body data model including rotation and rigid body data in accordance with aspects of the present disclosure. FIG. 4 illustrates an example method for iterative methods for approximating body data in accordance with aspects of the present disclosure. FIG.5 illustrates an example customer data in accordance with aspects of the present disclosure. FIG. 6 illustrates an example method for simulating physics of an object in accordance with aspects of the present disclosure. FIG.7 is a block diagram illustrating example physical components of a computing device with which aspects of the disclosure may be practiced. FIG.8A is a simplified diagram of a mobile computing device with which aspects of the present disclosure may be practiced. FIG.8B is another simplified block diagram of a mobile computing device with which aspects of the present disclosure may be practiced. DETAILED DESCRIPTION Various aspects of the disclosure are described more fully below with reference to the accompanying drawings, which form a part hereof, and which show specific example aspects. However, different aspects of the disclosure may be implemented in many different ways and should not be construed as limited to the aspects set forth herein; rather, these aspects are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the aspects to those skilled in the art. Practicing aspects may be as methods, systems, or devices. Accordingly, aspects may take the form of a hardware implementation, an entirely software implementation or an implementation combining software and hardware aspects. The following detailed description is, therefore, not to be taken in a limiting sense. This Detailed Description section of Specification references the following document as an Appendix: “Mathematics of the Physics Engine.” The present disclosure relates to simulating movement of an object (e.g., a rigid body) along a trajectory as specified by input. For example, the object includes an aircraft, a car, and the like. In particular, the disclosed technology generates a rigid body model based on sets of customer data that represent the object. The customer data may include a set of exemplary data associated with a type of an object, including a mass, an inertia tensor, and force/moment co-efficient associated with the type of the object. For example, a type object as represented by the customer data may be a particular model of an aircraft (e.g., Airbus A320). In aspects, the rigid body model includes a combination of three rotational data of the object expressed in Euler angles (e.g., 312 Intrinsic Euler) and three translational data of the object (e.g., x, y, z) composing a Six Degrees of Freedom dynamics system with inertial properties associated with the object. In aspects, the rigid body data accommodate data that represent customizable external force (e.g., Newton-Euler Body Equations of Motion). The rigid body model further indicates a constant inertial tensor for constant inversion calculation. FIG.1 illustrates an overview of an example physics simulation system 100 in which a motion of a three-dimensional object may be implemented in accordance with the aspects of the present disclosure. The example physics simulation system 100 includes input 101, an input receiver 102, an NDI tracking controller 104, an NDI slow rotational dynamics controller 106, an NDI fast rotational dynamics controller 108, a rigid body model 110, and a display/transmitter 114. In aspects, input 101 includes global position of the object, global linear velocity of the object, and orientation of the object. In aspects, the input receiver 102 receives input data associated with simulating physical aspects of an object. In aspects, the input receiver 102 (ControlLawD) may receive input data from an artificial intelligence processor that generates a trajectory data based on a predefined and/or a predicted set of data for simulating a motion of the object. Additionally, or alternatively, the input receiver 102 may receive input data as a command for generating a body state at a specific time-tick. In aspects, the input receiver 102 provides data including a velocity 120, an altitude 122, a heading 124 of the object, and a time-tick 126. The velocity may represent a linear velocity of the object (e.g., an aircraft flying at Mach 2.0). The altitude 122 may indicate an altitude of the flying object at a particular time (e.g., at 10,000 meters). The heading 124 may indicate a heading direction of the object (e.g., North). The time-tick 126 specifies a time tick (e.g., a time instance and/or i-th time) of the physical aspects (e.g., physics) of the object in motion. In aspects, an interval of the time- tick may be 1/30 second (e.g., thirty time ticks per second). The example physics simulation system 100 may generate a body state of an object for a given time-tick based on a given set of the velocity 120, the altitude 122, the heading 124 and the current state of the rigid body model 110. The NDI tracking controller 104 (AirplaneControllerD) determines attitude commands associated with the object based on a combination of input, the velocity 120, the altitude, the heading 124, and the time-tick 125. In aspects, the attitude commands as output includes an angle of attack (AoA), a sideslip 132, a roll 134, and object body forces (including thrust, lift, drag and local weight) 128 of the object (e.g., the aircraft). The object body forces may be in a vector expression. In aspects, the NDI tracking controller 104 is based on a non-linear dynamic inversion (NDI) of data. The NDI slow rotational dynamics controller 106 (AirplaneControllerD) determines a set of angular rate commands (e.g., P 140, Q 142, and R 144) based on the given angle of attack, sideslip and roll commands. Processing by the NDI slow rotational dynamics controller 106 includes a non-linear dynamic inversion (NDI) of data. The NDI fast rotational dynamics controller 108 (AirplaneControllerD) determines moment commands of the object in motion based on the received angular rate commands (e.g., P 140, Q 142, R 144). The moment commands may include a roll torque 152, a pitch torque 154, and a yaw torque 156 applied to the flying object. Processing by the NDI fast rotational dynamics controller 108 includes a non-linear dynamic inversion (NDI) of data. The rigid body model 110 (RigidBodyD) represents physical characteristics of the object and its instantaneous dynamics conditions (e.g., states). In aspects, the rigid body model 110 includes a mass of the object, an inertia tensor of the object, force / moment vectors associated with the object. In aspects, a combination of the NDI tracking controller 104, the NDI slow rotational dynamics controller 106, and the NDI fast rotational dynamics controller 108 determines forces applied to the object for determining an updated body date of the object. In some traditional systems, a controller system that determines external forces may be distinct from a system that models a rigid body of the object. The present disclosure integrates the combination of three controllers (e.g., the NDI tracking controller 104, the NDI slow rotational dynamics controller 106, and the NDI fast rotational dynamics controller 108) with the rigid body model 110 for determining a body state at a specified time based on the time-tick. The body state 112 represents output physics data associated with the object at a time as specified as the i-th time tick. In aspects, the body state 112 indicates a position of the object, a velocity of the object, and an acceleration of the object at the specified time. In aspects, an interval between two time ticks may be 1/30 seconds. The display/transmitter 114 displays and/or transmits the body state 112 of the object. In aspects, the display/transmitter 114 displays and/or transmits the body state 112 in numerical expressions and/or graphs. In some other aspects, the display/transmitter 114 displays and/or transmits the object (e.g., the aircraft) in a three-dimensional graphics that represents the body state. In aspects, the display/transmitter 114 transmits the body state 112 to another display (not shown) or a database over a network (not shown) for storage. In aspects, the system 100 may be implemented as in-core, real time simulator, without relying upon data from external computational resources in iteratively determining body states of the object according to the rigid body model. The system 100 may iteratively determine body states based on a time-tick that is incremented as input without interpolating data. Some traditional systems interpolated body states based on a trajectory, which often result in accumulating error in body states as time lapses and the object moves. The present disclosure determines a particular body state that is associated with a time-tick as specified, thereby minimizing errors caused by interpolation in the traditional systems. As will be appreciated, the various methods, devices, applications, features, etc., described with respect to FIG. 1 are not intended to limit the system 100 to being performed by the particular hardware and/or software described. Accordingly, additional controller configurations may be used to practice the methods and systems herein and/or features and applications described may be excluded without departing from the methods and systems disclosed herein. FIG. 2 illustrates an overview of an example rigid data model in accordance with the present disclosure. As illustrated, a system 200 includes a rigid body data model 206. The rigid body data model 206 receives a command to initialize 202 and a time-tick 204. Receiving the command to initialize causes the rigid body data model 206 to initialize itself by resetting data in the rigid body data model 206. In aspects, the rigid body data model 206 (RigidBodyD) is a numerical solver including a rotation formation data 212 and translation formation data 214. In aspects, the rotation formation data 212 use the 312-intrinsic Euler for determining rotational data of the object using unity quatrain Euler angles. . The translation formation data 214 is represented based on three linear states. The combination of the rotation formation data 212 and the translation formation data 214 form the six-degrees of freedom system. In aspects, the rigid body data model 206 works in co-simulation with the AirplaneControllerD (e.g., a combination of the NDI tracking Controller 104, the NDI slow rotational dynamics controller 106, and the NDI fast rotational dynamics controller 108 as shown in FIG.1) to determine at each tick the proper forces and moments which need to be applied to the rigid body to follow the trajectory from input received by the input receiver 102. In aspects, these forces and moments include a lag response filter for dynamic response of control interfaces, thereby improving process performance and generalizing modelling of the controller from object data that are specific to a particular type of object (e.g., generalizing from a rigid body model that is specific to aircraft). The rigid body data model 206 determines a body state of the object at a specified time-tick (e.g., at i-th time tick). The body state 210 may include a position of the object, a velocity of the object, and an acceleration of the object. In aspects, the disclosed technology receives an object data 208 for generating the rigid body data model 206. The object data 208 (e.g., customer data) includes a set of configuration data that describes physical characteristics of the object. For example, the set of configuration data may include a mass information associated with the object (e.g., the aircraft), an inertia tensor of the object, and one or more coefficients associated with an external force and/or moment upon the object. In some aspects, the configuration data may be expressed using Java Script Object Notation. In the disclosed technology includes force / moment coefficients as a part of the object data 208. In aspects, the rigid body data model 206 includes a set of inertia properties of the object. The rigid body data model 206 may further includes the Newton-Euler Body Equations of Motion (EOM) with additional parameters that represent a customizable external force onto the object. The rigid body data model 206 may include conversion of data between the Newtonian reference system and the Euler reference system. In aspects, the rigid body data model 206 maintains a value of inertia tensor at a predetermined constant value for maintaining constant inversion to reduce computational steps needed to determine an updated body state as affected by forces. In some aspects, the object data 208 includes data that represent the integral total effect of the moment used for rotating the object, without data that are specific to a particular type of object (e.g., surface deflection data of an aircraft). As will be appreciated, the various methods, devices, applications, features, etc., described with respect to FIG. 2 are not intended to limit the system 200 to being performed by the particular hardware and/or software described. Accordingly, additional controller configurations may be used to practice the methods and systems herein and/or features and applications described may be excluded without departing from the methods and systems disclosed herein. FIG. 3 illustrates an overview of an example process containers in accordance with the present disclosure. As illustrated, a list of process containers 300 includes a plurality of containers for processing data for collectively determining a body state of an object based on a set of specific conditions at a specific time. EulerAngleD 304 converts rotational data into Euler angles in a reference system. Containers Vec3D 306, Mat3D 308, and QuatD 310 collectively generate vectors and transform data between the Cartesian and unitary quaternion reference systems. RigidBodyD 314 (e.g., the rigid body model 110 as shown in FIG. 1) represents a rigid body associated with the object (e.g., an aircraft) in six degrees of freedom and rotational information in quaternion. ControllerBaseD 316 (e.g., the base framework as shown in FIG.1) executes as a base controller framework for interfacing the rigid body model 110 with the three NDI controllers (e.g., the NDI tracking Controller 104, the NDI slow rotational dynamics controller 106, and the NDI fast rotational dynamics controller 108 as shown in FIG.1) for generating motion data associated with the object. AirplaneControllerD 318 (e.g., the three NDI controllers: the NDI tracking Controller 104, the NDI slow rotational dynamics controller 106, and the NDI fast rotational dynamics controller 108 as shown in FIG.1) processes motion data that are specific to an airplane as an object. ControlLawD 320 (e.g., the input receiver 102 as shown in FIG. 1) processes control laws. In aspects, the control laws may specify how to follow a trajectory and specifies forces applied to the object during its motion. In aspects, program instruction codes used in the respective process containers 302 may be optimized for execution using specific compilers for a particular programming language (e.g., C++ programming language). The optimization may include but not limited to matrix multiplication, quaternion multiplication, exponential multiplication, and the like. In aspects, the optimized programming instruction codes helps determining forces and body states of the object at a specific time-tick in real-time. Some traditional systems determine global positions of an object by integrating velocity in the global reference system. Both the linear and the angular velocity are related to the linear and angular momentum of the dynamic equations. The disclosed technology uses the local body reference system with a local linear and angular velocity data. The disclosed technology determines the global position of the body using rotational orientation based on quaternion. The integration processing data using the local referencing system and the translation into the global referencing system enables the disclosed technology to efficiently determine the global positions of the object. As will be appreciated, the various methods, devices, applications, features, etc., described with respect to FIG.3 are not intended to limit the list of process containers 300 to being performed by the particular hardware and/or software described. Accordingly, additional controller configurations may be used to practice the methods and systems herein and/or features and applications described may be excluded without departing from the methods and systems disclosed herein. FIG. 4 illustrates an overview of an example computational processing in accordance with the present disclosure. As illustrated, a system 400 includes a set of processors for determining derivative values associated with state data of an object for determining a body state of an object based on a set of specific conditions at a specific time. The system 400 includes a current body state 402 of an object (e.g., an aircraft) for simulating physics of the object at a specific time. The current body state 402 includes a RK4 state 404 and a RKMK4 state 406 of the object. In aspects, the RK4 state 404 includes a combination of a global position 410, a local body linear velocity 412, and a local body angular velocity 414. The RK4 state 404 (Euler integration) collects three states of the rigid body. The global position 410 indicates a position of the body in a reference system in a global frame of reference which is inertial. The local body linear velocity 412 indicates a velocity of the body in a linear reference system attached to the body in motion (e.g., lateral x, horizontal y, vertical z). The local body angular velocity 414 indicates a velocity of the body in an angular reference system attached with the body in motion (e.g., pitch along x-axis, roll along y-axis, and yaw along z-axis). In aspects, the RKMK4 state 406 indicates the orientation of the local body reference system with respect to the global inertial reference system. Rotational data that indicate the RKMK4 state 406 include a quaternion orientation 416. The quaternion orientation 416 indicates a rotational state of the body of the object in unit quaternions, representing spatial orientation and rotation of the body in a three-dimensional space. In aspects, the present disclosure uses a combination of a first order Euler 420 and an exponential integration 422 for determining numerical results of a body state at a specified time (e.g., an i-th time tick). The combined operation that processes linear data and exponential data in parallel in distinct methods improves the overall performance of determining a body state in simulating physics of the object. In particular, the disclosed technology processes the RK4 state 404 of the body based on the first order Euler 420. In some aspects, the disclosed technology uses the 4th-order Runge-Kutta (RK4) for processing the RK4 state 404 (e.g., the combination of the global position 410, the local body linear velocity 412, and the local body angular velocity 414) of the body. In aspects, the disclosed technology processes the RKMK4 state 406 of the body based on the exponential integration 422. For example, the exponential integration 422 includes the 4th-order Runge-Kutta Munthe-Kaas (RKMK4) for processing the quaternion orientation 416 of the body. In some aspects, the operation using RKMK4 is more process intensive than an operation using RK4 because of computing exponential data. In aspects, the disclosed technology limits use of RKMK4 to the RKMK4 state 406 and use RK4 for processing the RK4 state 404, thereby optimizing computing resources needed to determine numerical results of the body state at a specific time. In aspects, the disclosed technology maintains unitary of quaternion vector for representing a rotation of the object. In aspects, the disclosed technology may use the Hamiltonian quaternion algebraic operations for processing vectors and quaternion. In some aspects, a Lie group algebraic operation may be used for computing exponential updates for RKMK4. In aspects, the disclosed technology may determine forces in the local reference frame, not in the global reference frame and avoid needs for computing inversion of data matrices and reduce a burden in computational resources. As will be appreciated, the various methods, devices, applications, features, etc., described with respect to FIG.4 are not intended to limit the list of process containers 300 to being performed by the particular hardware and/or software described. Accordingly, additional controller configurations may be used to practice the methods and systems herein and/or features and applications described may be excluded without departing from the methods and systems disclosed herein. FIG.5 illustrates an example object data (e.g., customer data) in accordance with aspects of the present disclosure. The object data 500 includes an exemplary data associated with a particular type of an aircraft (e.g., A-320 aircraft). In aspects, the object data includes data that represents a mass of the object (e.g., 73500 kg), an altitude (e.g., 10972.8 meters), a speed (e.g., 231.4998 meters per second), a wing area (e.g., 318.2421852 square meters), a density (e.g., 0.371845818 kg per square meter), a speed sound (e.g., 295.9081888 meters/second), Mach (e.g., 0.782336579), a gravity (e.g., 9.81 meters per second squared), CL (e.g., 0.227386383), and an Angle of Attack (1.000001091 degrees). In aspects, the object data 500 may further include data that represent a lift-draft curve with retracted and un-retracted landing gear (flaps not deflected). In some aspects, the object data may exclude data associated with surface deflection. In aspects, the lift-draft curve of the object may be expressed in a two-dimensional mapping between lift coefficient values and drag coefficient values. The lift-drag curve data may further indicate an angle of attack for zero lift, an optimum angle of attack, and a critical angle of attack. In aspects, the object data 500 may be specific to a particular type of an object (e.g., IL-14 aircraft, A320 aircraft, and the like). The disclosed technology in the present disclosure generates a rigid body model based on the object data and determine a body state at a specified time based on the input parameters (e.g., a velocity, an altitude, a heading, and a time-tick) by minimizing errors and computational resources needed for processing in a real time. In aspects, the object data 500 includes data that represent three linear momentum (e.g., a set of thrust, lift, and draft, and/or a set of yaw, roll, and pitch), without data associated with surface deflection. In aspects, the present disclosure enables simulating motion of an object without processing actuator and engine dynamics logics. Accordingly, the disclosed technology determines the integral total effect of the moment associated with the object, which rotates the object. The lack of needs for processing logics associated with the actuator and engine dynamics enables generalizing the rigid body model from models that are specific to particular objects (e.g., a model that is specific to an A320 aircraft). In aspects, the disclosed technology receives and use variations of the object data 500 that represents an object for simulating physical aspects of the object. As will be appreciated, the various methods, devices, applications, features, etc., described with respect to FIG.5 are not intended to limit the list of process containers 300 to being performed by the particular hardware and/or software described. Accordingly, additional controller configurations may be used to practice the methods and systems herein and/or features and applications described may be excluded without departing from the methods and systems disclosed herein. FIG.6 illustrates an example method for simulating physics of a body of an object in accordance with aspects of the present disclosure. A general order of the operations for the method 600 is shown in FIG. 6. Generally, the method 600 begins with start operation 602 and ends with end operation 622. The method 600 may include more or fewer steps or may arrange the order of the steps differently than those shown in FIG.6. The method 600 can be executed as a set of computer- executable instructions executed by a computer system and encoded or stored on a computer readable medium. Further, the method 600 can be performed by gates or circuits associated with a processor, an ASIC, an FPGA, a SOC or other hardware device. Following start operation 602, a generate rigid body model operation 604 generates a rigid body model. In aspects, the generate rigid body model operation 604 generates the rigid body model based on an object model that represents physical characteristics of an object. For example, a rigid body model (e.g., the rigid body model 110 as shown in FIG.1) may be based on the object data (e.g., the object data 116 as shown in FIG.1 and/or the object data 500 as shown in FIG.5). In aspects, the generate rigid body model operation 604 may include vectorization and/or transformation of parameter values from the object data from one reference system to another reference system to reduce needs for the overall computational resources to determine a body state of the object at a specified time. In aspects, the generate operation further includes initializing the rigid body model. Initializing the rigid body model may include setting an initial position and, if any, an initial movement of the object. In aspects, determining and updating a body state of the object may be based on the initial position and movement of the object. A receive operation 606 receives a set of input commands including time-tick. In aspects, the input commands include a velocity, an altitude, a heading, and a time as conditions for determining a body state of the object. In aspects, the time may be expressed as a time-tick. For example, an interval of the time-tick may be 1/30 of a second. A generate attitude commands operation 608 generates attitude commands based on the received set of input commands for determining a body state of the object. In aspects, the attitude commands (e.g., angular commands) include a combination of an angle of attack, a sideslip, and a roll associated with tracking the object in motion. In aspects, an NDI tracking controller (e.g., the NDI tracking controller 104 as shown in FIG.1) generates the attitude commands based on a longitudinal reference system. In aspects, the generated attitude commands operation 608 further determines a thrust (e.g., the object body forces 128) associated with the body of the object. In aspects, the tracking controller determines forces needed to perform movement of the object from one position to another position according to a trajectory. A generated angular rate commands operation 610 generates angular rate commands based on the altitude commands. In aspects, the angular rate commands include a combination of a P value, a Q value, and an R value (e.g., the P 140, the Q 142, and the R 144 as shown in FIG. 1). Additionally, or alternatively, the angular rate commands may include a roll rate, a pitch rate, and a yaw rate). In aspects, the angular rates represent the Euler angular velocity. In aspects, the generate angular rate commands operation 610 may include a non-linear dynamic inversion (NDI) operation. A generate moment commands operation 612 generates moment commands based on the angular rate commands. In aspects, the moment commands include roll, pitch, and yaw (e.g., the roll torque 152, the pitch torque 154, and the yaw torque 156 as shown in FIG. 1) of the object. In aspects, the generate moment commands operation 612 may include an NDI operation. A determine body state operation 614 determines a body state of the body of the object. In particular, the determine body state operation 614 uses the received body forces data (e.g., the object body forces 128 as shown in FIG. 1; the object body forces 128 may include thrust, lift, drag and weight local) and the generated moment commands on the rigid body model (e.g., the rigid body model 110 as shown in FIG.1) of the object. The rigid body model may include a mass, an inertial tensor, and force/moment coefficient data for determining a body state at a specified time. The body state may include a position, a velocity, and an acceleration of the body of the object. An update parameters operation 616 updates parameters associated with the rigid body model based on the determined body state. The rigid body model with the updated parameters may be used for the next iteration of determining a body state by feeding back the values to one or more controllers (e.g., the NDI tracking controller 104, the NDI slow rotational dynamics controller 106, and the NDI fast rotational dynamics controller 108 as shown in FIG.1). A display operation 618 displays the body state of the object at the time-tick as specified by the input. In aspects, the display operation 618 may display the body state in a list of numerical values and/or a graph, which depict the body state. In some other aspects, the display operation 618 may display the body state of the object by rendering the object (e.g., an aircraft) in a three-dimensional graphics based on the body state. Additionally, or alternatively, the display operation 618 may transmit the determined body state of the object by storing in a database and/or other devices for rendering via a network. A decision operation 620 decides whether this is the last time-tick (and/or whether there a subsequent time-tick) for simulating physics of the object. It the present iteration is not the last time-tick and/or there is a subsequent time-tick for the simulation, the operational step proceeds to the receive operation 606 and executes the next iteration of determining a body state of the object. In aspects, the disclosed technology completes an iteration of determining and updating a body state in 1/30 of a second with a time-tick interval of the 1/30 of a second. In aspects, the disclosed technology determines a body state for a specified time without interpolating a trajectory based on a previous computation of a body state. Rather, the disclosed technology determines values associated with a body state associated with a specified time-tick through the iteration of the method 600 as detailed above. Unlike some traditional systems with interpolating motion data based on a trajectory, the disclosed technology may provide a higher level of accuracy in a real-time (or a near real-time) operation of 1/30 of a second per iteration. As should be appreciated, operations 602-622 are described for purposes of illustrating the present methods and systems and are not intended to limit the disclosure to a particular sequence of steps, e.g., steps may be performed in different order, additional steps may be performed, and disclosed steps may be excluded without departing from the present disclosure. FIG.7 is a block diagram illustrating physical components (e.g., hardware) of a computing device 700 with which aspects of the disclosure may be practiced. The computing device components described below may have computer executable instructions for implementing a rigid body model 726 on a computing device (e.g., a physics simulation system 100), including computer executable instructions for a trajectory tracking controller 720 (e.g., the NDI tracking controller 104 as shown in FIG.1), a slow rotational controller (e.g., the NDI slow rotational dynamics controller 106 as shown in FIG. 1), a fast rotational dynamics controller (e.g., the NDI fast rotational dynamics controller 108 as shown in FIG.1), and the rigid body model 726 (the rigid body model 110 as shown in FIG.1), and the like, that can be executed to implement the methods disclosed herein. In a basic configuration, the computing device 700 may include at least one processing unit 702 and a system memory 704. Depending on the configuration and type of computing device, the system memory 704 may comprise, but is not limited to, volatile storage (e.g., random access memory), non-volatile storage (e.g., read-only memory), flash memory, or any combination of such memories. The system memory 704 may include an operating system 705 and one or more program modules 706 suitable for running the physics simulator. The operating system 705, for example, may be suitable for controlling the operation of the computing device 700. Furthermore, embodiments of the disclosure may be practiced in conjunction with a graphics library, other operating systems, or any other application program and is not limited to any particular application or system. This basic configuration is illustrated in FIG. 7 by those components within a dashed line 708. The computing device 700 may have additional features or functionality. For example, the computing device 700 may also include additional data storage devices (removable and/or non-removable) such as, for example, magnetic disks, optical disks, or tape. Such additional storage is illustrated in FIG.7 by a removable storage device 709 and a non-removable storage device 710. As stated above, a number of program modules and data files may be stored in the system memory 704. While executing on the processing unit 702, the program modules 706 (e.g., the trajectory tracking controller 720, the slow rotational dynamics controller 722, the fast rotational dynamics controller724, and a rigid body model 726) may perform processes including, but not limited to, the aspects, as described herein. Furthermore, embodiments of the disclosure may be practiced in an electrical circuit comprising discrete electronic elements, packaged or integrated electronic chips containing logic gates, a circuit utilizing a microprocessor, or on a single chip containing electronic elements or microprocessors. For example, embodiments of the disclosure may be practiced via a system-on- a-chip (SOC) where each or many of the components illustrated in FIG.7 may be integrated onto a single integrated circuit. Such an SOC device may include one or more processing units, graphics units, communications units, system virtualization units and various application functionality all of which are integrated (or “burned”) onto the chip substrate as a single integrated circuit. When operating via an SOC, the functionality, described herein, with respect to the capability of client to switch protocols may be operated via application-specific logic integrated with other components of the computing device 700 on the single integrated circuit (chip). Embodiments of the disclosure may also be practiced using other technologies capable of performing logical operations such as, for example, AND, OR, and NOT, including but not limited to mechanical, optical, fluidic, and quantum technologies. In addition, embodiments of the disclosure may be practiced within a general-purpose computer or in any other circuits or systems. The computing device 700 may also have one or more input device(s) 712 such as a keyboard, a mouse, a pen, a sound or voice input device, a touch or swipe input device, etc. The output device(s) 714 such as a display, speakers, a printer, etc. may also be included. The aforementioned devices are examples and others may be used. The computing device 700 may include one or more communication connections 716 allowing communications with other computing devices 750. Examples of suitable communication connections 716 include, but are not limited to, radio frequency (RF) transmitter, receiver, and/or transceiver circuitry; universal serial bus (USB), parallel, and/or serial ports. The term computer readable media as used herein may include computer storage media. Computer storage media may include volatile and nonvolatile, removable and non-removable media implemented in any method or technology for storage of information, such as computer readable instructions, data structures, or program modules. The system memory 704, the removable storage device 709, and the non-removable storage device 710 are all computer storage media examples (e.g., memory storage). Computer storage media may include RAM, ROM, electrically erasable read-only memory (EEPROM), flash memory or other memory technology, CD-ROM, digital versatile disks (DVD) or other optical storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other article of manufacture which can be used to store information and which can be accessed by the computing device 700. Any such computer storage media may be part of the computing device 700. Computer storage media does not include a carrier wave or other propagated or modulated data signal. Communication media may be embodied by computer readable instructions, data structures, program modules, or other data in a modulated data signal, such as a carrier wave or other transport mechanism, and includes any information delivery media. The term “modulated data signal” may describe a signal that has one or more characteristics set or changed in such a manner as to encode information in the signal. By way of example, and not limitation, communication media may include wired media such as a wired network or direct-wired connection, and wireless media such as acoustic, radio frequency (RF), infrared, and other wireless media. FIGS. 8A and 8B illustrate a mobile computing device 800, for example, a mobile telephone, a smart phone, wearable computer (such as a smart watch), a tablet computer, a laptop computer, and the like, with which embodiments of the disclosure may be practiced. In some aspects, the client may be a mobile computing device. With reference to FIG. 8A, one aspect of a mobile computing device 800 for implementing the aspects is illustrated. In a basic configuration, the mobile computing device 800 is a handheld computer having both input elements and output elements. The mobile computing device 800 typically includes a display 805 and one or more input buttons 810 that allow the user to enter information into the mobile computing device 800. The display 805 of the mobile computing device 800 may also function as an input device (e.g., a touch screen display). If included, an optional side input element 815 allows further user input. The side input element 815 may be a rotary switch, a button, or any other type of manual input element. In alternative aspects, mobile computing device 800 may incorporate more or less input elements. For example, the display 805 may not be a touch screen in some embodiments. In yet another alternative embodiment, the mobile computing device 800 is a portable phone system, such as a cellular phone. The mobile computing device 800 may also include an optional keypad 835. Optional keypad 835 may be a physical keypad or a “soft” keypad generated on the touch screen display. In various embodiments, the output elements include the display 805 for showing a graphical user interface (GUI), a visual indicator 820 (e.g., a light emitting diode), and/or an audio transducer 825 (e.g., a speaker). In some aspects, the mobile computing device 800 incorporates a vibration transducer for providing the user with tactile feedback. In yet another aspect, the mobile computing device 800 incorporates input and/or output ports, such as an audio input (e.g., a microphone jack), an audio output (e.g., a headphone jack), and a video output (e.g., a HDMI port) for sending signals to or receiving signals from an external device. FIG. 8B is a block diagram illustrating the architecture of one aspect of a mobile computing device. That is, the mobile computing device 800 can incorporate a system (e.g., an architecture) 802 to implement some aspects. In one embodiment, the system 802 is implemented as a “smart phone” capable of running one or more applications. In some aspects, the system 802 is integrated as a computing device, such as an integrated personal digital assistant (PDA) and wireless phone. One or more application programs 866 may be loaded into the memory 862 and run on or in association with the operating system 864. The system 802 also includes a non-volatile storage area 868 within the memory 862. The non-volatile storage area 868 may be used to store persistent information that should not be lost if the system 802 is powered down. The application programs 866 may use and store information in the non-volatile storage area 868, such as email or other messages used by an email application, and the like. A synchronization application (not shown) also resides on the system 802 and is programmed to interact with a corresponding synchronization application resident on a host computer to keep the information stored in the non- volatile storage area 868 synchronized with corresponding information stored at the host computer. As should be appreciated, other applications may be loaded into the memory 862 and run on the mobile computing device 800, including the instructions for implementing the physics simulator. The system 802 has a power supply 870, which may be implemented as one or more batteries. The power supply 870 may further include an external power source, such as an AC adapter or a powered docking cradle that supplements or recharges the batteries. The system 802 may also include a radio interface layer 872 that performs the function of transmitting and receiving radio frequency communications. The radio interface layer 872 facilitates wireless connectivity between the system 802 and the “outside world,” via a communications carrier or service provider. Transmissions to and from the radio interface layer 872 are conducted under control of the operating system 864. In other words, communications received by the radio interface layer 872 may be disseminated to the application programs 866 via the operating system 864, and vice versa. The visual indicator 820 may be used to provide visual notifications, and/or an audio interface 874 may be used for producing audible notifications via an audio transducer 825 (e.g., audio transducer 825 illustrated in FIG.8A). In the illustrated embodiment, the visual indicator 820 is a light emitting diode (LED) and the audio transducer 825 may be a speaker. These devices may be directly coupled to the power supply 870 so that when activated, they remain on for a duration dictated by the notification mechanism even though the processor 860 and other components might shut down for conserving battery power. The LED may be programmed to remain on indefinitely until the user takes action to indicate the powered-on status of the device. The audio interface 874 is used to provide audible signals to and receive audible signals from the user. For example, in addition to being coupled to the audio transducer 825, the audio interface 874 may also be coupled to a microphone to receive audible input, such as to facilitate a telephone conversation. In accordance with embodiments of the present disclosure, the microphone may also serve as an audio sensor to facilitate control of notifications, as will be described below. The system 802 may further include a video interface 876 that enables an operation of peripheral device 830 (e.g., on-board camera) to record still images, video stream, and the like. A mobile computing device 800 implementing the system 802 may have additional features or functionality. For example, the mobile computing device 800 may also include additional data storage devices (removable and/or non-removable) such as, magnetic disks, optical disks, or tape. Such additional storage is illustrated in FIG.8B by the non-volatile storage area 868. Data/information generated or captured by the mobile computing device 800 and stored via the system 802 may be stored locally on the mobile computing device 800, as described above, or the data may be stored on any number of storage media that may be accessed by the device via the radio interface layer 872 or via a wired connection between the mobile computing device 800 and a separate computing device associated with the mobile computing device 800, for example, a server computer in a distributed computing network, such as the Internet. As should be appreciated such data/information may be accessed via the mobile computing device 800 via the radio interface layer 872 or via a distributed computing network. Similarly, such data/information may be readily transferred between computing devices for storage and use according to well-known data/information transfer and storage means, including electronic mail and collaborative data/information sharing systems. As should be appreciated, FIGS.8A and 8B are described for purposes of illustrating the present methods and systems and is not intended to limit the disclosure to a particular sequence of steps or a particular combination of hardware or software components. The description and illustration of one or more aspects provided in this application are not intended to limit or restrict the scope of the disclosure as claimed in any way. The claimed disclosure should not be construed as being limited to any aspect, for example, or detail provided in this application. Regardless of whether shown and described in combination or separately, the various features (both structural and methodological) are intended to be selectively included or omitted to produce an embodiment with a particular set of features. Having been provided with the description and illustration of the present application, one skilled in the art may envision variations, modifications, and alternate aspects falling within the spirit of the broader aspects of the general inventive concept embodied in this application that do not depart from the broader scope of the claimed disclosure. The present disclosure relates to systems and methods for simulating physics of an object, according to at least the examples provided in the sections below. The system comprises a processor; and a memory storing instructions that, when executed by the processor, cause the system for simulating motion of an object to execute a method comprising: generating, based on object data of the object, a rigid body model associated with the object; receiving a set of input commands over time, wherein an input commands in the set of input commands is associated with a time-tick, the time-tick increments over time; generating, based on the set of input commands, a set of moment commands over time; determining, based on a combination of the rigid body model and the set of moment commands, a body state of the object at the time-tick; and causing an iterative display of the body state of the object at the time-tick incrementing over time. The set of input commands further comprises a velocity of the object at the time-tick; an altitude of the object at the time-tick; and a heading of the object at the time-tick. The rigid body model further comprises rigid body data based on six degrees of freedom; and rotation formulation data based on 312-intrinsic Euler. The determined body state further comprises a position of the object at the time-tick; a velocity of the object at the time-tick; and an acceleration of the object at the time- tick. The determining the set of moment commands uses a combination of a first order Euler based on 4th-order Runge-Kutta on linear state data and an exponential integration based on the 4th- order Runge-Kutta Munthe-Kaas for rotational state data, the linear state data includes at least a global position, a local body linear velocity, and a local body angular velocity of the object, and the rotational state data indicates a quaternion orientation of the object. The processor further causes the system to execute a method comprising determining, based on the received set of input commands including the incremented time-tick, a set of attitude commands; determining, based on the set of attitude commands, a set of angular rate commands; and determining, based on the set of angular rate commands, the set of moment commands. The object data includes at least a type of the object, a mass of the object, an inertia tensor of the object, or a moment co-efficient associated with the type of the object. The processor further causes the system to execute a method comprising iteratively updating, based on the determined body state, a parameter associated with the rigid body model. The processor further causes the system to execute a method comprising: transmitting the body state of the object at the incremented time-tick. Another aspect of the technology relates to a computer-implemented method for simulating motion of an object. The method comprises receiving a set of input commands over time, wherein an input commands in the set of input commands is associated with a time-tick, the time-tick increments over time; generating, based on the set of input commands, a set of moment commands over time; determining, based on a combination of a rigid body model and the set of moment commands, a body state of the object at the time-tick; and causing an iterative display of the body state of the object at the time-tick incrementing over time. The set of input commands further comprises a velocity of the object at the time-tick; an altitude of the object at the time-tick; and a heading of the object at the time-tick. The method further comprises generating, based on object data of the object, the rigid body model associated with the object, wherein the rigid body model includes rigid body data based on six degrees of freedom, and rotation formulation data based on 312-intrinsic Euler, and the object data includes at least: a type of the object, a mass of the object, an inertia tensor of the object, or a moment co-efficient associated with the type of the object. The determined body state includes a position of the object at the time-tick, a velocity of the object at the time-tick, and an acceleration of the object at the time-tick. The determining the set of moment commands uses a combination of a first order Euler based on 4th-order Runge-Kutta on linear state data and an exponential integration based on the 4th-order Runge-Kutta Munthe-Kaas for rotational state data, the linear state data include at least a global position, a local body linear velocity, and a local body angular velocity of the object, and the rotational state data indicate a quaternion orientation of the object. The method further comprises determining, based on the received set of input commands including the incremented time-tick, a set of attitude commands; determining, based on the set of attitude commands, a set of angular rate commands; and determining, based on the set of angular rate commands, the set of moment commands. In still further aspects, the technology relates to a device for simulating motion of an object. The device comprises a processor configured to execute a method comprising: generating, based on object data of the object, a rigid body model associated with the object, wherein the object includes an aircraft; receiving a set of input commands over time, wherein an input commands in the set of input commands is associated with a time-tick, the time-tick increments over time; determining, based on the received set of input commands including the incremented time-tick, a set of attitude commands; determining, based on the set of attitude commands, a set of angular rate commands; determining, based on the set of angular rate commands, a set of moment commands; determining, based on a combination of the rigid body model and the set of moment commands, a body state of the object at the time-tick; and causing an iterative display of the body state of the object at the time-tick incrementing over time. The set of input commands further comprises a velocity of the object at the time-tick; an altitude of the object at the time-tick; and a heading of the object at the time-tick. The rigid body model further comprises rigid body data based on six degrees of freedom; and rotation formulation data based on 312-intrinsic Euler. The determined body state comprises a position of the object at the time-tick; a velocity of the object at the time-tick; and an acceleration of the object at the time-tick. The determining the set of moment commands uses a combination of a first order Euler based on 4th-order Runge-Kutta on linear state data and an exponential integration based on the 4th-order Runge-Kutta Munthe-Kaas for rotational state data, the linear state data include at least a global position, a local body linear velocity, and a local body angular velocity of the object, and the rotational state data indicate a quaternion orientation of the object. Any of the one or more above aspects in combination with any other of the one or more aspect. Any of the one or more aspects as described herein.

APPENDIX

“Mathematics of the Physics Engine.”

Mathematics of the Physics Engine

Fabio Croce, PhD Microsoft Corporation

January 2022

Abstract

This paper describes the mathematical derivation used in the physics engine, producing the novel formulation of a six- degrees- of- freedom (6DOFs) rigid body in the real-time strategy (RTS) reference system. The body receives input forces from a nonlinear inverse dynamics (NDI) controller to follow a trajectory imposed by an external user, either artificial intelligence (Al) or a real user. The set of ordinary differential equations (ODE) produced by the rigid body (Aircraft) and controller (Autopilot) are integrated with a mixed explicit algorithm aiming to preserve the consistency of states during the simulation. The rigid body has been parameterized using the commercial aircraft Airbus A320. The input trajectory to follow is designed with a climb from 5,000 m to 10,000 m, a coordinate turn of 120° heading direction, and a subsequent descent back to 5,000 m. Finally, Simulink® Aerospace Blockset™ with MATLAB® Functions™ are used to numerically validate the dynamics of this system by carrying out a comparative simulation between the physics engine and Simulink. This results in a confirmed optimal accuracy of the Aircraft/ Autopilot modeling for provided inertial, geometric, and aerodynamics properties of military and commercial aircraft.

Introduction

The material in this paper is organized into five sections which describe the areas of development of the engine. For each section, the current base technology is described and the novelties/non obviousness of the new formulation are derived in detail. The five areas are referenced as paper sections:

1. RTS Dynamic Simulation Formulation

2. Nonlinear Inversion Dynamics Controller

3. Mixed Euler /Exponential Numerical Integration

4. Airbus A320 Aircraft Input Data

5. Accuracy Validation of Dynamic Response 1 RTS Dynamic Simulation Formulation

This section first describes how the 3D rotation library is implemented in this engine and what parameters are used to describe the orientation of the body in space. Second, the governing equations for the six-degrees-of-freedom (6DOFs) rigid body are established within the reference systems developed in the rotation library. Finally, the externally driven time step synchronization will be discussed from an architectural point of view, while the numerical solution technique adopted to solve the ordinary differential equations (ODE) generated by the rigid body will be discussed.

1.1 Rotation Library

In this derivation, we will be using a right-handed Cartesian reference system defined as follows: x to the right side of the body, y to the forward direction of the body, and z to the upward direction of the body. This choice is consistent with the standard real-time strategy (RTS) reference system used in gaming. Two main Cartesian reference systems will be used during this derivation, that is, a global reference system with capital X, Y, and Z and a local reference system with lower case x, y, and z. Figure 1 depicts those two basic reference systems:

Figure 1: Global and local reference systems.

In the same figure 1, three basic vectors are shown representing the typical reference kinematics of multibody dynamics bodies in space - see [1] “Reference Kinematics” and the references therein. In this derivation, the orientation of the body is composed using three sequential rotations for transforming a vector expressed in global reference to local reference. This sequence, named as 3-1-2 Intrinsic Euler, is formed by a first rotation about the z-axis of the yaw angle ψ , followed by a second rotation about the instantaneous x-axis of the pitch angle 0, followed by a third rotation about the instantaneous y-axis of the roll angle Φ . Thus, the rotation matrix representing the orientation of the body in the global reference system is the following equation (1.1):

(1-1)

Having defined the orientation of the body with respect to the global reference system, one can characterize each point of the body using the translation/rotation formulation of multibody systems (see figure 1), that is in formula:

(1-2)

Formula (1.2) will be used to translate local states into global states within the solver solution process. Note that the nomenclature indicates an over bar for vector/matrix quantities expressed in the local reference system, e.g. u indicates the position of a point P of the body in local coordinates. Thus, equation (1.2) combines the translation and rotation of the rigid body, providing the coordinates of the point P depicted in figure 1 in global coordinates.

Other important relations which need to be derived in the RTS reference system are the so- called G matrices expressed in terms of the Euler angles. These matrices depend upon the sequence that has been chosen for representing the global-to-local body orientation, thus existing [1] cannot be applied. Since infinitesimal rotations are vector quantities, we can reconstruct the angular velocity vector by using the three infinitesimal yaw, pitch, and roll rotations to compute the relationship between Euler angles time derivatives and the global/local angular velocity vector. These equations will be used in the nonlinear dynamics inversion controller formulation:

Once the G matrix is defined, one can easily derive its inverse, the Those are key relationships between the global angular velocity vector ω, the local angular velocity vector and the time derivatives of the Euler angles θ, that is:

Note that the nomenclature of the global angular velocity vector is as a>i while the one for the local angular velocity vector is as p, q, r, which is the traditional solver states notation for this vector.

The rigid body solver utilizes the traditional Euler parameters, known as quaternion, as orientation states, while, as mentioned before, the aircraft controller utilizes Euler angles to describe the body orientation. Note that there is a subtle nomenclature difference between three Euler angles and four Euler parameters - the two sets are mathematically very different quantities, thus refer to [1] to clarify the relationship between them. The reason for this different choice in body orientation description will be evident after reading the controller and numerical integration sections of this paper.

When Euler parameters are used as orientation states, we can formally write analogous relation- ships to (1.3), (1.4), (1.5), (1.6), by substituting the expressions for G matrices and θ accordingly. The advantage of using Euler parameters versus a particular Euler angles sequence is that the rotation matrix A is unique and can be expressed as follows:

By using equations (1.7), (1.8), and (1.9), one can derive the following relationships between angular velocity and Euler parameters time derivatives - this can be helpful to implement a first order Euler-based update in the quanternion states:

Note that the relationship between equations (1.10) - (1.11) and equations (1.12) - (1.13) is not a trivial matrix inversion, but instead requires using E ^T E = E ^T E = I4 — 00 ^T and noticing that

= 0 because the dot product of a constant length vector and its derivative is zero since the two vectors must be perpendicular - the proof is left to the reader for exercise.

It is possible to move from one formulation, i.e. rotation matrix expressed in terms of Euler angles (1.1), to the other, i.e. rotation matrix expressed in terms of Euler parameters (1.7), by using the relationships defined in [2] (Appendix A, 9). Note that we are using a Jet Propulsion Laboratory (JPL) quaternion convention for unit rotations, as defined in [3],

Another important part of this derivation is the interface with standard aerospace conventions to represent the orientation of an airplane. Traditionally, the aerospace industry adopts Euler angles sequence defined as 3-2-1 Intrinsic Euler - see [2] (Appendix A, 10), hence the yaw, pitch, roll angles as well as the rotation matrices or inertia tensor properties defined in the two reference systems are not directly translatable without proper transformations. This second reference system is usually named as North, East, Down (NED) or Tait-Bryan angles (TB), and it is defined as follows: x to the forward direction of the body, y to the right side of the body, and z to the downward direction of the body. As can be seen, this is also a right-handed reference system, since there exists x X y = z.

Figure 2 provides the relative orientation between the two reference systems TB and RTS. It also shows how Euler angles map to different axes depending on which sequence we are adopting.

Figure 2: Real-time strategy (RTS) and Tait-Bryan (TB) reference systems.

Thus, the vector, matrix and orientation need to be properly translated from one reference system to the other.

We begin the derivation noticing that one can find a transformation from a vector expressed in TB to the corresponding in RTS and vice-versa by simply using the relative rotation matrix of one system versus the other, e.g.: where the columns of A _TB→RTS are the direction cosines of TB axis expressed in RT reference system. We should also note that matrix A _TB→RTS is symmetric and orthonormal, hence the following properties exist:

Now we consider a generic vector transformation in TB from local to global, e.g.: which can be expressed in RTS by means of (1.14), that is: and thus:

Hence, starting from a rotation matrix expressed in TB, e.g. RT _B , the corresponding rotation matrix in RTS is expressed as:

In the particular case of TB/RTS transformation, this quadratic form shows symmetric matrices pre and post multiplying. Nonetheless, given the generality of the derivation, this can be applied for non- symmetric, yet orthonormal matrices. We recall that a rotation matrix is by definition orthonormal, e.g. its inverse is equal to its transpose Q ^T Q = QQ ^T = I.

Equation (1.21) has a very important application to the inertia tensor. Because customer data are usually provided in TB reference system, we will be using this same transformation to retrieve the RTS inertia tensor which is necessary for the equations of motion (EOM). Hence, we can generalize the above as follows:

Now we have all the pieces for computing Euler angles RTS from Euler angles TB, as follows:

• create the rotation matrix R _BE from Euler angles TB, by using [2] (Appendix A, 10)

• convert to the corresponding rotation matrix R _RTS using equation (1-21)

• extract Euler angles RTS from R _RTS using [2] (Appendix A, 9)

The same approach can be used to go from RTS to TB. Again, we leave it to the reader to prove this derivation as an exercise.

1.2 Rigid Body Library

This section describes how the equations of motion for the rigid body have been written using the RTS rotation formulation developed. We will develop the 6DOFs equations of motion expressed in the local body reference frame previously described. The details of the numerical integration will be covered in the section describing the mixed Euler /Exponential approach, while in this section we will formally derive the equations to be integrated.

Since Euler parameters, e.g. quaternion, will be used as orientation parameters for the body, one could argue that the correct nomenclature for this rigid body is 7DOFs, three linear states plus four angular states. For consistency with literature and historic nomenclature, we will be using the 6DOFs naming convention. It should be noted that the extra state for the quaternion orientation allows for resolving singularities such as Gimbal lock for pitch angles of 90° positive or negative, by exploiting the fact that: where 0 is the angle of rotation of the rigid body along the unit vector Hence, while the numerical integration of the rotation states will exploit the more robust quaternion implementation, at converged step the graphical client may require to extract the Euler angles for rendering the orientation of the body. If that is the case, again [2] (Appendix A, 9) can be used for this purpose.

In case the graphical client uses left-handed reference systems with sequences different from RTS, an approach similar to the one described for converting Euler angles between RTS and TB can be exploited. Thus, the procedure previously described is general as long as one properly defines the transformation matrix (1.14).

We will be writing the equations of motion (EOM) according to the Newton- Euler approach, which provides three equations for the linear DOFs and three equations for the rotational DOFs. Those equations will be written in compact matrix form to allow the integration with the RTS rotation library for position and orientation of the body.

The second aspect of this derivation is the choice of writing the equations in the global or local reference system. We chose the latter for multiple reasons. First, the expression of forces (especially in the case of the aerodynamics thrust, lift, drag and rotation moments of an aircraft) have a more natural formulation in the body reference system. This aspect will be described at length in the controller portion of this paper. Second, the inertial tensor has a constant expression within the body reference system, equal to the value provided by customers. If the equations were to be written in the global reference system, a re-evaluation of the inertia tensor (and therefore its inverse) should be carried out at each time step integration, according to the following where A is the rotation matrix representing the instantaneous orientation of the body. Finally, the effect of fictitious inertia forces F» due to the fact we are using a rotating frame, can be taken into account as follows - see [1] “Analytic Techniques” and the references therein:

The dynamic equations of motion produce second order ordinary differential equations (ODEs) which need to be integrated numerically. Traditionally, the state vector is defined using a statebased approach, where we lower a degree of differentiation by duplicating the number of states, that is, we introduce the global position x (t) as integration of body linear velocity and the quaternion orientation θ (t) as integration of body angular velocity - note the mixed approach global/local states, where we differentiate the vectors expressed in local reference system by using the over bar. The “3 4 3 3” state vector is therefore defined as follows: and the derivatives with respect to time of equation (1.24) provide the following first order differential equations, whose expression is based on the linear momentum (for linear states) and angular momentum (for angular states):

The equation of linear momentum in body reference coordinate system can be written as follows:

The equation of angular momentum L (t) J co (t) in body reference coordinate system can be written as follows:

At this point, we have all the information for carrying over four numerical integrations over the four vector expressions in equation (1.25). The reader will note that the four expressions are mathematically linked together, therefore we have to compute these numerical integrations sequentially; that is, in a multistage solver approach, within each step the values of the current stage is dependent on the result from the previous stage.

The first equation provides the derivatives with respect to time of the global position:

The second equation provides the derivatives with respect to time of the quaternion orientation: The third equation provides the derivatives with respect to time of the local linear velocity:

The fourth equation provides the derivatives with respect to time of the local angular velocity:

The F/ _oea / and T/,,,,,,/ are also functions of time, and depending on the scenario simulated may have simple or complex expressions. For example, the constant expression of gravity force in global RTS frame needs to be transposed to the local body coordinate system by using the instantaneous body orientation, e.g.:

Other forces such as drag depend on the local body velocity of the aircraft via proper aerodynam- ics coefficients, thus they do not require a transformation via the rotation matrix, but only the knowledge of the body velocity:

A third most important application is the one where the rigid body needs to follow a prescribed trajectory. In such case, a controller needs to compute the forces needed to properly drive the body in the desired trajectory and most likely will work with a co-simulation approach with the rigid body, since the computation of such forces depends upon the current state of the body. This is the use case addressed later on in this paper by implementing an NDI controller trajectory tracking.

However, the derivation developed in this section is modular enough that if forces are known a priori or prescribed according a predefined law, the rigid body solver can be simulated in a standalone fashion. Most importantly, this derivation can be applied to a variety of rigid bodies of physics. For example, by properly customizing the expressions of Fi _ocai and Ti _OC ai , effects of weather conditions, fluid interaction, and electromagnetism can be taken into account in the rigid body solution process.

For more detailed applications to the aerospace industry, two sources can be reviewed. [4] provides Tait-Bryan formulations of aircraft dynamics. [5] overviews simulation techniques used in Microsoft Flight Simulator.

1.3 Externally Driven Time Step Synchronization

The above developed solver provides the dynamic conditions of a body, e.g. its position, velocity and acceleration, which undergoes applied external forces at a certain time. The update state routine receives as input a time step increment to compute the dynamic evolution of states starting from the values at current time. At the very beginning, the process of numerical integration requires initial conditions (ICs). Mathematically speaking, ODE requires a number of initial conditions equal to the number of states which need to be integrated. Therefore, the following set of ICs are provided at the beginning of the simulation:

• global initial position x ₀ (t)

• quaternion initial orientation θ ₀ (t)

• body initial linear velocity

• body initial angular velocit

At this point, assuming that the forces are known at time to, the system is able to compute the dynamic conditions of the body at time ti = t ₀ + dt. This time step could be internally generated for a standalone simulation from the time t _initial to t _final . Traditional computer aided design

(CAE) packages provide the user with simulation results within those extremes of integration. This is usuall ^v known as offline standalone simulation.

Two main problems arise with this type of approach. First, the simulation can be very lengthy depending on the complexity of the model and most importantly cannot be modified while it is running. A prescribed system of forces interacting with the rigid body will determine the timetrace of the states values. This type of simulation is useful for design study or validation for virtual prototyping.

Within the real-time strategy (RTS) requirements, the simulation needs to satisfy two conditions:

1. Real-time execution: ability of compute the t-th state update for the time step dt before the simulation time = t ₀ + dt occurs;

2. Environment interaction: ability to modify the evolution of states at time t > ti following the occurring event at time f;.

In order to achieve both requirements, the code computing the i-th state update needs to run faster than the wall-clock time ti = to + dt. Thus, very complex physics modeling needs to be carefully evaluated and benchmarked against the maximum frame rate engine tick. Second, the forces acting on the body need to update according to the event. Both conditions require therefore a time step synchronization between the game engine tick and the physics engine time step integration process.

Some approaches execute offline simulation and then import data which will be interpolated between two known and already solved simulation times. This approach requires knowledge of the system before an event occurs, hence it loses the flexibility of real-timeness. Moreover, interpolation techniques introduce numerical approximation because they do not solve the dynamics at the required t _i , but compute an interpolated value for t _i between the prior and next solved step.

The architectural design of this solver instead allows for satisfying the above requirements. A current 10 frame-per-second simulation has been conducted with an externally driven time step. This means that the physics engine works on demand with a provided time step from the game engine tick. At a 10 frame per second tick rate, the solver was able to run faster than real time, avoiding any lagging in the simulation output. Moreover, the NDI controller computed the forces per tick in co-simulation with the current state of the rigid body, by evaluating the forces necessary for trajectory tracking with respect to the previous time step state values.

Figure 3: Externally driven time step synchronization.

Figure 3 depicts a schematic of the physics engine interaction with external tick. Based on this schematic, the flow of operations and time step synchronization work as follows:

1. A rigid body is created based on the .ISON file configuration which specifies the geometric, inert ial and physical properties of the body, e.g. aerodynamics coefficients for a known aircraft type.

2. The game engine provides initialization for the four set of states previously specified, e.g. global position, body orientation, local linear velocity, and local angular velocity.

3. The game engine requests a step update providing the delta time step (function of its tick rate).

4. The rigid body integrates the internal states based on the previous converged step (or ICs if first step).

5. The rigid body sets the relevant dynamic states to the owning body, e.g. position, velocity. both linear and angular.

It should be noted that the numerical integration scheme will be chosen according to this time step synchronization. Among implicit, and explicit schemes, an explicit is preferred for several reasons. First, the time execution of an update has a constant wall-time execution because it is not based on a Newton-Raphson integration scheme, but more traditional Euler schemes, where the first derivatives are used to compute the next step.

Without going into detail about the two schemes (the reader can further research the differ- ences) , it is clear that an implicit scheme provides an unpredictable performance execution update depending on the smoothness of the external forces. It also does not guarantee (at least in its original design) a constant time step execution because it is designed to accelerate and slow down depending on the numerical stability of the simulation. Rather, the explicit schemes fit very well in a requirement for real-timeness driven by an external time step, which is exactly the design of figure 3.

The disadvantages of explicit methods - such as divergence of integration error - are greatly mitigated by smaller time step execution dictated by the game engine, which typically runs at 30 frames per second. Moreover, the mixed Euler /Exponential approach discussed later on as integration technique preserves the quaternion unit norm, a basic requirement for it to represent an accurate body rotation.

2 Nonlinear Inversion Dynamics Controller

In this section, we will derive the nonlinear dynamics inversion (NDI) formulation for a trajectory tracking of an aircraft under input commands from an artificial intelligence (Al) module or user inputs. This derivation has as background the work listed in [6], [7], [8]. However, there are some key differentiators from prior literature:

• Existing formulation has been developed in typical aerospace industry Tait-Bryan reference systems. As such, it is not applicable to the rigid body solver derived in the previous section.

• One could potentially use transformations similar to the ones developed in the rotation library to re-use existing derivations. However, this would lead to highly unoptimized algorithms which would miss the real-timeness requirement for this engine. Moreover, the amount of floating point conversions at each solver time step would provide numerical noise in the simulation results.

• The system has been developed to fulfill a data-driven approach for this engine. To fulfill this requirement, we decided to not implement “control mixer/surface deflection” as stated in the abov^mentioned literature, since that approach leads to a non-generic type of algorithm and requires customization per aircraft model. It would also add computation cost, which is again a key differentiator for the real-timeness of this engine.

• The dynamics response provided by surface deflections has been mathematically implemented via lag response filters which are easy tunable to achieve a characterization of the dynamic response of the real aircraft.

• Finally, the goal is to obtain a controller which does not require gain optimization before running, since again this would fight against the real-timeness of the engine requirements. Thus, a set of fixed bandwidth factors between fast, medium and slow dynamics have been used for the controller.

2.1 Derived Formulation in RTS

We begin the derivation by considering the body reference system attached to the aircraft from figure 1. We separate the longitudinal dynamics to the rotational dynamics as follows: the longi- tudinal dynamics is computed in steady-state dynamic equilibrium of forces for a point mass at the center of gravity of the aircraft; the rotation dynamics is computed as rigid body rotation for the single body composing the aircraft. Since the generic motion of a body in space is defined as translation plus rotation, this approximation fulfills the computation of forces necessary to follow a trajectory with translation plus rotation of the aircraft. It must be emphasized that this reduced order modeling (ROM) is only used for providing steady-state controller forces to the rigid body differential equations. The output dynamics of the aircraft remains fully computed as 6DOFs with- out approximation. The following section will provide a detailed description of the co-simulation between controller forces and rigid body solver.

Following [11] and [12], the free body diagram of the point mass system in RTS reference system is depicted in figure 4, with the relevant longitudinal dynamics forces acting on it.

Figure 4: RTS point mass aircraft model with linear dynamics.

The relevant forces responsible for the longitudinal dynamics of the aircraft are:

• Thrust T in the direction of the body linear velocity

• Drag D opposite to the body linear velocity

• Lift L perpendicular to the body linear veloc on a plane rotated of roll angle 0

• Weight W

One can write the longitudinal and vertical equilibrium equation in the RTS y-z plane of the aircraft as follows:

Additionally, we can exploit the constant radius turn equation for planar motion - see [11] (3.1) for the derivation - for a simplified relationship between the yaw angle and roll angle Φ of an aircraft executing a coordinate turn with a climb angle 7:

Within the NDI formulation we need to define the control variables CV (see [12]) which will be the input of our controller. This is an arbitrary choice of the controller designer, and since typically an external Al module or user would want to control an aircraft in terms of speed, an altitude and a heading angle, we will select the following control variables:

(2-3)

We will then select the following state variables from the 6DOFs body which affect the longitudinal dynamics:

(2-4)

The next step is writing a Jacobian matrix able to transform time derivative of the state variables (2.4) into time derivative of the control variables (2.3). By doing so, we will be able to express equations (2.1) and (2.2) in terms of the control variables, which is our ultimate goal. The order of derivation of this Jacobian matrix needs to be sufficient to show an order of derivation equal to desired control variable derivatives

Hence, we will write the following: The reader can verify the above expression, noticing that:

Now we can carry out the following multiplication:

By substituting the equations (2.1) and (2.2), we can derive the expressions of the body reference system forces which need to be applied to the 6DOFs system for input commands as a function of the control variable derivatives.

The first equation for lift L can be is as follows, where W = mg:

The drag D can now be computed as function of lift L, using customer provided polar aerody- namics coefficients (for their expression, refer to the next section of this paper):

Finally, the thrust T can now be computed as follows:

At this point, we can compose the full body force vector expressed in RTS reference system which are the input to the rigid body without any control mixer or surface deflection:

The remaining quantities from this first controller will then be passed in the rotation dynamics section comprised of other two controllers, whose output will provide the three remaining moments in body local RTS reference system.

Following [8] and [10] nomenclature, we can now provide inputs to the slow and fast rotation dynamics. The difference here is that the control derivatives herein derived are meaningless within RTS formulation, hence we need to rc-dcrivc them according to the 6DOF developed prior. Also, as for the forces computed in (2.10), we will use provided moments in the body local reference system which can be directly applied to the rigid body rather than to a control mixer or surface deflection system. Please refer to the previous paragraphs for an explanation of this reasoning. Second, the dynamic response of the aircraft regarding rotation dynamics will be taken into account with a fast lag filtering as explained in [9] . This has been proven sufficient to capture the response dynamics of the aircraft within this modeling technique.

For the second controller, we select the control variables angle of attack α, slip angle β and roll angle 0, which are outputs or computable from the first controller: while the body states selected are the three Euler angles representing the orientation of the body and the three components of the linear body velocity v (t):

The first, step is mapping the output from first controller to equation 12,111, We notice that a kinematic coordinate turn requires maintaining a slip angle null. Second, the angle of attack can be obtained again from the customer data characteristic airfoil as in figure 5 - this can be obtained by either linear interpolation or look up table of the linear region from this plot (m and n coefficients), limiting its validity from —2° to 18°:

Different aircraft may have steeper characteristics or different stall limits. Finally, the roll angle is directly computed as output from the coordinate turn inversion equation (2.2):

Figure 5: Relationship between lift coeflicieut and angle of attack.

Given the provided by the RTS rigid body, the angle of attack α and the side-slip angle β are defined as follows: a = tan" ¹ (~^)

Hence the Jacobian matrix here is expresses] as follows:

We can define the transfer function matrix. G _x from the desired output from this controller, which is the body angular velocity to the derivative of states (2.12). We will recognize here the matrix (1.6) which provides the derivative of Euler angles if applied to the body angular velocity vector, and the equation (1.26) with zero input force

Now we have all the pieces to compose the NDI law of the second controller: which can be inverted to provide the following final transfer function for second controller: where the transfer function T has the following three column vectors: The output from equation (2.19) is then fed into a third controller for fast rotation dynamics. The Jacobian of this final NDI controller is the Identity matrix 3x3, since the control variables are the body angular velocity (computed by second controller) and the output are the moments in body reference frame. Such a relationship is the Euler equation of motion for the RTS rigid body, equation (1.27), in which we use the derivative of control variables as angular accelerations:

2.2 Co-simulation with Rigid Body

The following figure 6 summarizes the series of operations in cascades on the 3 NDI controllers working in co-simulation with the RTS rigid body.

Figure 6: Schematic of the three NDI controllers in cascade.

It is clear from equation (2.6), (2.19), and (2.20) that the three NDI blocks take as input the derivatives of the control variables and not the control variables (CV) themselves. Thus, it is necessary to provide a relationship between the commands CV1 (desired trajectory provided as input to 1 ^st controller), CV2 (1 ^st controller output provided as input to 2 ^nd controller), CV3 (2 ^nd controller output provided as input to 3 ^rd controller), and their respective derivatives, as required by equations (2.6), (2.19), and (2.20).

As we can see from the diagram, the controller receives feedback from the previous iteration states of the rigid body (global position, linear global and local velocity, orientation, angular veloc- ity), to compute current iteration step forces which allow the integration of current iteration state within the rigid body. This is known as co-simulation approach, because of the tight interaction between the history of the two components. Suffice it to say that the time histories of rigid body and controller need to be carefully tied together as described in this paragraph.

Additionally, as in every control strategy, we will add a feedback loop control responsible for closing any error between the current rigid body state iteration and the desired input value at that particular time. This will be part of the input to each NDI controller, and its function is correcting any discrepancy that the model has with the desired input. We will not use any gain optimization technique, but a simple unit proportion gain, since for NDI controllers the dynamic inversion provides a sufficiently accurate guidance to the system, and therefore we need to rely less on proportional- integral-derivative (PID) strategies - see [10], [11] , [12] for comparative results and differentiators from the gain selection choices therein. Following the figure 6 from left, the three inputs are provided as step commands from an external entity, such as Al or user module. A step command is not realistically achievable by any mechanical devices since inertia plays a fundamental role in the dynamic evolution of the system. Thus, the CV1 needs to be translated into smooth and continuous control laws to the controller. Moreover, the derivation of NDI shows that at least a second order differentiation needs to be achieved by the altitude. Hence, we develop an analytical polynomial piecewise expression which should at least support C2 requirements (continuous and derivable up to the second order). This will guarantee one can obtain from a given command the first and second derivative over time. The time execution of the whole command is a characteristic of the aircraft and therefore must be provided as input data, e.g. a commercial aircraft has a turning and climb rate much lower than that of a fighter jet. The knowledge of this rate allows us to compute the maneuver time, e.g. the time from the initial state when the external command started and the final time needed to reach the external command.

We will consider an example of input climb command to show how to create signal command, first derivative of the command and second derivative of the command. Velocity and heading inputs can be mathematically obtained in the same way, although the second derivative is not strictly needed in equation (2.6). The time to execute the total maneuver is computable as:

Knowing the total time of maneuver, we are going to divide it into three segments (they could be equally or differently spaced depending on the style of the autopilot). As mentioned, we need to make sure that a piece-wise definition of the control law still produces continuous derivatives for the signals up to the second order. This will avoid discontinuity on forces definition which would produce spikes in the acceleration of the rigid body integrated states. In order to obtain that, we define the piecewise function as follows:

1. A first smoothing phase between signal _initial and the linear value provided by signal _rate at the end of the first sector. The first and second derivative will be both smooth polynomial during this phase.

2. A second linear phase in which the command is expressed by the linear function with slope equal to signal _rate . The first derivative will be constant, and the second derivative will be null.

3. A third smoothing phase from the linear value provided by signal _rate at the end of the second sector and the value signal _final .he first and second derivative will be both smooth polynomial during this phase.

To guarantee continuity at the connection phase, we define a Hermitian polynomial for the first derivative of the command during phase 1 and 3. In this way, the signal command will be the integral of this function and the second derivative of the command will be the derivative of this polynomial.

We leave to the reader to derive the expression of the linear function during phase 2 - it is a simple linear equation based on simulation time of type y = m * x + n. The analytical expression of the smooth signals used during phase 1 and 3 are expressed below.

For the below example, we assume we want to reach a rate of climb of mi = 15.24 m/s from mo = 0, in an amount of time equal to 1/3 of the total maneuver time. The actual simulation time needs to be linearly mapped between zero and one, as required by the Hermitian interpolation function. The first derivative of the command is expressed by equation (2.21) and represented in figure 7:

Figure 7: First derivative of the command.

By integrating with respect to time equation (2.21), we can then find the corresponding signal command which brings our slope from mo = 0 to m ₁ = 15.24 m/s. The signal command is expressed by equation (2.22) and represented in figure 8:

Figure 8: Signal command.

By deriving with respect to time equation (2.21), we can then find the corresponding second derivative of the command which brings our slope from = 15.24m/s. The second derivative of the command is expressed by equation (2.23) and represented in figure 9:

Figure 9: Second derivative of the command.

We leave to the reader the exercise of implementing phase 3, where we can reuse equations (2.22), (2.21), and (2.23) to compute respectively the signal command, the first derivative of the command, and the second derivative of the command which brings OUT slope from m± = 15.24 m/s to mo = 0.

At this point we have all the pieces to compose the input to first controller (CV1), second controller (CV2), and third controller (CV3). For CV1, the control laws already provide the derivatives of control variables as analytical functions - see equations (2.22), (2.21), and (2.23), so we will be adding only the contribution of the feedback loop for closing the error. That is, with reference to (2.6), CV1 is expressed as: where the subscript cmd indicates the signal command values from (2.22) and the subscript cur indicates the value from the previous converged i-th states from RTS rigid body. One can notice in the expression both the closed loop error contribution with proportional gain equal to 1 and the analytical command derivatives.

For the second controller, with reference to equation (2.19), no analytical expression is available for CV2. Thus, according to standard procedure in control theory, we approximate the derivatives as cmd _value — cmd _state multiplied by a bandwidth. Since the second controller is responsible for the slow rotation dynamics, we will use a lag response system as designed in [9] to capture the dynamic response of “control mixer/surface deflection”. We will use bandwidth gains of 2rad/s for a and /3, and 1.5 rad/s for (b as advised in [8]:

Similarly, for the third controller, CV3 for equation (2.20) will be computed with a bandwidth gain of 10 rad/s to differentiate the fast dynamics of the body angular velocity vector from the slow

In conclusion, we now have all the inputs to run the co-simulation between a 6DOFs rigid body and a series of NDI controllers, receiving the input as control laws from external modules such as Al or user commands. While no external command is acting, the controller still tracks the current states of the rigid body versus the last known desired commands. Thus, this co-simulation package provides a complete way for aircraft autopilots both to maintain a route under effect of random external excitations, such as wind force applied, or to change the relevant flying conditions, e.g. altitude, speed, and/or heading.

For completeness, we will show some results where the three control laws for velocity change (figure 10), altitude change (figure 11), and heading change (figure 12) have been implemented by coding equations (2.22), (2.21), and (2.23) into MATLAB® Functions™. These inputs will enable validation against the results of Simulink® Aerospace Blockset™.

Figure 10: Velocity Control Laws.

Figure 11: Altitude Control Laws.

Figure 12: Heading Control Laws. 3 Mixed Euler /Exponential Numerical Integration

In this section, we will describe the numerical integration scheme which has been adopted for integrating the equations of motion of the rigid body, once the controller has computed the forces from the previous converged step (or from initial condition states if it is the first integration step). We previously recommended the choice of the numerical integration scheme in consideration of the externally driven time step, as 4 ^th order Rung^Kutta. This scheme is widely accepted for its performance and robustness. The novelty to follow will be the approach of mixing a traditional 1 ^st Euler approach with the novel Exponential quaternion integration exploiting the Lie Group algebra. References to these numerical integration schemes can be found in [13], [14], [15], [16].

3.1 Mixed Euler RK4 and Exponential RKMK4

The typical approach to integrate differential equations from time with the explicit algorithm RK4 - Runge-Kutta 4 ^th order - can be summarized as follows: where the coefficients ayj, bj and c,- can be found from the Butcher table for 4 ^th order - see figure 13. This approach is usually implemented for continuous states which have no restriction to their own domain, thus will be applied to solve equations (1.28), (1.30) and (1.31).

It should be clear to the reader that a brute-force application of first Euler integration method to equation (1.29), would violate the integrity condition for quaternion states, which requires the following equation to be satisfied in order to represent a rotation:

Hence, a simple Euler approach for integration of quaternion is mathematically incorrect because it violates the previous condition:

Figure 13 : Explicit coefficients from Butcher table for 4 ^th order Runge-Kutta.

Some approaches try to correct the equation (3.2) by either re- normalizing the quaternion after each update or by forcing a small time step integration which would slightly deviate the quaternion itself from having unit norm. This renormalization can be expensive, and it is not geometrically correct because it provides a projection of the non-normalized quaternion back on the unit circle, deviating therefore from the correct position of the incremented quaternion. Results are shown to be acceptable (see [14] ) only if the time step is sufficiently small:

Instead of being limited by small-time step or losing accuracy by projecting non-proper Euler parameters back on the unit circle, we chose to implement a mixed approach for integrating equation (1.29). The body angular velocity from equation (1.31) , integrated with MK4, it is used within a RKMK4 algorithm described in [13] , This second integration scheme preserves the unit of the quaternion via the exponential update^ and although computationally more intensive, it provides mathematically accurate Euler parameters without resorting to patch normalization or relying on a small time steps to hide the numerical error therein.

Two main operators are needed to implement the RKMK4 algorithm: the exponential mapping of a Cartesian vector to quaternion and the inverse right Jacobian associated with the body angular velocity - see [13] and [16], The exponential mapping of a vector u is expressed as follows: while the inverse right Jacobian associated with the body angular velocity of a vector u is expressed as follows:

(3-4)

Once we have these operators, we can implement the mixed RK4-RKMK4 approach. We begin integrating equations (1.28), (1.30) and (1.31) for the i-th stage using the current orientation quaternion - either IC value or previous converged value, using the equations set (3.1) for stage i-th. At this point, we compute the corresponding i-th stage for the quaternion using the RKMK4 outlined in detail below. This will use the just calculated i-th stage value of body angular velocity to compute the corresponding i-th stage value of the quaternion for equation (1.29). The process then repeats for the i-th + 1 stage, by using first RK4 - for global position, body linear velocity and body angular velocity - and then RKMK4, for quaternion orientation.

Equation (3.1) is well known in literature, thus we leave to the reader its expansion for the stages 1, The corresponding equation set for RKMK4 - Runge-K utta-M unth^Kaas is detailed as follows, using the exponential mapping (3.3) and the inverse right Jacobian (3.4):

We can now expand equation (3.5) for each of the four stages using the corresponding computed angular velocity from (3.1) - the formulas show the moment when we must send information to RK4 (#2 providing the current orientation) and when we must receive information from RK4 (#3 retrieving the corresponding body angular velocity). This synchronization is important to maintain state consistency between the two algorithms working in parallel.

• The first stage of RKMK4 for s = 1 is:

• The second stage of RKMK4 for s = 2 is:

• The third stage of RKMK4 for s — 3 is:

• The fourth stage of RKMK4 for s = 4 is:

• The summation stage or RKMK4 which provides the dt quaternion update is:

Finally, we have all the detailed stage expansions for integrating the quaternion from value θ _k to value for a time step dt, as well as the parallel integration of equations (1.28), (4.30) and (1.31), based on the i-th stage quaternion value. As is evident, no quaternion normalization is required, and it can be verified by numerical assertions that the above operations maintain quaternion state integrity to correctly represent a rigid body rotation.

3.2 Specialized Containers for Mixed Approach

The above exponential integration requires more floating point operations than the standard Euler approach (RK4). Thus, specialized containers for high performance matrix computation should be developed by the implementation of this code to maintain real-time performance even with higher computational costs.

From our testing, a vectorialized approach in 3x3 matrix multiplication using a modern compiler was able to achieve high computational speed. Given the sequential requirements of the sets of operations between RK4 and RKMK4, parallel implementation did not provide any performance improvement. Dynamic memory also showed a consistent slowdown in computation.

Compiler-based vectorization from Visual Studio 2019 provides about 10% speed improvement in a standard N ² loop needed for 3x3 matrix multiplication with respect to the open source Eigen library for the standard RK4 algorithm.

By using these containers over the RKMK4 algorithm, we measured a 5% slowdown with respect to the Eigen library for the RKMK4. This delta of 15% gap between RK4 and RKMK4 is consistent with experiments shown in literature - see [13]. However, the higher quality of results on a fast package composed by RTS RigidBody, NDI Controller, and Mixed RK4/RKMK4 integrator, is the preferred choice for a real-time physics engine.

Figure 14: Operators for Hamiltonian and Lie Group algebras.

From a mathematical standpoint, the above-mentioned specialized containers must implement both the Hamiltonian algebra to efficiently handle quaternion operations, see (3.2), and the Lie Group algebra to handle the RKMK4 operations, see (3.3) and (3,4) . The reader can review the mathematical details of both algebras in [13] and in [16]. Figure 14 shows the operators which implement those two algebras, while figure 15 shows the class declarations for the components described in the previous sections.

Figure 15: Class declarations for physics engine components.

4 Airbus A320 Aircraft Input Data

In this section, we will describe the parameterization of both rigid body and controller to simulate Airbus A320 commercial aircraft. The major differentiator from existing literature modeling is the “data-driven” design approach used for the controller. A set of NDI controllers such as the ones described in the previous section are generally model dependent, in the sense that properties of the real model which is simulated need to be carefully replicated in the controller to provide accurate forces to the rigid body. This is the approach followed by [10], [7], [8], [6]. Although this approach is widely adopted in the industry, it clashes with the “data-driven” approach which is a requirement of this physics engine, particularly in the modeling of the “control mixer/surface deflection” componentization. For example, an Airbus A320 has wings with six ailerons, one rudder and one stabilizer. A Lockheed Martin F-117 Nighthawk has a completely different wing profile. Therefore, a control mixer developed for the Airbus would need to be re-written for the Nighthawk, providing a “model-driven” design.

4.1 Integral Effect of Force/Moments

Figure 16 based on the design suggested in [8] shows how our novel approach differs from existing literature:

Figure 16: Implemented “data-driven” approach in NDI controller.

Both Engine Dynamics and Actuators (control mixer) have been replaced with the body forces written in the local rigid body reference frame. An integral effect of the three forces (2.10) and three moments (2.20) has been computed for the rigid body local frame, and inertial correction considering the fictitious forces has been added to the Newton/Euler equations of motion (1.23).

The customer data which are necessary to compute these forces, as well as the rigid body parameters needed to characterize the rigid body, can be divided into four sections: 1. Inertial - used by RigidBody and Controller:

• mass of the aircraft

• moment of inertia of the aircraft

2. Geometrical - used by Controller:

• surface reference for the airfoil

3. Aerodynamics - used by Controller:

• airfoil lift coefficient versus angle of attack

• polar diagram lift versus drag coefficient

4. Performance - used by Control Law:

• climb and descent rate for the aircraft

• heading change rate for the aircraft

Let us examine how these data are used with reference to the equations developed in the previous sections. The mass and the moment of inertia of the vehicle are used both in the differential equations for rigid body - Newton/Euler ODEs (1.30) and (1.31) - as well in the NDI controller #1 and #3 - equations (2.7), (2.9), (2.10), (2.20). The moment of inertia needs to be computed in RTS from Tait-Bryan using equation (1.22).

Equation (2.8) uses this parameter within the lift and drag coefficients and the dynamic pressure

Q, defined as:

(4-1)

The density of the air must be computed according to the geopotential altitude of the aircraft along with the speed of sound to compute the Mach number (speed) of the aircraft, according to the following look up table in figure 17:

Figure 17: Density and speed of sound as function of the geopotential altitude.

The airfoil lift coefficient versus angle of attack has already been discussed in equation (2.13) and relative figure (5) .

The polar diagram lift versus drag coefficient is another key parameter which ties together the lift and drag depending on the profile of the airfoil. These coefficients are necessary for computing equation (2.8) and a typical customer example can be seen in figure 18.

Figure 18: Lift versus Drag polar coefficients as function of multiple angles of attack.

Climb rate (2000 feet per minute for A320) and heading rates (120° in 600 s for A320) are used to compute the correct control laws as explained in the co- simulation section. Their contribution is used in equations (2.21), (2.22), and (2.23) and their examples are used in figures 11 and 12. For the velocity control lave, we tune a change in velocity with the same maneuver time of the climb rate in order to maintain an angle of attack of 1° which enables a steady flight at the beginning and ending of maneuver - the relevant control law can be seen in figure 10.

These steady flight conditions can be calculated by means of a performance spreadsheet with the above information. Below are reported two calculation examples, with the performance parameters necessary to maintain 1° angle of attack at different altitudes, for Airbus A320 (figure 19) and Lockheed F-117 Nighthawk (figure 20). It is noteworthy that, in order to climb from 5,000 m to 10,000 m, the commercial aircraft takes 492 s versus the 164 s of the military aircraft.

5 Accuracy Validation of Dynamic Response

In this section, we will provide a comparison of the package developed in C++ through the imple- mentation of the formulations derived in this paper versus Simulink® Aerospace Blockset™. This commercial package provides a similar 6DOFs rigid body system implemented in Tait-Bryan for- mulation and uses Euler angles as orientation parameters. The Simulink model has been simulated with a standard RK4 solver - named ode4 - and fixed time step, in order to provide similar inte- gration parameters to the C++ physics engine. This same product has also been used in [10] et al. to validate typical NDI controllers with “control mixer/surface deflection”.

In addition to the standard Simulink 6DOFs block, we implemented in the same block diagram a version of our controller using MATLAB® Functions™. This was necessary to provide a comparison with our C++ physics engine comprised of the classes shown in figure 15.

It also should be noted that Simulink® Aerospace Blockset™, along with other MathWorks products, are the gold standard in the aerospace industry for validating real aircraft dynamics - see Lockheed Martin testimonial in figure 21. Thus, validating the output of the new C++ physics engine against Simulink® Aerospace Blockset™, we confirm the quality of the results of our physics engine.

Figure 21: Lockheed Martin uses MathWorks products to predict F-35 fleet performance.

5.1 In-core Real Time Simulation

The following are the key differentiators between the new C++ physics engine and the Simulink® Aerospace Blockset :

• We built a model in Simulink Aerospace Blockset using the Euler angles formulation for the 6DOFs model - this was done on purpose to verify our implementation with quaternion orientation used within the C++ Rigid Body.

• The mixed algorithm RK4/RKMK4 differs from the standard ode4 implemented in Simulink, which is similar to a fixed time step RK4.

• Simulink works on Tait-Bryan reference system formulation, while our package has been entirely written in RTS formulation.

• Customers use Simulink models to run offline simulations and then post-process results. We instead use a game engine tick to drive the rigid body in a real-time simulation.

The comparative model in Simulink is shown in figure 22.

Figure 22: Simulink® Aerospace Blockset™ with MATLAB® Functions™.

Comparing figure 22 with our C++ class implementation shown in figure 15, we can make the following functional correlation:

• “Simulink 6DOF (Euler Angles)” corresponds to the classes Matrix4D, RigidBodyD, Vec3D, Mat 3D, and QuatD

• “Al Control Laws” corresponds to the class ControlLawD • “NDI Controller” corresponds to the classes ControllerBaseD, AirplaneControllerD, and

EulerAnglesD

The block masks “Al Controller Laws” and “NDI Controller” contain a set of MATLAB Func- tions used to implement in m-code the control laws and the NDI controller logic implemented in C++ classes to allow a one-to-one comparison between an offline simulation over Simulink model and an in-core simulation of the physics engine.

Figure 23 and figure 24 show the MATLAB function blocks containing the source code for the NDI controller and the control laws. It should be clear why we need a conversion in/out to the NDI RTS controller block. Simulink 6DOFs rigid body accepts only states conditions in Tait-Bryan formulation, hence for the sake of this comparison, a conversion TB → RTS must be the input of this RTS controller, and a conversion RTS → TB must be the output of this RTS controller. This guarantees state consistency within the Simulink platform - see equations (1.14) through (1.22) for details.

Figure 23: NDI controller implemented with MATLAB Functions.

Figure 24: Control laws implemented with MATLAB Functions.

This model setup can now be used for an accurate comparison against the C++ physics engine.

5.2 Constant Accuracy of l.e-6 (SI Physics Units)

We can now carry out a detailed state comparison at each time step between the Simulink output and the C++ physics engine. Given a fixed time step equal to 0.1 s with ode4 and a simulation time of 2,500 s in Simulink, we can compare the error after 25,000 solver evaluations, each evaluation composed by four stages. It is clear that any error in the integration state would accumulate over time, hence claiming an accuracy over 1 s is very easy, while it is much harder to obtain a constant accuracy after 41.67 min of flying time, for example.

It should be evident that by evaluating such a long simulation time, we are testing a potential high level of numerical inaccuracy - a small error would cumulatively sum-up over time, providing a divergence of results at the end of the simulation. Thus, it is numerically very challenging to achieve such level of accuracy by two different numerical solvers which integrate differential equations produced by such different formulations - see the formulas of this paper for reference.

Thus, our achieved results of a constant accuracy across all the states of the rigid body - global position, global linear velocity, body linear velocity, body angular velocity, Euler angles - guarantees without a doubt that the output of the two software are consistent without a divergence in integration error. The dynamics of two aircraft simulated by different software - Simulink versus C++ physics engine - is exactly the same. With that in mind, we can claim an accuracy of l.e-6 in SI Physics Units over a traditional offline simulation, while implementing a real-time in-core simulation driven by the game tick.

The above claims have been proven by internal unit tests over the developed libraries and can be seen on the following series of plots, where only a zoom below l.e-6 reveals a differentiation in the floating point values - see figure 25, 26, and 27. Note also that, in figure 27, the derivative of the traces of the error are parallel, which proves a constant offset below l.e-6 with no accumulation of error through the dynamic simulation. This means that the two software produce outputs which are consistent at the micron level (1 pm). Moreover, as seen in figure 25, the C++ physics engine is able to simulate 41.67 min in 0.38 s.

Figure 25: 2,500 s comparative simulation with Simulink Aerospace Blockset.

Figure 26: Details of second order dynamics response.

Figure 27: Parallel traces between Simulink and Rigid Body, differing at l.e-6.

For completeness, we will conclude by showing the full results of the simulation for the C++ physics engine during the test maneuver: climb from 5,000 m to 10,000 m, change in heading of 120°, descend to 5,000 m. The internal states for the C++ physics engine have been exported in mat file via MathWorks MATFile facilities, to allow post-processing. The plots show the output in standard aerospace NED reference frame.

Figure 28: Trajectory summary for 2,500 s simulation on Airbus A320.

Figure 29: Global position - C++ physics engine transformed into NED.

Figure 30: Global linear velocity - C++ physics engine transformed into NED.

Figure 31: Body linear velocity - C++ physics engine transformed into NED.

Figure 32: Body angular velocity - C++ physics engine transformed into NED.

Figure 33: Euler angles orientation - C++ physics engine transformed into NED.

Figure 34: Forces computed by controller - C++ physics engine transformed into NED.

Figure 35: Torques computed by controller - C++ physics engine transformed into NED.

Figure 36: Angle of attack and side-slip - C++ physics engine transformed into NED.

References

[1] A. A. Shabana, Dynamics of Multibody Systems, Cambridge University Press (2013).

[2] NASA Mission Planning and Analysis Division, Euler Angles, Quaternions, and Transformation Matrices, NASA Lyndon B. Johnson Space Center (Jul. 1977).

[3] W. G. Breckenridge, Quaternions proposed standard conventions, NASA Jet Propulsion Lab- oratory, Technical Report (Oct. 1979).

[4] R. F. Stengel, Flight Dynamics, Princeton University Press (2004).

[5] M. Zyskowski, Aircraft Simulation Techniques Used in Low Cost, Commercial Software, AIAA Modeling and Simulation Technologies Conference and Exhibit (2003).

[6] D. Enns et al., Dynamic inversion: An evolving methodology for flight control design Interna- tional Journal of Control, vol. 59, no. 1, pp. 71-91 (1994).

[7] R. Adams et al., Design of nonlinear control laws for high angle of attack flight, Journal of Guidance, Control, and Dynamics, vol. 17, no. 4, pp. 737-746 (1994).

[8] S. A. Snell et al., Nonlinear Inversion Flight Control for a Supermaneuverable Aircraft, Journal of Guidance, Control, and Dynamics, vol. 15, no. 4, pp. 976-984 (1992).

[9] J. Beaty, Introduction to Digital Signal Processing (DSP), Systems Analysis & Concepts Di- rectorate, NASA Langley Research Center, Hampton, VA (2011).

[10] D. Dierker, An Automated Controller Design Methodology for Six Degree of Freedom Aircraft Models, Master’s thesis, Wright State University, Dayton, Ohio (2017).

[11] K. Shimmin, An architecture for rapid modeling and simulation of an air vehicle system, Master’s thesis, University of Dayton, Dayton, Ohio (2016).

[12] J. Brendlinger, Development of Guidance Laws for a Reduced Order Dynamic Aircraft Model, Master’s thesis, Wright State University, Dayton, Ohio (2017).

[13] A. Sveier, et al., Applied Runge-Kutta-Munthe-Kaas Integration for the Quaternion Kinematics, Journal of Guidance, Control, and Dynamics, vol. 42, no. 12, pp. 2747-2754 (2019).

[14] M. S. Andrle, and J. L. Crassidis, Geometric Integration of Quaternions, Journal of Guidance, Control, and Dynamics, Vol. 36, No. 6, pp. 1762-1767 (2013).

[15] P. E. Crouch, and R. Grossman, Numerical Integration of Ordinary Differential Equations on Manifolds, Journal of Nonlinear Science, Vol. 3, No. 1, pp. 1-33 (1993).

[16] A. Iserles, H. Munthe-Kaas, S. Nprsett, and A. Zanna, Lie-Group Methods, Acta Numerica, Vol. 9, pp. 215-365 (2000).