Flight dynamic model

FMS is one of the freeware high-quality flight simulator for model airplanes. It is written in Delphi language with some third party components. You can download it from http://n.ethz.ch/student/mmoeller/fms/index_e.html
My question is about the "flight dynamic model" of FMS and/or other similar simulators. Do you know any Delphi and/or C/C++ flight dynamics model code which takes the rudder, elevator, throttle, aileron, airplane's position, velocity and attitude (roll, pitch, yaw angles) as inputs and calculates the airplane's next position, velocity and attitude?
Flight Dynamics Model : =================INPUTS: Time increment : dT (delta T in milli second) Control inputs : Rudder Elevator Throttle Aileron Airplane's current state : Position (in 3D. ie X,Y,Z) Velocity (in 3D, Vx, Vy,Vz) Attitude (in 3D, Roll, Pitch, Yaw angles) Airplane's constants: Wing span Weight ???? ???? etc.. OUTPUT: Airplane's next state after dT : Position (in 3D. ie X,Y,Z) Velocity (in 3D, Vx, Vy,Vz) Attitude (in 3D, Roll, Pitch, Yaw angles)

| My question is about the "flight dynamic model" of FMS and/or other similar | simulators. | Do you know any Delphi and/or C/C++ flight dynamics model code which takes | the rudder, elevator, throttle, aileron, airplane's position, velocity and | attitude (roll, pitch, yaw angles) as inputs and calculates the airplane's | next position, velocity and attitude?
There's a few open source flight simulators out there that obviously must do what you're referring to.
They include --
Slope Soaring Simulator -- http://www.rowlhouse.co.uk/sss / CRRCsim -- http://groups.yahoo.com/group/crrcsim / Flight Gear Flight Simulator -- http://seneca.flightgear.org/flightgear/links.html
And there's more, just google for 'em.
This mention on the FGFS page looks like it might be exactly what you want --
LaRCsim - A set of ANSI C routines that implement a full set of equations of motion for a rigid-body aircraft in atmospheric and low-earth orbital flight, suitable for pilot-in-the-loop simulations on a workstation-class computer.

