Working on small flight sim game, need help with equations

124

Author

March 17, 2017 09:53 PM

Hi, just registered because have a problem with my hobby flight-sim game project (written in JavaScript). I'm not sure if the equations in update() method are correct or more important - their order (which one should be first, second, etc) and if I messed something in using metric system (unit conversion). My doubts are caused by that the bot is somehow not capable of going through test map's checkpoints, it crashes at second checkpoint (for now I'm drawing in canvas with 512x512 px dimensions)

The initial values of bot's characteristis are:

100 dragpoints, 15000 kgs thrust, 5000kgs mass, 5 meters height, 15 meters wingspan, 70 sq. meters wingarea and 1000 liters of fuel

This is the sim properties and bot's update() method, I'm starting with fuel level check and continue to expand to calculations needed for setting current position and altitude. It is obvious that sim variables are const (not UPPERCASEd for now)

---


var PI = 3.141592;
var TAU = 2 * PI;


var sim = {
  state: 1, 
  startTime: 0, playedTime: 0, tempPlayedTime: 0, currentTime: 0, lastTime: 0, dTime: 0, timeProduct: 0, updateTime: 0, worldUpdates: 0, 
  distPerFrame: 0, 
  timeScale: 1, /*< 1 speeds up, > 1 slows down the game, must be fixed*/ 
  timeStep: 0.01667, //in s
  timeDiff: 0, temp: 0, 
  referenceGravity: 9.80665, //reference gravity in m/sq. s at sea level
  currentGravity: this.referenceGravity, //in m/sq. s
  earthRadius: 6371000, //in m
  earthPowRadius: 40589641000000, //in m
  geopotentialHeight: 0, //in m
  distance: 0, x: 0, y: 0, 
  racetrackLength: 0, //in m
  airReferenceTemperature: 15, //reference air temp in Celsius at sea level, or 288.15 K 
  airReferenceViscosity: 0.01827, //0.00001789 reference air viscosity in centipoise at reference temperature in Celsius, or 0.01827 at reference temperature in Rankine
  airReferencePressure: 101325, //reference air pressure in Pa at sea level and 15'C air temperature
  airReferenceDensity: 1.225, //reference air density at sea level in kg/cub.m, or 0.0764734 lb/cub. ft
  RankineRefTemp: 518.76, //reference air temp in Rankine degree == 15'C
  SutherlandConstant: 120, //for air
  gasConstant: 8.3144598, //reference gas constant at sea level, in J/mol Kelvin, or 1545.35 ft lb/lbmol Rankine
  airMolarMass: 0.0289644, // in kg/mol
  airLapseRate: 0.0065, //0.0098 reference air lapse rate in Celsius per meter altitude, 0.0065 reference air lapse rate in Kelvin per meter, 0.0019812 in Kelvin per foot
  airViscosity: this.airReferenceViscosity, 
  airPressure: this.airReferencePressure, 
  airDensity: this.airReferenceDensity, 
  kinematicViscosity: 0, 
  ReynoldsNumber: 0, 
  airCurrentTemp: this.airReferenceTemperature, 
  counter: 0, 
  showFPSMem: true, currentFPS: 0, currentMS: 0, currentMem: 0, showTraj: false, showUnitStats: true, 
  botIdCounter: 0, bots: [], runningBots: 0, 
  targetSize: 10
};


var bot = function() {
  this.position = {x: 0, y: 0, z: 0}; //z is the ALTITUDE!!!
  this.destination = {x: 0, y: 0, z: 0}; //z is the ALTITUDE!!!
  this.distance = 0; //in m
  this.distanceTravelled = 0; //in m
  this.altitude = this.position.z; //in m
  this.pitch = 0; //angle of attack ? in degrees?
  this.roll = 0; //in degrees?
  this.yaw = 0; //in degrees?
  this.heading = 0; //in degrees, 0 - East
  this.dragCoeff = dragpoints * 0.0001;
  this.thrust = thrust; //in kgf
  this.fuel = fuel; //in liters
  this.mass = mass; //in kg
  this.loadedMass = this.fuel + this.mass; //in kg
  this.height = height; //in m
  this.acceleration = this.thrust / this.loadedMass; //acceleration in m/s
  this.oldAcceleration = 0;
  this.velocity = {x: 0, y: 0, z: 0}; //maybe needed for proper integration
  this.speed = 0; //current speed in m/s
  this.wingSpan = wingspan; //in m
  this.wingArea = wingarea; //in sq.m
  this.wingLoading = this.loadedMass / this.wingArea; //in kg/sq.m
  this.wingType = wingtype; //1 - flat bottom, 2 - semi-symmetric, 3 - symmetric
  this.frontalArea = this.wingSpan * this.height; //in sq.m
  this.liftCoeff = 0;
  this.liftForce = 0; //in kN
  this.fuelConsumption = (this.thrust / 100) * 0.035; //0.035 litres per second
  this.trajectoryLine = traj;
  this.rotationSpeedX = TAU / 100; //100 FPS
  this.rotationSpeedY = TAU / 100;
  this.rotationSpeedZ = TAU / 100;
  this.size = this.wingArea / 10; //size of the bot on the canvas
  this.id = sim.botIdCounter;
  this.target = 0;
  this.hasTarget = false;
  this.time = 0;
  this.startTime = 0;
  this.isTimed = false;
  this.state = true;
  this.distanceBeforeFirstTarget = 0;
  this.isLanding = false;
  sim.botIdCounter++;
  sim.runningBots++;

  this.update = function() {
    if (this.fuel >= this.fuelConsumption * sim.timeStep) { 
       
  this.fuel -= this.fuelConsumption * sim.timeStep;
        this.loadedMass -= this.fuelConsumption * sim.timeStep;
        this.acceleration = (this.thrust*sim.timeStep) / this.loadedMass;
      } else {
        if (this.target <= targetCount) {
          this.fuel = 0; //directly set fuel to 0 when it is < 0 or has some very low values ~0.050l
          if (this.acceleration > 0) this.acceleration -= (this.thrust*sim.timeStep) / this.loadedMass;
          if (this.acceleration < 0) this.acceleration = 0;
        } else {
            this.fuel = 0;
            this.isTimed = false; //stop the stopwatch for our bot
            this.state = false; //set state of our bot to inactive
            sim.runningBots--; //remove our bot from the list with active bots
            return; //maybe it's better to early exit from the function, execution of the conditions below is useless
        }
      }
      
      //after fuel check we can calculate all characteristics of our bot
      this.wingLoading = this.loadedMass / this.wingArea;
      this.liftCoeff = ((this.loadedMass * sim.currentGravity) / ((sim.airDensity / 2) * pow(this.acceleration, 2) * this.wingArea)) * sim.timeStep;
      this.liftCoeff = isNaN(this.liftCoeff) ? 0 : this.liftCoeff;
      this.liftForce = (this.liftCoeff * (0.5*(sim.airDensity * pow(this.acceleration, 2))) * this.wingArea) * sim.timeStep;
      this.liftForce = isNaN(this.liftForce) ? 0 : this.liftForce;
    
      //set destination coords for our bot only once for each destination
      if (this.hasTarget === false && this.state === true) {
        this.hasTarget = true;
        this.destination.x = targets[this.target].x;
        this.destination.y = targets[this.target].y;
        this.destination.z = targets[this.target].z;
      }
      
      //calculate the difference between current position and destination for each axis
      if (this.state === true) {
        var diffX = this.destination.x - this.position.x;
        var diffY = this.destination.y - this.position.y;
        var diffZ = this.destination.z - this.position.z;
     
        this.distance = hypot(diffX, diffY, diffZ);

        //this normalizes the vector, so our calculations for direction and speed in Cartesian system are not skewed
        if (this.distance > 0) {
          diffX = diffX / this.distance;
          diffY = diffY / this.distance;
          diffZ = diffZ / this.distance;
        }

        if (this.distance >= this.size) {
          //finds direction to the target in radians and convert it to degrees, y BEFORE x!!!
          var targetAngleXY = atan2(diffY, diffX) * (180 / PI);
          var targetAngleZ = atan(diffX/diffY) * (180 / PI);          
          
          //controls bot's altitude to be always bigger than its size, this is a crash prevention
          if (this.altitude <= this.size) {
            this.pitch += this.rotationSpeedZ + this.size;
          }
          
          //controls direction and turning speed of our bot, turning speed 0.1667 is one minute or 1/60 degree
          if (this.heading > targetAngleXY) {
            this.heading -= (this.rotationSpeedX + this.rotationSpeedY) + 1;
          } else if (this.heading < targetAngleXY) {
            this.heading += (this.rotationSpeedX + this.rotationSpeedY) + 1;
          }
          //this.heading = targetAngleXY;
          
          if (this.pitch > targetAngleZ) {
            this.pitch -= this.rotationSpeedZ + 1;
          } else if (this.pitch < targetAngleZ) {
            this.pitch += this.rotationSpeedZ + 1;
          }
          
          //calculate velocities for each axis
          this.velocity.x += diffX * this.acceleration * sim.timeStep;
          this.velocity.y += diffY * this.acceleration * sim.timeStep;
          this.velocity.z += diffZ * this.acceleration * sim.timeStep;
          
          //FIRST calculate Geopotential height
          //this is gravity-adjusted altitude of our bot, using variation of the gravity with latitude and elevation
          //based on https://en.wikipedia.org/wiki/Barometric_formula#Source_code
          sim.geopotentialHeight = (sim.earthRadius * this.altitude / (sim.earthRadius + this.altitude)) / 1000; //in kilometers

          //calculate current air temperature, in Kelvin
          //based on https://en.wikipedia.org/wiki/Barometric_formula#Source_code
          if (sim.geopotentialHeight <= 11) { sim.airCurrentTemp = 288.15 - (6.5 * sim.geopotentialHeight); } // Troposphere
          else if (sim.geopotentialHeight <= 20) { sim.airCurrentTemp = 216.65 - (6.5 * sim.geopotentialHeight); } // Stratosphere starts
          else if (sim.geopotentialHeight <= 32) { sim.airCurrentTemp = 196.65 + sim.geopotentialHeight; }
          else if (sim.geopotentialHeight <= 47) { sim.airCurrentTemp = 228.65 + 2.8 * (sim.geopotentialHeight - 32); }
          else if (sim.geopotentialHeight <= 51) { sim.airCurrentTemp = 270.65 - (6.5 * sim.geopotentialHeight); }// Mesosphere starts
          else if (sim.geopotentialHeight <= 71) { sim.airCurrentTemp = 270.65 - 2.8 * (sim.geopotentialHeight - 51); }
          else if (sim.geopotentialHeight <= 84.85) { sim.airCurrentTemp = 214.65 - 2 * (sim.geopotentialHeight - 71); }
          //geopotHeight must be less than 84.85 km
          var KelvinToRankine = sim.airCurrentTemp * 1.8; //in Rankine
          
          //calculate current air density and pressure         
          //in Pascals
          //based on https://en.wikipedia.org/wiki/Barometric_formula#Source_code
          if (sim.geopotentialHeight <= 11) { sim.airPressure = 101325 * pow(288.15 / sim.airCurrentTemp, -5.255877); }
          else if (sim.geopotentialHeight <= 20) { sim.airPressure = 22632.06 * exp(-0.1577 * (sim.geopotentialHeight - 11)); }
          else if (sim.geopotentialHeight <= 32) { sim.airPressure = 5474.889 * pow(216.65 / sim.airCurrentTemp, 34.16319); }
          else if (sim.geopotentialHeight <= 47) { sim.airPressure = 868.0187 * pow(228.65 / sim.airCurrentTemp, 12.2011); }
          else if (sim.geopotentialHeight <= 51) { sim.airPressure = 110.9063 * exp(-0.1262 * (sim.geopotentialHeight - 47)); }
          else if (sim.geopotentialHeight <= 71) { sim.airPressure = 66.93887 * pow(270.65 / sim.airCurrentTemp, -12.2011); }
          else if (sim.geopotentialHeight <= 84.85) { sim.airPressure = 3.956420 * pow(214.65 / sim.airCurrentTemp, -17.0816); }
          //altitude must be less than 86 km
          sim.airDensity = (sim.airPressure * sim.airMolarMass) / (sim.gasConstant * sim.airCurrentTemp);
          
          //calculate current air pressure at bot's position
          sim.airPressureX = 0.5 * sim.airDensity * pow(this.velocity.x, 2);
          sim.airPressureY = 0.5 * sim.airDensity * pow(this.velocity.y, 2);
          sim.airPressureZ = 0.5 * sim.airDensity * pow(this.velocity.z, 2);
          
          //calculate current air viscosity
          sim.airViscosity = sim.airReferenceViscosity*(((0.555*sim.RankineRefTemp) + sim.SutherlandConstant)/((0.555*KelvinToRankine) + sim.SutherlandConstant))*pow(KelvinToRankine/sim.RankineRefTemp, 3.2); //in centipose

          //calculate the Reynolds number so to use the correct drag law
          sim.kinematicViscosity = sim.airViscosity / sim.airDensity;
          sim.ReynoldsNumber = (this.acceleration * this.size) / sim.kinematicViscosity;
          
          //choose which drag law to use
          if (sim.ReynoldsNumber < 1) {
            //for low velocity, linear drag or laminar flow
            var dragX = 6 * PI * sim.airViscosity * this.size * this.velocity.x;
            var dragY = 6 * PI * sim.airViscosity * this.size * this.velocity.y;
            var dragZ = 6 * PI * sim.airViscosity * this.size * this.velocity.z;
          } else {
            //for high velocity, quadratic drag or turbulent flow
            var dragX = sim.airPressureX * this.dragCoeff * this.frontalArea;
            var dragY = sim.airPressureY * this.dragCoeff * this.frontalArea;
            var dragZ = sim.airPressureZ * this.dragCoeff * this.frontalArea;
          }
          dragX = (isNaN(dragX) ? 0 : dragX*sim.timeStep);
          dragY = (isNaN(dragY) ? 0 : dragY*sim.timeStep);
          dragZ = (isNaN(dragZ) ? 0 : dragZ*sim.timeStep);
          
          //calculate actual Earth's gravity
          var distanceFromEarthCenter = sim.earthRadius + this.position.z;
          sim.currentGravity = sim.referenceGravity * (sim.earthPowRadius / pow(distanceFromEarthCenter, 2));
          sim.currentGravity = parseFloat(sim.currentGravity).toFixed(6);
          
          //move our bot, maybe a Verlet version?
          this.acceleration *= 0.5*sim.timeStep;
          sim.currentGravity *= 0.5*sim.timeStep*sim.timeStep;
          this.velocity.x += diffX * ((this.acceleration - this.oldAcceleration)/2) * sim.timeStep;
          this.velocity.y += diffY * ((this.acceleration - this.oldAcceleration)/2) * sim.timeStep;
          this.velocity.z += diffZ * ((this.acceleration - this.oldAcceleration)/2) * sim.timeStep;
          var stepX = (this.velocity.x - dragX + this.acceleration);
          var stepY = (this.velocity.y - dragY + this.acceleration);
          var stepZ = (this.velocity.z - dragZ - sim.currentGravity + this.liftForce);// + this.acceleration);
          this.position.x += stepX;
          this.position.y += stepY;
          this.position.z += stepZ;
          this.oldAcceleration = this.acceleration;
          
          //calculate current speed and travelled distance of our bot
          var tmp = hypot(stepX, stepY, stepZ);
          this.distanceTravelled += tmp;
          this.speed = (tmp / sim.timeStep); //V = S / t, in m/s
          this.altitude = this.position.z;
};

Hope someone will show me where I did it wrong

_WeirdCat_

1,512

March 19, 2017 02:06 PM

you did almost everything wrong there are maybe drag equations for 1 or two things that MAY be correct, rest is not physically based sim,

at least you need to define, points of pressure for each element that causes lift or drag, you need skin drag and form drag + lift for each element, dont bother retards over internet telling you you need drag lift coefficients - you don't, you just define sin or cos for an angle and max/min angle when shape doesn't produce anything yeah sure thing that could be done by coefficient tables, but whatif you';re making your own aeroplane.

After that defined welcome to rigid body physics, the easiest way would be to define Ixx Iyy and Izz which you divide by torque generated by forces, anyway to achieve that you should define a constant mass for each object (wing, hull, flaps, airlerons, elevators, rudders etc.)

not quite sure center of mass will change over time when you fly (i guess yes)

theres alot of sh(Y)t roaming about reynolds number which actually will be needed


//---------------------------------------------------------------------------

#ifndef aeroH
#define aeroH
//---------------------------------------------------------------------------
//Mix of bullshitness
//don't even think it is correct in 0.0001%
#include <math.h>
#include "DxcMath.h"

const float AIR_STD_KINEMATIC_VISCOSITY = 1.460 * pow(10.0, -5.0); // m^2/s
const float WATER_STD_KINEMATIC_VISCOSITY = 1.3083 * pow(10.0, -6.0); // m^2/s             at temp = 10*C

inline float ReynoldsNum(float V, float choord_len, float kinematic_viscosity_of_fluid)
{
return (V*choord_len) / kinematic_viscosity_of_fluid;
}

				 //0.894

inline vec3 FindTriangleAerodynamicCenter(vec3 A, vec3 B, vec3 C)
{
}

inline float getChoordLength(vec3 velocity, vec3 surf_normal, vec3 pC, vec3 * verts, int vert_count, bool project, vec3 point_on_surf)
{
//Project velocity vector onto surface
vec3 projected = Normalize(ProjectVectorOnPlane(surf_normal, velocity));

//find the biggest side and square it - this will ensure us that ray will always hit two sides
float adst = -999.0;
for (int i=0; i < vert_count; i++)
{
int next = i + 1;
if (next == vert_count) next = 0;
adst = Max(adst, n3ddistance(verts[i],verts[next]));
}
adst = adst*adst;

//compute ray
vec3 rA = pC - projected*adst;
vec3 rB = pC + projected*adst;

vec3 cp;

for (int i=0; i < vert_count; i++)
cp = cp + verts[i];

cp = cp / float(vert_count);

vec3 n  = surf_normal*1000.0;

vec3 colp;
//we compute center point to determine whenever a side face is at front or back to the center point (because i don't need to waste time and run math on paper i can just check that and use it for further processing)
int cpside;
for (int i=0; i < vert_count; i++)
{
int next = i + 1;
if (next == vert_count) next = 0;

vec3 normal = Normal(verts[i], verts[next], verts[next] + n, true);
float dst = getplaneD(normal, verts[i]);
//cpside will now determine the
cpside = classifyapointagainstaplane(cp, normal, dst); // which side faces 'the inside' of polygon  (always either front or back never on plane) - solid convex polygon
if (cpside != isBack) //when center point is not behind side plane
{
normal = -normal;
dst = getplaneD(normal, verts[i]);
}

colp = rA;
if (SegmentPlaneIntersection(normal, dst, rA, rB, colp))
{
	//if ray and face normal 'look' (face) at each other
if (dot(Normalize(vectorAB(rA,rB)), normal) < 0) rA = colp;
else rB = colp;
}

}

return n3ddistance(rA, rB);
}



inline float getTriangleChoordLength(vec3 vel, vec3 surf_n, vec3 pC, vec3 A, vec3 B, vec3 C, bool project, vec3 pop)
{
vec3 ahue[3];
ahue[0] = A;
ahue[1] = B;
ahue[2] = C;
return getChoordLength(vel, surf_n, pC, ahue, 3, project, pop);
}

inline vec3 ShearStress(vec3 Force, float Area)
{
return Force / Area;
}

//1/ (pressure*v^2*Area) * ShearStress * dot(surf_normal, tangent_surf_dir)
inline float SkinFrictionDrag( float dens, float V, float Area, float ShearStress, t3dpoint<double> surf_normal, t3dpoint<double> surf_tangent_relative_to_airflow)
{
return (dens*V*V*Area) * ShearStress * absnf(dot(surf_normal, surf_tangent_relative_to_airflow));
}

inline float BodyDrag(float V, float Area, float fluidDens, t3dpoint<double> SurfNormal, t3dpoint<double> FluidFlowDir, bool front_only)
{

			  //now im not sure if its only the percentage of the rotation or if body is rotated the area increases
			  //anyways im not increasing the area.
float ratio = -dot(FluidFlowDir, SurfNormal);
if(!front_only) ratio = absnf(ratio);


return (0.5*(V*V)*fluidDens*Area*ratio);
}

#endif

i cant even say if all functions there are correct but reynolds num and getchoordlen should be fine.

now take a look at this for each element that causes drag or lift (i use only drag here since its a part of mu sail simulation)



for (int i=0; i < submerged_tri_count; i++)
	//now iterate through form_percentage_coefff
if (ACTUAL_FRAME[i].surface > 0.1) //do not coun't anything that is less than 10 cm^2
{
	whole_sub_surf = whole_sub_surf + ACTUAL_FRAME[i].surface;
	vec3 element_form_drag = vec3(0.0, 0.0, 0.0); //somehow i dont trust that compiler due to optimization will zero it
	vec3 element_skin_drag = vec3(0.0, 0.0, 0.0);
	vec3 cog2pC 	= vectorAB(pos + ROTATION_MAT * CENTER_OF_GRAVITY, ACTUAL_FRAME[i].pC);
	vec3 local_vel = vel + AngVel * cog2pC;
	vec3 vn = Normalize(local_vel);
	float vv = VectorLength(local_vel);
	float v2 = vv;
	vv = vv * vv;
float form_percentage_coeff = absnf( acos(dot(ACTUAL_FRAME[i].normal, vn)) / float(pi) );
if (dot(ACTUAL_FRAME[i].normal, vn) <= 0.0) form_percentage_coeff = 0; //suction force

element_form_drag = (-vn) * (form_percentage_coeff * (999.1026 * 0.5 * vv) * ACTUAL_FRAME[i].surface); //form drag




float skin_percentage_coeff = absnf(1.0 - form_percentage_coeff); //due to float inaccuracy
vec3 surf_n = Normal(ACTUAL_FRAME[i].V[0], ACTUAL_FRAME[i].V[1], ACTUAL_FRAME[i].V[2], true);
float surf_choordlen = getTriangleChoordLength(
-local_vel, surf_n, ACTUAL_FRAME[i].pC,
ACTUAL_FRAME[i].V[0], ACTUAL_FRAME[i].V[1], ACTUAL_FRAME[i].V[2], true, ACTUAL_FRAME[i].pC);

if (surf_choordlen > 20.0) surf_choordlen = 1.0;

//calc reynolds number
float Rn = ReynoldsNum(VectorLength(local_vel), surf_choordlen, 0.894);
bool ah = false;
if (absnf(Rn) <= 0.001 ) {Rn = 1.0; ah= true;}

float CFRn = 1.328 / (sqrt(Rn));
if (ah) CFRn = 0.0;
//float CFRn = 0.075 / (lg * lg); //1.328 / (sqrt(Rn));//

vec3 Vfi = Normalize(ProjectVectorOnPlane(ACTUAL_FRAME[i].normal, -local_vel));
element_skin_drag = (-vn) * (CFRn * ((999.1026 * vv) / 2.0 ) * ACTUAL_FRAME[i].surface * skin_percentage_coeff);



skin_drag = skin_drag + element_skin_drag;
form_drag = form_drag + element_form_drag;

cog2cob_vec = vectorAB(pos + ROTATION_MAT * CENTER_OF_GRAVITY, ACTUAL_FRAME[i].pC);

ETorque = cog2cob_vec	*  ( element_skin_drag + element_form_drag );

elements_torque = elements_torque + ETorque;

}


drag_force = form_drag + skin_drag;

ACTUAL_FRAME.pC pC means pressure center.

although youll need to check if equations apply to speeds above 1 mach.

after finding all forces that act on your aeroplane you should apply them to center of mass and torque

then you can do

vec3 accel = result_force / mass;

vel = vel + accel*dt;

pos = pos + vel*dt;

after all that theres another almost the most important thing DAMPING - without it simulations tend to be very fkn unstable like a unicron without that di(k on its head


void damp( vec3 * ptr, float amount, float epsilon)
{
	(*ptr) = (*ptr) - (*ptr)/amount;
	vec3 p = vec3(absnf<float>((*ptr).x),absnf<float>((*ptr).y),absnf<float>((*ptr).z));
	if (IsNan(p.x)) p.x = 0.0;
	if (IsNan(p.y)) p.y = 0.0;
	if (IsNan(p.z)) p.z = 0.0;

	if ( (p.x > 0) && (p.x < epsilon) ) (*ptr).x = 0.0;
	if ( (p.y > 0) && (p.y < epsilon) ) (*ptr).y = 0.0;
	if ( (p.z > 0) && (p.z < epsilon) ) (*ptr).z = 0.0;
}



inline void damp( float * ptr, float amount, float epsilon)
{
	(*ptr) = (*ptr) - (*ptr)/amount;

	if (IsNan((*ptr))) (*ptr) = 0.0;
	if ( ((*ptr) > 0) && ((*ptr) < epsilon) ) (*ptr) = 0.0;
}

where you damp things before physics frame

damp(&vel, 10.0, 0.007);

damp(&AngVel, 10.0, 0.0005);

this is important because you'll have to check if you get nans or numbers close to zero otherwise you get nans or extremely big numbers and it will screw your simulation. that 10 in function parameter is just a constant that should be fine, you can change it.

Anyway to obtain Ixx, Iyy and Izz you will need a fancy function which is actually avail for free in cpp but i forgot the author and source, anyway it says that your model trangles have to be counterclockwise, its easy to set it up

you find a center of each triangle defining an element then you move outwards ( face normal), now using matrix math (projection and view matrix) you need to put that triangle to fit your view matrix, then you sort your vertices to be counterclockwise

sraps of code:



void MakeModelTriCounterClockWise(Matrix44<T> proj_mat, T fovy, T wh_aspect, float SCRW, float SCRH)
{
Matrix44<T> world_mat;
T vmp[16];
for (int i=0; i < FaceLength; i++)
{
Matrix44<T> view_mat;
view_mat.LoadIdentity();

int index = VBO_BE[i].INDEX_START;

t3dpoint<T> up_vec = Normalize(vectorAB(FACE_CENTER_POINT[i], AOS[index].v));


T side_dst = -1;

side_dst = Max(VectorLength(vectorAB(AOS[index].v, AOS[index+1].v)), VectorLength(vectorAB(AOS[index+1].v, AOS[index+2].v)));
side_dst = Max(VectorLength(vectorAB(AOS[index+2].v, AOS[index].v)), side_dst);
side_dst = side_dst * side_dst;

T fov = fovy * wh_aspect;
T alpha = fov / 2.0;

//tan(alpha) = (side_dst / 2.0) / x; //we need to find x
T tg = tan(alpha*imopi);

if (absnf(tg) <= 0.0001) tg = 0.001;
T cam_distance = (side_dst / 2.0) / tg;

t3dpoint<T> eyepos = FACE_CENTER_POINT[i] + FACE_N[i]*cam_distance;


glLookAt(view_mat, eyepos, FACE_CENTER_POINT[i], up_vec);

Matrix44<T> mvp = (world_mat * view_mat) * proj_mat;


vec4 A = mvp % vec4(AOS[index].v, 1.0);
A.x = A.x / A.w;
A.y = A.y / A.w;
A = A * 0.5 + 0.5;
A.x = A.x * SCRW;
A.y = A.y * SCRH;


vec4 B = mvp % vec4(AOS[index+1].v, 1.0);
B.x = B.x / B.w;
B.y = B.y / B.w;
B = B * 0.5 + 0.5;
B.x = B.x * SCRW;
B.y = B.y * SCRH;


vec4 C = mvp % vec4(AOS[index+2].v, 1.0);
C.x = C.x / C.w;
C.y = C.y / C.w;
C = C * 0.5 + 0.5;
C.x = C.x * SCRW;
C.y = C.y * SCRH;



if (IsTriClockwise(vec3(A.x, A.y, A.z), vec3(B.x, B.y, B.z), vec3(C.x, C.y, C.z))) //swap it to be counter clockwise
{
	TTGLVertex<T,float,float> v1 = AOS[index+2];
	TTGLVertex<T,float,float> v3 = AOS[index];
AOS[index]     = v1;
AOS[index+2]   = v3;
}

}
}










template <class T> bool IsTriClockwise(t3dpoint<T> A, t3dpoint<T> B, t3dpoint<T> C)
{
    T sum = T(0.0);
	  sum = sum + (B.x - A.x) * (B.y + A.y);
	  sum = sum + (C.x - B.x) * (C.y + B.y);
	  sum = sum + (A.x - C.x) * (A.y + C.y);
    return (sum > 0.0);
}

if you get anything out of it then you can say something about flight sim. but you use java i hope your computer burns with all your code. best of no luck

and you are only at half or even less of things you need to define.

Master & Mentor
https://sites.google.com/site/customprog/

common_swift

124

Author

March 19, 2017 05:10 PM

Whoa, thanks ...

Yeah, you're right about skin drag and form drag + points of pressure but this will be added later (at least hard-coded classic aircraft designs), I just started from somewhere to have some base and improve it over time - still not reached the moment for adding real elevators, ailerons, etc, at first need a properly working simulation

common_swift

124

Author

March 20, 2017 01:37 AM

Purely theoretical, what will looks a step-by-step calculation list with all physical quantities required for an accurate flight simulation? For the question lets assume we are starting from finding object's current mass, and we know we are at 100% throttle (always max thrust) ... ie at step 1 I'm finding current mass, then at step 2 the wing loading/lift force, at step 3 current AoA/roll/yaw/heading/altitude of the plane and angle between origin and destination point, etc ...

Working on small flight sim game, need help with equations

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Working on small flight sim game, need help with equations

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Reticulating splines