Hi Lora, Inspired by FMS I started writing a R/C simulator last year. Never got around to finishing it though. But I figured out a working flightmodel. Mostly by directly implementing some basic physics. Of course I haven't documented my code well :-) but maybe it's a start for you. If I remember correctly it worked pretty well, but some tweaking is needed when pitch approaches -90 degrees. Sorry about the long post.
Good luck!! Davy
Some declarations :
struct AirplaneParameters { Float wingspan; Float wingsurface; Float Cl0; // lift at alpha=0 Float Cla; // lift coefficient (=lift curve slope) Float Cd0; // drag at alpha=0 Float Cda; // drag coefficient (=drag curve slope) Float Clb; // lift coefficient of fuselage (f(beta)) Float Cdb; // drag coefficient of fuselage (f(beta)) Float Cdx; // drag coefficient x-axis Float Cdy; // drag coefficient y-axis Float Cdz; // drag coefficient z-axis Float maxThrust; // maximum thrust Float thrustAngle; // angle of thrust vector and fuselage. Float weight; // weight Float aileron; // aileron effect Float rudder; // rudder effect Float elevator; // elevator effect Float yawDamp; // yaw damping Float pitchDamp; // pitch damping Float rollDamp; // roll damping };
class FlightModel { /* .... */
private: AirplaneParameters m_par; // Position Float m_dX; Float m_dY; Float m_dZ; // Speed Float m_U; Float m_V; Float m_W; // Acceleration Float m_aU; Float m_aV; Float m_aW; // Attitude Float m_phi; // roll Float m_theta; // pitch Float m_psi; // yaw // Worldattitude Float m_wphi; // roll Float m_wtheta; // pitch Float m_wpsi; // yaw // Angular velocity Float m_P; Float m_Q; Float m_R; // Angular acceleration Float m_aP; Float m_aQ; Float m_aR; // Moments Float m_L; Float m_M; Float m_N; // Forces Float m_Fx; Float m_Fy; Float m_Fz; // Cumulated movements Float m_dXcumul; Float m_dYcumul; Float m_dZcumul; DirectX::Vector3 m_velocity;
private: Float m_speedsq; Float m_alpha; Float m_beta; Float m_thrust; Int m_throttle; /* 0..100 */ Int m_yokeHor; /*-100...+100*/ Int m_yokeVer; /*-100...+100*/ Int m_pedals; /*-100...+100*/ };
Implementation of the main function (needs to be run about every 5 ms). Some debug code is still present in there
const float _rho = 10.0f; const float _grav = 29.0f;
DWORD WINAPI runFlightModel(void* flightModel) { FlightModel* fm = reinterpret_cast<FlightModel*>(flightModel); DirectX::Timer timer; timer.microElapsed( ); CRITICAL_SECTION cs; ::InitializeCriticalSection(&cs);
while (fm->running( )) { ::EnterCriticalSection(&cs ); Huge helapsed = timer.microElapsed( ); Float elapsed = helapsed / 1000000.0f; if (elapsed > 0.5f) elapsed = 0.0f; fm->updateInput( ); // initialize some basic numbers fm->m_speedsq = fm->m_U*fm->m_U + fm->m_V*fm->m_V + fm->m_W*fm->m_W; Float speed = (Float)(sqrt(fm->m_speedsq)); if (absval<Float>(fm->m_U) < 0.01f) { fm->m_alpha = 0; //DirectX::Pi / 2; fm->m_beta = 0; //DirectX::Pi / 2; } else { fm->m_alpha = atanf(fm->m_W/fm->m_U); fm->m_beta = atanf(fm->m_V/fm->m_U); }
fm->m_thrust = fm->m_par.maxThrust*0.1f*(fm->m_throttle);
// calculate lift and drag Float Dragx = fm->m_par.Cdx*fm->m_U*fm->m_U*sign<Float>(fm->m_U); Float Dragy = fm->m_par.Cdy*fm->m_V*fm->m_V*sign<Float>(fm->m_V); Float Dragz = fm->m_par.Cdz*fm->m_W*fm->m_W*sign<Float>(fm->m_W); Float Lift = (fm->m_par.Cl0 + fm->m_par.Cla*fm->m_alpha)*fm->m_U*fm->m_U;
if (absval<Float>(fm->m_alpha) > 0.2f) Lift = Lift - Lift*(sinf(fm->m_alpha - 0.2f)); if (fm->m_U < 0.0f) Lift = fm->m_par.Cl0;
Float FLift = (fm->m_par.Clb*fm->m_beta)*fm->m_U*fm->m_U;
fm->m_Fx = -_grav *sinf(fm->m_theta)*fm->m_par.weight - Dragx + fm->m_thrust*cosf(fm->m_par.thrustAngle); fm->m_Fy _grav*sinf(fm->m_phi)*cosf(fm->m_theta)*fm->m_par.weight - Dragy - FLift; fm->m_Fz _grav*cosf(fm->m_phi)*cosf(fm->m_theta)*fm->m_par.weight - Dragz - Lift - fm->m_thrust*sinf(fm->m_par.thrustAngle);
fm->m_aU = (fm->m_V*fm->m_R - fm->m_W*fm->m_Q)/1.0f + fm->m_Fx/(fm->m_par.weight); fm->m_aV = (fm->m_W*fm->m_P - fm->m_U*fm->m_R)/1.0f + fm->m_Fy/(fm->m_par.weight); fm->m_aW = (fm->m_U*fm->m_Q - fm->m_V*fm->m_P)/1.0f + fm->m_Fz/(fm->m_par.weight);
if (absval<Float>(fm->m_aU) > 10000) { fm->m_aU = 0.0f; fm->m_U 0.0f; } if (absval<Float>(fm->m_aV) > 10000) { fm->m_aV = 0.0f; fm->m_V 0.0f; } if (absval<Float>(fm->m_aW) > 10000) { fm->m_aW = 0.0f; fm->m_W 0.0f; }
// calculate new speeds fm->m_U += elapsed*fm->m_aU; fm->m_V += elapsed*fm->m_aV; fm->m_W += elapsed*fm->m_aW;
// calculate new position fm->m_dX = elapsed*fm->m_U; fm->m_dY = elapsed*fm->m_V; fm->m_dZ = elapsed*fm->m_W;
Float ailDefl = fm->m_yokeHor / 200.0f; Float elDefl = fm->m_yokeVer / 300.0f; Float rdDefl = fm->m_pedals / 300.0f; fm->m_L = (fm->m_par.rollDamp*fm->m_P/(fm->m_U + 0.1f) + fm->m_par.aileron*ailDefl)*fm->m_U*fm->m_U; fm->m_M = (fm->m_par.pitchDamp*fm->m_Q/(fm->m_U + 0.1f) + fm->m_par.elevator*elDefl)*fm->m_U*fm->m_U - 0.2f*fm->m_alpha*sign<Float>(fm->m_U)*fm->m_speedsq; fm->m_N = (fm->m_par.yawDamp*fm->m_R/(fm->m_U + 0.1f) + 0.2f*fm->m_beta*sign<Float>(fm->m_U) + fm->m_par.rudder*rdDefl)*fm->m_U*fm->m_U;
fm->m_aP = fm->m_L*elapsed; fm->m_aQ = fm->m_M*elapsed; fm->m_aR = fm->m_N*elapsed;
if (absval<Float>(fm->m_aP) > 1000) { fm->m_aP = 0.0f; fm->m_P 0.0f; } if (absval<Float>(fm->m_aQ) > 1000) { fm->m_aQ = 0.0f; fm->m_Q 0.0f; } if (absval<Float>(fm->m_aR) > 1000) { fm->m_aR = 0.0f; fm->m_R 0.0f; }
fm->m_P += elapsed*fm->m_aP; fm->m_Q += elapsed*fm->m_aQ;// - sign<Float>(fm->m_U)*fm->m_alpha*(0.0003f/(elapsed*elapsed)); fm->m_R += elapsed*fm->m_aR;// - fm->m_beta*(0.0003f/(elapsed*elapsed));
Float phi = fm->m_phi; Float theta = fm->m_theta; Float psi = fm->m_psi;
if (absval<Float>(cosf(theta)) < 0.00001f) if (cosf(theta) < 0.0f) { theta = acosf(-0.00001f); fm->m_theta = theta; } else { theta = acosf(0.00001f); fm->m_theta = theta; }
fm->m_phi += (fm->m_P + (fm->m_Q*sinf(phi) + fm->m_R*cosf(phi))*tanf(theta))*elapsed; fm->m_theta += (fm->m_Q*cosf(phi) - fm->m_R*sinf(phi))*elapsed; fm->m_psi += ((fm->m_Q*sinf(phi) + fm->m_R*cosf(phi))/cosf(theta))*elapsed;
fm->m_phi = modulo<Float>(fm->m_phi, 2*DirectX::Pi); fm->m_theta = modulo<Float>(fm->m_theta,2*DirectX::Pi); fm->m_psi = modulo<Float>(fm->m_psi, 2*DirectX::Pi);
// update cumulative movements Float X = 0.01f*fm->m_dX, Y = 0.01f*fm->m_dY, Z = 0.01f*fm->m_dZ; Float dX = cosf(fm->m_theta)*sinf(fm->m_psi)*X + (sinf(fm->m_phi)*sinf(fm->m_theta)*sinf(fm->m_psi) - cosf(fm->m_phi)*cosf(fm->m_psi))*Y + (cosf(fm->m_phi)*sinf(fm->m_theta)*sinf(fm->m_psi) - sinf(fm->m_phi)*cosf(fm->m_psi))*Z; Float dY = -sinf(fm->m_theta)*X + sinf(fm->m_phi)*cosf(fm->m_theta)*Y + cosf(fm->m_phi)*cosf(fm->m_theta)*Z; Float dZ = cosf(fm->m_theta)*cosf(fm->m_psi)*X + (sinf(fm->m_phi)*sinf(fm->m_theta)*cosf(fm->m_psi) - cosf(fm->m_phi)*sinf(fm->m_psi))*Y + (cosf(fm->m_phi)*sinf(fm->m_theta)*cosf(fm->m_psi) + sinf(fm->m_phi)*sinf(fm->m_psi))*Z;
fm->m_dXcumul -= dX; fm->m_dYcumul -= dY; fm->m_dZcumul -= dZ;
fm->m_velocity DirectX::Vector3(-dX/elapsed, -dY/elapsed, -dZ/elapsed);
::LeaveCriticalSection(&cs); ::Sleep(5);
} ::DeleteCriticalSection(&cs); return S_OK; }
Some example parameters for an aircraft :
AirplaneParameters _parsExtra300 { 1.4f, 0.4f, 0.0003f, // lift at alpha=0 0.0180f, // lift coefficient (=lift curve slope) 0.02f, // drag at alpha=0 0.04f, // drag coefficient (=drag curve slope) 0.0180f, // lift coefficient of fuselage (f(beta)) 0.10f, // drag coefficient of fuselage (f(beta)) 0.002f, // drag coefficient x-axis 0.002f, // drag coefficient y-axis 0.050f, // drag coefficient z-axis 10.0f, // maximum thrust 0.0175f, // thrustangle 2.2f, // weight 0.150f, // aileron effect 0.150f, // rudder effect 0.650f, // elevator effect -4.0f, // yaw damping -4.0f, // pitch damping -2.0f // roll damping };
  Quoted Text Here  
http://n.ethz.ch/student/mmoeller/fms/index_e.html
similar

Dear Davy,
Thank you very much for sharing your experience and code with me/us. I'm sure I will learn a lot from your experience.
THANK YOU
Lora Lee
  Quoted Text Here  

Lora Lee wrote:
  Quoted Text Here  
Very cool, gang, I've never done anything in Delphi, only in Fortran, but it looks good. Looks like you stopped short of making the graphical part (am I right here or not?). Seems like that would be a challenge, but maybe not as hard as I think for someone familiar with it. It looks like you're applying all the forces to the plane as a whole, and airflow is regarded as uniform over the whole plane. Note that without using a velocity-panel type model, where airflow is calculated individually at each location, that you will not be able to do any non-linear maneuvers, such as spins, snap rolls,lomcevaks, etc. but then again, this is how many commercial flight sims are made. Good job.     Paul

Hi Paul and Lora, For the record, the code is in c++. You're correct : the flight model is simple and limited, but it was the best I could do in a couple of weeks without taking an advanced course in aerodynamics :-) About the graphical side : I did some work on it (using DirectX). It actually had animated flight surfaces, smoke, shadows, sun flare, engine noise with Doppler effect... I have put some screenshots up on http://users.pandora.be/playinginthedark/airdancers/airdancers1.jpg
http://users.pandora.be/playinginthedark/airdancers/airdancers2.jpg
It was about a month and a half of work to get to this point. I never finished it because - well - the fun part was done and there was still A LOT of work left : collision model, scenery modeling, configurable controls, network code, ... Lora, can you let me know if you ever find a use for the code? Thanks!
Davy
  Quoted Text Here  
remember

Davy,
Do you also some documentation explaining the math that you used in your program?
R.R.
  Quoted Text Here  

Well, it's been a while, but I got a lot out of the first part of the following paper: http://www.movesinstitute.org/~zyda/pubs/Presence.1.4.pdf
There are two parts in the calculations : position and attitude. Basically you calculate the forces on the aircraft (gravity, drag, lift and thrust), say F for each of the aircraft axes. Then you calculate the acceleration using Newton F=m*a or a = F/m. This will give you the new position. The hard part is the rotations. Looking back at the code, I honestly have no idea where I got these but I remember it took me a couple of weeks to come to these. I can't guarantee the physical soundness, but the result was acceptable to me. The last calculations of dX, dY and dZ are simply translation from the aircraft axis to the world axis (some matrix math). There is also still some debug code in there (the "if ( ) > 50") because in the beginning I had some asymptotic behavior in some cases. Also, never forget to take the elapsed time since the previous run into account in the right order and places.
It's not simple and my implementation is far from perfect or finished, but it was great fun figuring it out!
Davy
  Quoted Text Here  
remember