This was masked because I didn't wipe src/FDM/libFlight.a and recreate it.
+++ /dev/null
-// IO360.cxx - a piston engine model currently for the IO360 engine fitted to the C172
-// but with the potential to model other naturally aspirated piston engines
-// given appropriate config input.
-//
-// Written by David Luff, started 2000.
-// Based on code by Phil Schubert, started 1999.
-//
-// This program is free software; you can redistribute it and/or
-// modify it under the terms of the GNU General Public License as
-// published by the Free Software Foundation; either version 2 of the
-// License, or (at your option) any later version.
-//
-// This program is distributed in the hope that it will be useful, but
-// WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-// General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
-
-#include <simgear/compiler.h>
-
-#include <math.h>
-
-#include STL_FSTREAM
-#include STL_IOSTREAM
-
-#include "IO360.hxx"
-#include "LaRCsim/ls_constants.h"
-
-#include <Main/fg_props.hxx>
-
-//*************************************************************************************
-// Initialise the engine model
-void FGNewEngine::init(double dt) {
-
- // These constants should probably be moved eventually
- CONVERT_CUBIC_INCHES_TO_METERS_CUBED = 1.638706e-5;
- CONVERT_HP_TO_WATTS = 745.6999;
-
- // Properties of working fluids
- Cp_air = 1005; // J/KgK
- Cp_fuel = 1700; // J/KgK
- calorific_value_fuel = 47.3e6; // W/Kg Note that this is only an approximate value
- rho_fuel = 800; // kg/m^3 - an estimate for now
- R_air = 287.3;
-
- // environment inputs
- p_amb_sea_level = 101325; // Pascals
-
- // Control inputs - ARE THESE NEEDED HERE???
- Throttle_Lever_Pos = 75;
- Propeller_Lever_Pos = 75;
- Mixture_Lever_Pos = 100;
-
- //misc
- IAS = 0;
- time_step = dt;
-
- // Engine Specific Variables that should be read in from a config file
- MaxHP = 200; //Lycoming IO360 -A-C-D series
-// MaxHP = 180; //Current Lycoming IO360 ?
-// displacement = 520; //Continental IO520-M
- displacement = 360; //Lycoming IO360
- displacement_SI = displacement * CONVERT_CUBIC_INCHES_TO_METERS_CUBED;
- engine_inertia = 0.2; //kgm^2 - value taken from a popular family saloon car engine - need to find an aeroengine value !!!!!
- prop_inertia = 0.05; //kgm^2 - this value is a total guess - dcl
- Max_Fuel_Flow = 130; // Units??? Do we need this variable any more??
-
- // Engine specific variables that maybe should be read in from config but are pretty generic and won't vary much for a naturally aspirated piston engine.
- Max_Manifold_Pressure = 28.50; //Inches Hg. An approximation - should be able to find it in the engine performance data
- Min_Manifold_Pressure = 6.5; //Inches Hg. This is a guess corresponding to approx 0.24 bar MAP (7 in Hg) - need to find some proper data for this
- Max_RPM = 2700;
- Min_RPM = 600; //Recommended idle from Continental data sheet
- Mag_Derate_Percent = 5;
- Gear_Ratio = 1;
- n_R = 2; // Number of crank revolutions per power cycle - 2 for a 4 stroke engine.
-
- // Various bits of housekeeping describing the engines initial state.
- running = false;
- cranking = false;
- crank_counter = false;
- starter = false;
-
- // Initialise Engine Variables used by this instance
- if(running)
- RPM = 600;
- else
- RPM = 0;
- Percentage_Power = 0;
- Manifold_Pressure = 29.96; // Inches
- Fuel_Flow_gals_hr = 0;
-// Torque = 0;
- Torque_SI = 0;
- CHT = 298.0; //deg Kelvin
- CHT_degF = (CHT_degF * 1.8) - 459.67; //deg Fahrenheit
- Mixture = 14;
- Oil_Pressure = 0; // PSI
- Oil_Temp = 85; // Deg C
- current_oil_temp = 298.0; //deg Kelvin
- /**** one of these is superfluous !!!!***/
- HP = 0;
- RPS = 0;
- Torque_Imbalance = 0;
-
- // Initialise Propellor Variables used by this instance
- FGProp1_RPS = 0;
- // Hardcode propellor for now - the following two should be read from config eventually
- prop_diameter = 1.8; // meters
- blade_angle = 23.0; // degrees
-}
-
-//*****************************************************************************
-// update the engine model based on current control positions
-void FGNewEngine::update() {
-
-/*
- // Hack for testing - should output every 5 seconds
- static int count1 = 0;
- if(count1 == 0) {
-// cout << "P_atmos = " << p_amb << " T_atmos = " << T_amb << '\n';
-// cout << "Manifold pressure = " << Manifold_Pressure << " True_Manifold_Pressure = " << True_Manifold_Pressure << '\n';
-// cout << "p_amb_sea_level = " << p_amb_sea_level << '\n';
-// cout << "equivalence_ratio = " << equivalence_ratio << '\n';
-// cout << "combustion_efficiency = " << combustion_efficiency << '\n';
-// cout << "AFR = " << 14.7 / equivalence_ratio << '\n';
-// cout << "Mixture lever = " << Mixture_Lever_Pos << '\n';
-// cout << "n = " << RPM << " rpm\n";
-// cout << "T_amb = " << T_amb << '\n';
-// cout << "running = " << running << '\n';
-// cout << "fuel = " << fgGetFloat("/consumables/fuel/tank[0]/level-gal_us") << '\n';
-// cout << "Percentage_Power = " << Percentage_Power << '\n';
-// cout << "current_oil_temp = " << current_oil_temp << '\n';
-// cout << "EGT = " << EGT << '\n';
- }
- count1++;
- if(count1 == 100)
- count1 = 0;
-*/
-
- // Check parameters that may alter the operating state of the engine.
- // (spark, fuel, starter motor etc)
-
- // Check for spark
- bool Magneto_Left = false;
- bool Magneto_Right = false;
- // Magneto positions:
- // 0 -> off
- // 1 -> left only
- // 2 -> right only
- // 3 -> both
- if(mag_pos != 0) {
- spark = true;
- } else {
- spark = false;
- } // neglects battery voltage, master on switch, etc for now.
- if((mag_pos == 1) || (mag_pos > 2))
- Magneto_Left = true;
- if(mag_pos > 1)
- Magneto_Right = true;
-
- // crude check for fuel
- if((fgGetFloat("/consumables/fuel/tank[0]/level-gal_us") > 0) || (fgGetFloat("/consumables/fuel/tank[1]/level-gal_us") > 0)) {
- fuel = true;
- } else {
- fuel = false;
- } // Need to make this better, eg position of fuel selector switch.
-
- // Check if we are turning the starter motor
- if(cranking != starter) {
- // This check saves .../cranking from getting updated every loop - they only update when changed.
- cranking = starter;
- crank_counter = 0;
- }
- // Note that although /engines/engine[0]/starter and /engines/engine[0]/cranking might appear to be duplication it is
- // not since the starter may be engaged with the battery voltage too low for cranking to occur (or perhaps the master
- // switch just left off) and the sound manager will read .../cranking to determine wether to play a cranking sound.
- // For now though none of that is implemented so cranking can be set equal to .../starter without further checks.
-
-// int Alternate_Air_Pos =0; // Off = 0. Reduces power by 3 % for same throttle setting
- // DCL - don't know what this Alternate_Air_Pos is - this is a leftover from the Schubert code.
-
- //Check mode of engine operation
- if(cranking) {
- crank_counter++;
- if(RPM <= 480) {
- RPM += 100;
- if(RPM > 480)
- RPM = 480;
- } else {
- // consider making a horrible noise if the starter is engaged with the engine running
- }
- }
- if((!running) && (spark) && (fuel) && (crank_counter > 120)) {
- // start the engine if revs high enough
- if(RPM > 450) {
- // For now just instantaneously start but later we should maybe crank for a bit
- running = true;
-// RPM = 600;
- }
- }
- if( (running) && ((!spark)||(!fuel)) ) {
- // Cut the engine
- // note that we only cut the power - the engine may continue to spin if the prop is in a moving airstream
- running = false;
- }
-
- // Now we've ascertained whether the engine is running or not we can start to do the engine calculations 'proper'
-
- // Calculate Sea Level Manifold Pressure
- Manifold_Pressure = Calc_Manifold_Pressure( Throttle_Lever_Pos, Max_Manifold_Pressure, Min_Manifold_Pressure );
- // cout << "manifold pressure = " << Manifold_Pressure << endl;
-
- //Then find the actual manifold pressure (the calculated one is the sea level pressure)
- True_Manifold_Pressure = Manifold_Pressure * p_amb / p_amb_sea_level;
-
- //Do the fuel flow calculations
- Calc_Fuel_Flow_Gals_Hr();
-
- //Calculate engine power
- Calc_Percentage_Power(Magneto_Left, Magneto_Right);
- HP = Percentage_Power * MaxHP / 100.0;
- Power_SI = HP * CONVERT_HP_TO_WATTS;
-
- // FMEP calculation. For now we will just use this during motored operation.
- // Eventually we will calculate IMEP and use the FMEP all the time to give BMEP (maybe!)
- if(!running) {
- // This FMEP data is from the Patton paper, assumes fully warm conditions.
- FMEP = 1e-12*pow(RPM,4) - 1e-8*pow(RPM,3) + 5e-5*pow(RPM,2) - 0.0722*RPM + 154.85;
- // Gives FMEP in kPa - now convert to Pa
- FMEP *= 1000;
- } else {
- FMEP = 0.0;
- }
- // Is this total FMEP or friction FMEP ???
-
- Torque_FMEP = (FMEP * displacement_SI) / (2.0 * LS_PI * n_R);
-
- // Calculate Engine Torque. Check for div by zero since percentage power correlation does not guarantee zero power at zero rpm.
- // However this is problematical since there is a resistance to movement even at rest
- // Ie this is a dynamics equation not a statics one. This can be solved by going over to MEP based torque calculations.
- if(RPM == 0) {
- Torque_SI = 0 - Torque_FMEP;
- }
- else {
- Torque_SI = ((Power_SI * 60.0) / (2.0 * LS_PI * RPM)) - Torque_FMEP; //Torque = power / angular velocity
- // cout << Torque << " Nm\n";
- }
-
- //Calculate Exhaust gas temperature
- if(running)
- Calc_EGT();
- else
- EGT = 298.0;
-
- // Calculate Cylinder Head Temperature
- Calc_CHT();
-
- // Calculate oil temperature
- current_oil_temp = Calc_Oil_Temp(current_oil_temp);
-
- // Calculate Oil Pressure
- Oil_Pressure = Calc_Oil_Press( Oil_Temp, RPM );
-
- // Now do the Propeller Calculations
- Do_Prop_Calcs();
-
-// Now do the engine - prop torque balance to calculate final RPM
-
- //Calculate new RPM from torque balance and inertia.
- Torque_Imbalance = Torque_SI - prop_torque; //This gives a +ve value when the engine torque exeeds the prop torque
- // (Engine torque is +ve when it acts in the direction of engine revolution, prop torque is +ve when it opposes the direction of engine revolution)
-
- angular_acceleration = Torque_Imbalance / (engine_inertia + prop_inertia);
- angular_velocity_SI += (angular_acceleration * time_step);
- // Don't let the engine go into reverse
- if(angular_velocity_SI < 0)
- angular_velocity_SI = 0;
- RPM = (angular_velocity_SI * 60) / (2.0 * LS_PI);
-
- // And finally a last check on the engine state after doing the torque balance with the prop - have we stalled?
- if(running) {
- //Check if we have stalled the engine
- if (RPM == 0) {
- running = false;
- } else if((RPM <= 480) && (cranking)) {
- //Make sure the engine noise dosn't play if the engine won't start due to eg mixture lever pulled out.
- running = false;
- EGT = 298.0;
- }
- }
-
- // And finally, do any unit conversions from internal units to output units
- EGT_degF = (EGT * 1.8) - 459.67;
- CHT_degF = (CHT * 1.8) - 459.67;
-}
-
-//*****************************************************************************************************
-
-// FGNewEngine member functions
-
-////////////////////////////////////////////////////////////////////////////////////////////
-// Return the combustion efficiency as a function of equivalence ratio
-//
-// Combustion efficiency values from Heywood,
-// "Internal Combustion Engine Fundamentals", ISBN 0-07-100499-8
-////////////////////////////////////////////////////////////////////////////////////////////
-float FGNewEngine::Lookup_Combustion_Efficiency(float thi_actual)
-{
- const int NUM_ELEMENTS = 11;
- float thi[NUM_ELEMENTS] = {0.0, 0.9, 1.0, 1.05, 1.1, 1.15, 1.2, 1.3, 1.4, 1.5, 1.6}; //array of equivalence ratio values
- float neta_comb[NUM_ELEMENTS] = {0.98, 0.98, 0.97, 0.95, 0.9, 0.85, 0.79, 0.7, 0.63, 0.57, 0.525}; //corresponding array of combustion efficiency values
- float neta_comb_actual = 0.0f;
- float factor;
-
- int i;
- int j = NUM_ELEMENTS; //This must be equal to the number of elements in the lookup table arrays
-
- for(i=0;i<j;i++)
- {
- if(i == (j-1)) {
- // Assume linear extrapolation of the slope between the last two points beyond the last point
- float dydx = (neta_comb[i] - neta_comb[i-1]) / (thi[i] - thi[i-1]);
- neta_comb_actual = neta_comb[i] + dydx * (thi_actual - thi[i]);
- return neta_comb_actual;
- }
- if(thi_actual == thi[i]) {
- neta_comb_actual = neta_comb[i];
- return neta_comb_actual;
- }
- if((thi_actual > thi[i]) && (thi_actual < thi[i + 1])) {
- //do linear interpolation between the two points
- factor = (thi_actual - thi[i]) / (thi[i+1] - thi[i]);
- neta_comb_actual = (factor * (neta_comb[i+1] - neta_comb[i])) + neta_comb[i];
- return neta_comb_actual;
- }
- }
-
- //if we get here something has gone badly wrong
- SG_LOG(SG_AIRCRAFT, SG_ALERT, "error in FGNewEngine::Lookup_Combustion_Efficiency");
- return neta_comb_actual;
-}
-
-////////////////////////////////////////////////////////////////////////////////////////////
-// Return the percentage of best mixture power available at a given mixture strength
-//
-// Based on data from "Technical Considerations for Catalysts for the European Market"
-// by H S Gandi, published 1988 by IMechE
-//
-// Note that currently no attempt is made to set a lean limit on stable combustion
-////////////////////////////////////////////////////////////////////////////////////////////
-float FGNewEngine::Power_Mixture_Correlation(float thi_actual)
-{
- float AFR_actual = 14.7 / thi_actual;
- // thi and thi_actual are equivalence ratio
- const int NUM_ELEMENTS = 13;
- // The lookup table is in AFR because the source data was. I added the two end elements to make sure we are almost always in it.
- float AFR[NUM_ELEMENTS] = {(14.7/1.6), 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, (14.7/0.6)}; //array of equivalence ratio values
- float mixPerPow[NUM_ELEMENTS] = {78, 86, 93.5, 98, 100, 99, 96.4, 92.5, 88, 83, 78.5, 74, 58}; //corresponding array of combustion efficiency values
- float mixPerPow_actual = 0.0f;
- float factor;
- float dydx;
-
- int i;
- int j = NUM_ELEMENTS; //This must be equal to the number of elements in the lookup table arrays
-
- for(i=0;i<j;i++)
- {
- if(i == (j-1)) {
- // Assume linear extrapolation of the slope between the last two points beyond the last point
- dydx = (mixPerPow[i] - mixPerPow[i-1]) / (AFR[i] - AFR[i-1]);
- mixPerPow_actual = mixPerPow[i] + dydx * (AFR_actual - AFR[i]);
- return mixPerPow_actual;
- }
- if((i == 0) && (AFR_actual < AFR[i])) {
- // Assume linear extrapolation of the slope between the first two points for points before the first point
- dydx = (mixPerPow[i] - mixPerPow[i-1]) / (AFR[i] - AFR[i-1]);
- mixPerPow_actual = mixPerPow[i] + dydx * (AFR_actual - AFR[i]);
- return mixPerPow_actual;
- }
- if(AFR_actual == AFR[i]) {
- mixPerPow_actual = mixPerPow[i];
- return mixPerPow_actual;
- }
- if((AFR_actual > AFR[i]) && (AFR_actual < AFR[i + 1])) {
- //do linear interpolation between the two points
- factor = (AFR_actual - AFR[i]) / (AFR[i+1] - AFR[i]);
- mixPerPow_actual = (factor * (mixPerPow[i+1] - mixPerPow[i])) + mixPerPow[i];
- return mixPerPow_actual;
- }
- }
-
- //if we get here something has gone badly wrong
- SG_LOG(SG_AIRCRAFT, SG_ALERT, "error in FGNewEngine::Power_Mixture_Correlation");
- return mixPerPow_actual;
-}
-
-// Calculate Cylinder Head Temperature
-// Crudely models the cylinder head as an arbitary lump of arbitary size and area with one third of combustion energy
-// as heat input and heat output as a function of airspeed and temperature. Could be improved!!!
-void FGNewEngine::Calc_CHT()
-{
- float h1 = -95.0; //co-efficient for free convection
- float h2 = -3.95; //co-efficient for forced convection
- float h3 = -0.05; //co-efficient for forced convection due to prop backwash
- float v_apparent; //air velocity over cylinder head in m/s
- float v_dot_cooling_air;
- float m_dot_cooling_air;
- float temperature_difference;
- float arbitary_area = 1.0;
- float dqdt_from_combustion;
- float dqdt_forced; //Rate of energy transfer to/from cylinder head due to forced convection (Joules) (sign convention: to cylinder head is +ve)
- float dqdt_free; //Rate of energy transfer to/from cylinder head due to free convection (Joules) (sign convention: to cylinder head is +ve)
- float dqdt_cylinder_head; //Overall energy change in cylinder head
- float CpCylinderHead = 800.0; //FIXME - this is a guess - I need to look up the correct value
- float MassCylinderHead = 8.0; //Kg - this is a guess - it dosn't have to be absolutely accurate but can be varied to alter the warm-up rate
- float HeatCapacityCylinderHead;
- float dCHTdt;
-
- // The above values are hardwired to give reasonable results for an IO360 (C172 engine)
- // Now adjust to try to compensate for arbitary sized engines - this is currently untested
- arbitary_area *= (MaxHP / 180.0);
- MassCylinderHead *= (MaxHP / 180.0);
-
- temperature_difference = CHT - T_amb;
-
- v_apparent = IAS * 0.5144444; //convert from knots to m/s
- v_dot_cooling_air = arbitary_area * v_apparent;
- m_dot_cooling_air = v_dot_cooling_air * rho_air;
-
- //Calculate rate of energy transfer to cylinder head from combustion
- dqdt_from_combustion = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
-
- //Calculate rate of energy transfer from cylinder head due to cooling NOTE is calculated as rate to but negative
- dqdt_forced = (h2 * m_dot_cooling_air * temperature_difference) + (h3 * RPM * temperature_difference);
- dqdt_free = h1 * temperature_difference;
-
- //Calculate net rate of energy transfer to or from cylinder head
- dqdt_cylinder_head = dqdt_from_combustion + dqdt_forced + dqdt_free;
-
- HeatCapacityCylinderHead = CpCylinderHead * MassCylinderHead;
-
- dCHTdt = dqdt_cylinder_head / HeatCapacityCylinderHead;
-
- CHT += (dCHTdt * time_step);
-}
-
-// Calculate exhaust gas temperature
-void FGNewEngine::Calc_EGT()
-{
- combustion_efficiency = Lookup_Combustion_Efficiency(equivalence_ratio); //The combustion efficiency basically tells us what proportion of the fuels calorific value is released
-
- //now calculate energy release to exhaust
- //We will assume a three way split of fuel energy between useful work, the coolant system and the exhaust system
- //This is a reasonable first suck of the thumb estimate for a water cooled automotive engine - whether it holds for an air cooled aero engine is probably open to question
- //Regardless - it won't affect the variation of EGT with mixture, and we can always put a multiplier on EGT to get a reasonable peak value.
- enthalpy_exhaust = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
- heat_capacity_exhaust = (Cp_air * m_dot_air) + (Cp_fuel * m_dot_fuel);
- delta_T_exhaust = enthalpy_exhaust / heat_capacity_exhaust;
-
- EGT = T_amb + delta_T_exhaust;
-
- //The above gives the exhaust temperature immediately prior to leaving the combustion chamber
- //Now derate to give a more realistic figure as measured downstream
- //For now we will aim for a peak of around 400 degC (750 degF)
-
- EGT *= 0.444 + ((0.544 - 0.444) * Percentage_Power / 100.0);
-}
-
-// Calculate Manifold Pressure based on Throttle lever Position
-float FGNewEngine::Calc_Manifold_Pressure ( float LeverPosn, float MaxMan, float MinMan)
-{
- float Inches;
-
- //Note that setting the manifold pressure as a function of lever position only is not strictly accurate
- //MAP is also a function of engine speed. (and ambient pressure if we are going for an actual MAP model)
- Inches = MinMan + (LeverPosn * (MaxMan - MinMan) / 100);
-
- //allow for idle bypass valve or slightly open throttle stop
- if(Inches < MinMan)
- Inches = MinMan;
-
- //Check for stopped engine (crudest way of compensating for engine speed)
- if(RPM == 0)
- Inches = 29.92;
-
- return Inches;
-}
-
-// Calculate fuel flow in gals/hr
-void FGNewEngine::Calc_Fuel_Flow_Gals_Hr()
-{
- //DCL - calculate mass air flow into engine based on speed and load - separate this out into a function eventually
- //t_amb is actual temperature calculated from altitude
- //calculate density from ideal gas equation
- rho_air = p_amb / ( R_air * T_amb );
- rho_air_manifold = rho_air * Manifold_Pressure / 29.6; //This is a bit of a roundabout way of calculating this but it works !! If we put manifold pressure into SI units we could do it simpler.
- //calculate ideal engine volume inducted per second
- swept_volume = (displacement_SI * (RPM / 60)) / 2; //This equation is only valid for a four stroke engine
- //calculate volumetric efficiency - for now we will just use 0.8, but actually it is a function of engine speed and the exhaust to manifold pressure ratio
- //Note that this is cylinder vol eff - the throttle drop is already accounted for in the MAP calculation
- volumetric_efficiency = 0.8;
- //Now use volumetric efficiency to calculate actual air volume inducted per second
- v_dot_air = swept_volume * volumetric_efficiency;
- //Now calculate mass flow rate of air into engine
- m_dot_air = v_dot_air * rho_air_manifold;
-
-//**************
-
- //DCL - now calculate fuel flow into engine based on air flow and mixture lever position
- //assume lever runs from no flow at fully out to thi = 1.3 at fully in at sea level
- //also assume that the injector linkage is ideal - hence the set mixture is maintained at a given altitude throughout the speed and load range
- thi_sea_level = 1.3 * ( Mixture_Lever_Pos / 100.0 );
- equivalence_ratio = thi_sea_level * p_amb_sea_level / p_amb; //ie as we go higher the mixture gets richer for a given lever position
- m_dot_fuel = m_dot_air / 14.7 * equivalence_ratio;
- Fuel_Flow_gals_hr = (m_dot_fuel / rho_fuel) * 264.172 * 3600.0; // Note this assumes US gallons
-}
-
-// Calculate the percentage of maximum rated power delivered as a function of Manifold pressure multiplied by engine speed (rpm)
-// This is not necessarilly the best approach but seems to work for now.
-// May well need tweaking at the bottom end if the prop model is changed.
-void FGNewEngine::Calc_Percentage_Power(bool mag_left, bool mag_right)
-{
- // For a given Manifold Pressure and RPM calculate the % Power
- // Multiply Manifold Pressure by RPM
- float ManXRPM = True_Manifold_Pressure * RPM;
-
-/*
-// Phil's %power correlation
- // Calculate % Power
- Percentage_Power = (+ 7E-09 * ManXRPM * ManXRPM) + ( + 7E-04 * ManXRPM) - 0.1218;
- // cout << Percentage_Power << "%" << "\t";
-*/
-
-// DCL %power correlation - basically Phil's correlation modified to give slighty less power at the low end
-// might need some adjustment as the prop model is adjusted
-// My aim is to match the prop model and engine model at the low end to give the manufacturer's recommended idle speed with the throttle closed - 600rpm for the Continental IO520
- // Calculate % Power for Nev's prop model
- //Percentage_Power = (+ 6E-09 * ManXRPM * ManXRPM) + ( + 8E-04 * ManXRPM) - 1.8524;
-
- // Calculate %power for DCL prop model
- Percentage_Power = (7e-9 * ManXRPM * ManXRPM) + (7e-4 * ManXRPM) - 1.0;
-
- // Power de-rating for altitude has been removed now that we are basing the power
- // on the actual manifold pressure, which takes air pressure into account. However - this fails to
- // take the temperature into account - this is TODO.
-
- // Adjust power for temperature - this is temporary until the power is done as a function of mass flow rate induced
- // Adjust for Temperature - Temperature above Standard decrease
- // power by 7/120 % per degree F increase, and incease power for
- // temps below at the same ratio
- float T_amb_degF = (T_amb * 1.8) - 459.67;
- float T_amb_sea_lev_degF = (288 * 1.8) - 459.67;
- Percentage_Power = Percentage_Power + ((T_amb_sea_lev_degF - T_amb_degF) * 7 /120);
-
- //DCL - now adjust power to compensate for mixture
- Percentage_of_best_power_mixture_power = Power_Mixture_Correlation(equivalence_ratio);
- Percentage_Power = Percentage_Power * Percentage_of_best_power_mixture_power / 100.0;
-
- // Now Derate engine for the effects of Bad/Switched off magnetos
- //if (Magneto_Left == 0 && Magneto_Right == 0) {
- if(!running) {
- // cout << "Both OFF\n";
- Percentage_Power = 0;
- } else if (mag_left && mag_right) {
- // cout << "Both On ";
- } else if (mag_left == 0 || mag_right== 0) {
- // cout << "1 Magneto Failed ";
- Percentage_Power = Percentage_Power * ((100.0 - Mag_Derate_Percent)/100.0);
- // cout << FGEng1_Percentage_Power << "%" << "\t";
- }
-/*
- //DCL - stall the engine if RPM drops below 450 - this is possible if mixture lever is pulled right out
- //This is a kludge that I should eliminate by adding a total fmep estimation.
- if(RPM < 450)
- Percentage_Power = 0;
-*/
- if(Percentage_Power < 0)
- Percentage_Power = 0;
-}
-
-// Calculate Oil Temperature in degrees Kelvin
-float FGNewEngine::Calc_Oil_Temp (float oil_temp)
-{
- float idle_percentage_power = 2.3; // approximately
- float target_oil_temp; // Steady state oil temp at the current engine conditions
- float time_constant; // The time constant for the differential equation
- if(running) {
- target_oil_temp = 363;
- time_constant = 500; // Time constant for engine-on idling.
- if(Percentage_Power > idle_percentage_power) {
- time_constant /= ((Percentage_Power / idle_percentage_power) / 10.0); // adjust for power
- }
- } else {
- target_oil_temp = 298;
- time_constant = 1000; // Time constant for engine-off; reflects the fact that oil is no longer getting circulated
- }
-
- float dOilTempdt = (target_oil_temp - oil_temp) / time_constant;
-
- oil_temp += (dOilTempdt * time_step);
-
- return (oil_temp);
-}
-
-// Calculate Oil Pressure
-float FGNewEngine::Calc_Oil_Press (float Oil_Temp, float Engine_RPM)
-{
- float Oil_Pressure = 0; //PSI
- float Oil_Press_Relief_Valve = 60; //PSI
- float Oil_Press_RPM_Max = 1800;
- float Design_Oil_Temp = 85; //Celsius
- float Oil_Viscosity_Index = 0.25; // PSI/Deg C
-// float Temp_Deviation = 0; // Deg C
-
- Oil_Pressure = (Oil_Press_Relief_Valve / Oil_Press_RPM_Max) * Engine_RPM;
-
- // Pressure relief valve opens at Oil_Press_Relief_Valve PSI setting
- if (Oil_Pressure >= Oil_Press_Relief_Valve) {
- Oil_Pressure = Oil_Press_Relief_Valve;
- }
-
- // Now adjust pressure according to Temp which affects the viscosity
-
- Oil_Pressure += (Design_Oil_Temp - Oil_Temp) * Oil_Viscosity_Index;
-
- return Oil_Pressure;
-}
-
-
-// Propeller calculations.
-void FGNewEngine::Do_Prop_Calcs()
-{
- float Gear_Ratio = 1.0;
- float forward_velocity; // m/s
- float prop_power_consumed_SI; // Watts
- double J; // advance ratio - dimensionless
- double Cp_20; // coefficient of power for 20 degree blade angle
- double Cp_25; // coefficient of power for 25 degree blade angle
- double Cp; // Our actual coefficient of power
- double neta_prop_20;
- double neta_prop_25;
- double neta_prop; // prop efficiency
-
- FGProp1_RPS = RPM * Gear_Ratio / 60.0;
- angular_velocity_SI = 2.0 * LS_PI * RPM / 60.0;
- forward_velocity = IAS * 0.514444444444; // Convert to m/s
-
- if(FGProp1_RPS == 0)
- J = 0;
- else
- J = forward_velocity / (FGProp1_RPS * prop_diameter);
- //cout << "advance_ratio = " << J << '\n';
-
- //Cp correlations based on data from McCormick
- Cp_20 = 0.0342*J*J*J*J - 0.1102*J*J*J + 0.0365*J*J - 0.0133*J + 0.064;
- Cp_25 = 0.0119*J*J*J*J - 0.0652*J*J*J + 0.018*J*J - 0.0077*J + 0.0921;
-
- //cout << "Cp_20 = " << Cp_20 << '\n';
- //cout << "Cp_25 = " << Cp_25 << '\n';
-
- //Assume that the blade angle is between 20 and 25 deg - reasonable for fixed pitch prop but won't hold for variable one !!!
- Cp = Cp_20 + ( (Cp_25 - Cp_20) * ((blade_angle - 20)/(25 - 20)) );
- //cout << "Cp = " << Cp << '\n';
- //cout << "RPM = " << RPM << '\n';
- //cout << "angular_velocity_SI = " << angular_velocity_SI << '\n';
-
- prop_power_consumed_SI = Cp * rho_air * pow(FGProp1_RPS,3.0f) * pow(float(prop_diameter),5.0f);
- //cout << "prop HP consumed = " << prop_power_consumed_SI / 745.699 << '\n';
- if(angular_velocity_SI == 0)
- prop_torque = 0;
- // However this can give problems - if rpm == 0 but forward velocity increases the prop should be able to generate a torque to start the engine spinning
- // Unlikely to happen in practice - but I suppose someone could move the lever of a stopped large piston engine from feathered to windmilling.
- // This *does* give problems - if the plane is put into a steep climb whilst windmilling the engine friction will eventually stop it spinning.
- // When put back into a dive it never starts re-spinning again. Although it is unlikely that anyone would do this in real life, they might well do it in a sim!!!
- else
- prop_torque = prop_power_consumed_SI / angular_velocity_SI;
-
- // calculate neta_prop here
- neta_prop_20 = 0.1328*J*J*J*J - 1.3073*J*J*J + 0.3525*J*J + 1.5591*J + 0.0007;
- neta_prop_25 = -0.3121*J*J*J*J + 0.4234*J*J*J - 0.7686*J*J + 1.5237*J - 0.0004;
- neta_prop = neta_prop_20 + ( (neta_prop_25 - neta_prop_20) * ((blade_angle - 20)/(25 - 20)) );
-
- // Check for zero forward velocity to avoid divide by zero
- if(forward_velocity < 0.0001)
- prop_thrust = 0.0;
- // I don't see how this works - how can the plane possibly start from rest ???
- // Hmmmm - it works because the forward_velocity at present never drops below about 0.03 even at rest
- // We can't really rely on this in the future - needs fixing !!!!
- // The problem is that we're doing this calculation backwards - we're working out the thrust from the power consumed and the velocity, which becomes invalid as velocity goes to zero.
- // It would be far more natural to work out the power consumed from the thrust - FIXME!!!!!.
- else
- prop_thrust = neta_prop * prop_power_consumed_SI / forward_velocity; //TODO - rename forward_velocity to IAS_SI
- //cout << "prop_thrust = " << prop_thrust << '\n';
-}
+++ /dev/null
-// IO360.hxx - a piston engine model currently for the IO360 engine fitted to the C172
-// but with the potential to model other naturally aspirated piston engines
-// given appropriate config input.
-//
-// Written by David Luff, started 2000.
-// Based on code by Phil Schubert, started 1999.
-//
-// This program is free software; you can redistribute it and/or
-// modify it under the terms of the GNU General Public License as
-// published by the Free Software Foundation; either version 2 of the
-// License, or (at your option) any later version.
-//
-// This program is distributed in the hope that it will be useful, but
-// WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-// General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
-
-
-#ifndef _IO360_HXX_
-#define _IO360_HXX_
-
-#include <simgear/compiler.h>
-
-#include <math.h>
-
-#include STL_IOSTREAM
-#include STL_FSTREAM
-
-SG_USING_STD(ofstream);
-
-class FGNewEngine {
-
-private:
-
- // These constants should probably be moved eventually
- float CONVERT_CUBIC_INCHES_TO_METERS_CUBED;
- float CONVERT_HP_TO_WATTS;
-
- // Properties of working fluids
- float Cp_air; // J/KgK
- float Cp_fuel; // J/KgK
- float calorific_value_fuel; // W/Kg
- float rho_fuel; // kg/m^3
- float rho_air; // kg/m^3
-
- // environment inputs
- float p_amb_sea_level; // Sea level ambient pressure in Pascals
- float p_amb; // Ambient pressure at current altitude in Pascals
- float T_amb; // ditto deg Kelvin
-
- // Control inputs
- float Throttle_Lever_Pos; // 0 = Closed, 100 = Fully Open
- float Propeller_Lever_Pos; // 0 = Full Course 100 = Full Fine
- float Mixture_Lever_Pos; // 0 = Idle Cut Off 100 = Full Rich
- int mag_pos; // 0=off, 1=left, 2=right, 3=both.
- bool starter;
-
- //misc
- float IAS;
- double time_step;
-
- // Engine Specific Variables that should be read in from a config file
- float MaxHP; // Horsepower
- float displacement; // Cubic inches
- float displacement_SI; //m^3 (derived from above rather than read in)
- float engine_inertia; //kg.m^2
- float prop_inertia; //kg.m^2
- float Max_Fuel_Flow; // Units??? Do we need this variable any more??
-
- // Engine specific variables that maybe should be read in from config but are pretty generic and won't vary much for a naturally aspirated piston engine.
- float Max_Manifold_Pressure; // inches Hg - typical manifold pressure at full throttle and typical max rpm
- float Min_Manifold_Pressure; // inches Hg - typical manifold pressure at idle (closed throttle)
- float Max_RPM; // rpm - this is really a bit of a hack and could be make redundant if the prop model works properly and takes tips at speed of sound into account.
- float Min_RPM; // rpm - possibly redundant ???
- float Mag_Derate_Percent; // Percentage reduction in power when mags are switched from 'both' to either 'L' or 'R'
- float Gear_Ratio; // Gearing between engine and propellor
- float n_R; // Number of cycles per power stroke - 2 for a 4 stroke engine.
-
- // Engine Variables not read in from config file
- float RPM; // rpm
- float Percentage_Power; // Power output as percentage of maximum power output
- float Manifold_Pressure; // Inches Hg
- float Fuel_Flow_gals_hr; // USgals/hour
- float Torque_lbft; // lb-ft
- float Torque_SI; // Nm
- float CHT; // Cylinder head temperature deg K
- float CHT_degF; // Ditto in deg Fahrenheit
- float Mixture;
- float Oil_Pressure; // PSI
- float Oil_Temp; // Deg C
- float current_oil_temp; // deg K
- /**** one of these is superfluous !!!!***/
- float HP; // Current power output in HP
- float Power_SI; // Current power output in Watts
- float RPS;
- float angular_velocity_SI; // rad/s
- float Torque_FMEP; // The component of Engine torque due to FMEP (Nm)
- float Torque_Imbalance; // difference between engine and prop torque
- float EGT; // Exhaust gas temperature deg K
- float EGT_degF; // Exhaust gas temperature deg Fahrenheit
- float volumetric_efficiency;
- float combustion_efficiency;
- float equivalence_ratio; // ratio of stoichiometric AFR over actual AFR
- float thi_sea_level; // the equivalence ratio we would be running at assuming sea level air denisity in the manifold
- float v_dot_air; // volume flow rate of air into engine - m^3/s
- float m_dot_air; // mass flow rate of air into engine - kg/s
- float m_dot_fuel; // mass flow rate of fuel into engine - kg/s
- float swept_volume; // total engine swept volume - m^3
- /********* swept volume or the geometry used to calculate it should be in the config read section surely ??? ******/
- float True_Manifold_Pressure; //in Hg - actual manifold gauge pressure
- float rho_air_manifold; // denisty of air in the manifold - kg/m^3
- float R_air; // Gas constant of air (287) UNITS???
- float delta_T_exhaust; // Temperature change of exhaust this time step - degK
- float heat_capacity_exhaust; // exhaust gas specific heat capacity - J/kgK
- float enthalpy_exhaust; // Enthalpy at current exhaust gas conditions - UNITS???
- float Percentage_of_best_power_mixture_power; // Current power as a percentage of what power we would have at the same conditions but at best power mixture
- float abstract_mixture; //temporary hack
- float angular_acceleration; //rad/s^2
- float FMEP; //Friction Mean Effective Pressure (Pa)
-
- // Various bits of housekeeping describing the engines state.
- bool running; // flag to indicate the engine is running self sustaining
- bool cranking; // flag to indicate the engine is being cranked
- int crank_counter; // Number of iterations that the engine has been cranked non-stop
- bool spark; // flag to indicate a spark is available
- bool fuel; // flag to indicate fuel is available
-
- // Propellor Variables
- float FGProp1_RPS; // rps
- float prop_torque; // Nm
- float prop_thrust; // Newtons
- double prop_diameter; // meters
- double blade_angle; // degrees
-
-
-// MEMBER FUNCTIONS
-
- // Calculate Engine RPM based on Propellor Lever Position
- float Calc_Engine_RPM(float Position);
-
- // Calculate Manifold Pressure based on throttle lever position
- // Note that this is simplistic and needs altering to include engine speed effects
- float Calc_Manifold_Pressure( float LeverPosn, float MaxMan, float MinMan);
-
- // Calculate combustion efficiency based on equivalence ratio
- float Lookup_Combustion_Efficiency(float thi_actual);
-
- // Calculate percentage of best power mixture power based on equivalence ratio
- float Power_Mixture_Correlation(float thi_actual);
-
- // Calculate exhaust gas temperature change
- float Calculate_Delta_T_Exhaust(void);
-
- // Calculate cylinder head temperature
- void Calc_CHT(void);
-
- // Calculate exhaust gas temperature
- void Calc_EGT(void);
-
- // Calculate fuel flow in gals/hr
- void Calc_Fuel_Flow_Gals_Hr(void);
-
- // Calculate current percentage power
- void Calc_Percentage_Power(bool mag_left, bool mag_right);
-
- // Calculate Oil Temperature
- float Calc_Oil_Temp (float oil_temp);
-
- // Calculate Oil Pressure
- float Calc_Oil_Press (float Oil_Temp, float Engine_RPM);
-
- // Propeller calculations.
- void Do_Prop_Calcs(void);
-
-public:
-
-// ofstream outfile;
-
- //constructor
- FGNewEngine() {
-// outfile.open("FGNewEngine.dat", ios::out|ios::trunc);
- }
-
- //destructor
- ~FGNewEngine() {
-// outfile.close();
- }
-
- // set initial default values
- void init(double dt);
-
- // update the engine model based on current control positions
- void update();
-
- inline void set_IAS( float value ) { IAS = value; }
- inline void set_Throttle_Lever_Pos( float value ) {
- Throttle_Lever_Pos = value;
- }
- inline void set_Propeller_Lever_Pos( float value ) {
- Propeller_Lever_Pos = value;
- }
- inline void set_Mixture_Lever_Pos( float value ) {
- Mixture_Lever_Pos = value;
- }
- // set the magneto switch position
- inline void set_Magneto_Switch_Pos( int value ) {
- mag_pos = value;
- }
- inline void setStarterFlag( bool flag ) {
- starter = flag;
- }
- // set ambient pressure - takes pounds per square foot
- inline void set_p_amb( float value ) {
- p_amb = value * 47.88026;
- // Convert to Pascals
- }
- // set ambient temperature - takes degrees Rankine
- inline void set_T_amb( float value ) {
- T_amb = value * 0.555555555556;
- // Convert to degrees Kelvin
- }
-
- // accessors
- inline float get_RPM() const { return RPM; }
- inline float get_Manifold_Pressure() const { return True_Manifold_Pressure; }
- // inline float get_Rho() const { return Rho; }
- inline float get_MaxHP() const { return MaxHP; }
- inline float get_Percentage_Power() const { return Percentage_Power; }
- inline float get_EGT() const { return EGT_degF; } // Returns EGT in Fahrenheit
- inline float get_CHT() const { return CHT_degF; } // Note this returns CHT in Fahrenheit
- inline float get_prop_thrust_SI() const { return prop_thrust; }
- inline float get_prop_thrust_lbs() const { return (prop_thrust * 0.2248); }
- inline float get_fuel_flow_gals_hr() const { return (Fuel_Flow_gals_hr); }
- inline float get_oil_temp() const { return ((current_oil_temp * 1.8) - 459.67); }
- inline bool getRunningFlag() const { return running; }
- inline bool getCrankingFlag() const { return cranking; }
- inline bool getStarterFlag() const { return starter; }
-};
-
-
-#endif // _IO360_HXX_
--- /dev/null
+// IO360.cxx - a piston engine model currently for the IO360 engine fitted to the C172
+// but with the potential to model other naturally aspirated piston engines
+// given appropriate config input.
+//
+// Written by David Luff, started 2000.
+// Based on code by Phil Schubert, started 1999.
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of the
+// License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+
+#include <simgear/compiler.h>
+
+#include <math.h>
+
+#include STL_FSTREAM
+#include STL_IOSTREAM
+
+#include <Main/fg_props.hxx>
+
+#include "IO360.hxx"
+#include "ls_constants.h"
+
+
+//*************************************************************************************
+// Initialise the engine model
+void FGNewEngine::init(double dt) {
+
+ // These constants should probably be moved eventually
+ CONVERT_CUBIC_INCHES_TO_METERS_CUBED = 1.638706e-5;
+ CONVERT_HP_TO_WATTS = 745.6999;
+
+ // Properties of working fluids
+ Cp_air = 1005; // J/KgK
+ Cp_fuel = 1700; // J/KgK
+ calorific_value_fuel = 47.3e6; // W/Kg Note that this is only an approximate value
+ rho_fuel = 800; // kg/m^3 - an estimate for now
+ R_air = 287.3;
+
+ // environment inputs
+ p_amb_sea_level = 101325; // Pascals
+
+ // Control inputs - ARE THESE NEEDED HERE???
+ Throttle_Lever_Pos = 75;
+ Propeller_Lever_Pos = 75;
+ Mixture_Lever_Pos = 100;
+
+ //misc
+ IAS = 0;
+ time_step = dt;
+
+ // Engine Specific Variables that should be read in from a config file
+ MaxHP = 200; //Lycoming IO360 -A-C-D series
+// MaxHP = 180; //Current Lycoming IO360 ?
+// displacement = 520; //Continental IO520-M
+ displacement = 360; //Lycoming IO360
+ displacement_SI = displacement * CONVERT_CUBIC_INCHES_TO_METERS_CUBED;
+ engine_inertia = 0.2; //kgm^2 - value taken from a popular family saloon car engine - need to find an aeroengine value !!!!!
+ prop_inertia = 0.05; //kgm^2 - this value is a total guess - dcl
+ Max_Fuel_Flow = 130; // Units??? Do we need this variable any more??
+
+ // Engine specific variables that maybe should be read in from config but are pretty generic and won't vary much for a naturally aspirated piston engine.
+ Max_Manifold_Pressure = 28.50; //Inches Hg. An approximation - should be able to find it in the engine performance data
+ Min_Manifold_Pressure = 6.5; //Inches Hg. This is a guess corresponding to approx 0.24 bar MAP (7 in Hg) - need to find some proper data for this
+ Max_RPM = 2700;
+ Min_RPM = 600; //Recommended idle from Continental data sheet
+ Mag_Derate_Percent = 5;
+ Gear_Ratio = 1;
+ n_R = 2; // Number of crank revolutions per power cycle - 2 for a 4 stroke engine.
+
+ // Various bits of housekeeping describing the engines initial state.
+ running = false;
+ cranking = false;
+ crank_counter = false;
+ starter = false;
+
+ // Initialise Engine Variables used by this instance
+ if(running)
+ RPM = 600;
+ else
+ RPM = 0;
+ Percentage_Power = 0;
+ Manifold_Pressure = 29.96; // Inches
+ Fuel_Flow_gals_hr = 0;
+// Torque = 0;
+ Torque_SI = 0;
+ CHT = 298.0; //deg Kelvin
+ CHT_degF = (CHT_degF * 1.8) - 459.67; //deg Fahrenheit
+ Mixture = 14;
+ Oil_Pressure = 0; // PSI
+ Oil_Temp = 85; // Deg C
+ current_oil_temp = 298.0; //deg Kelvin
+ /**** one of these is superfluous !!!!***/
+ HP = 0;
+ RPS = 0;
+ Torque_Imbalance = 0;
+
+ // Initialise Propellor Variables used by this instance
+ FGProp1_RPS = 0;
+ // Hardcode propellor for now - the following two should be read from config eventually
+ prop_diameter = 1.8; // meters
+ blade_angle = 23.0; // degrees
+}
+
+//*****************************************************************************
+// update the engine model based on current control positions
+void FGNewEngine::update() {
+
+/*
+ // Hack for testing - should output every 5 seconds
+ static int count1 = 0;
+ if(count1 == 0) {
+// cout << "P_atmos = " << p_amb << " T_atmos = " << T_amb << '\n';
+// cout << "Manifold pressure = " << Manifold_Pressure << " True_Manifold_Pressure = " << True_Manifold_Pressure << '\n';
+// cout << "p_amb_sea_level = " << p_amb_sea_level << '\n';
+// cout << "equivalence_ratio = " << equivalence_ratio << '\n';
+// cout << "combustion_efficiency = " << combustion_efficiency << '\n';
+// cout << "AFR = " << 14.7 / equivalence_ratio << '\n';
+// cout << "Mixture lever = " << Mixture_Lever_Pos << '\n';
+// cout << "n = " << RPM << " rpm\n";
+// cout << "T_amb = " << T_amb << '\n';
+// cout << "running = " << running << '\n';
+// cout << "fuel = " << fgGetFloat("/consumables/fuel/tank[0]/level-gal_us") << '\n';
+// cout << "Percentage_Power = " << Percentage_Power << '\n';
+// cout << "current_oil_temp = " << current_oil_temp << '\n';
+// cout << "EGT = " << EGT << '\n';
+ }
+ count1++;
+ if(count1 == 100)
+ count1 = 0;
+*/
+
+ // Check parameters that may alter the operating state of the engine.
+ // (spark, fuel, starter motor etc)
+
+ // Check for spark
+ bool Magneto_Left = false;
+ bool Magneto_Right = false;
+ // Magneto positions:
+ // 0 -> off
+ // 1 -> left only
+ // 2 -> right only
+ // 3 -> both
+ if(mag_pos != 0) {
+ spark = true;
+ } else {
+ spark = false;
+ } // neglects battery voltage, master on switch, etc for now.
+ if((mag_pos == 1) || (mag_pos > 2))
+ Magneto_Left = true;
+ if(mag_pos > 1)
+ Magneto_Right = true;
+
+ // crude check for fuel
+ if((fgGetFloat("/consumables/fuel/tank[0]/level-gal_us") > 0) || (fgGetFloat("/consumables/fuel/tank[1]/level-gal_us") > 0)) {
+ fuel = true;
+ } else {
+ fuel = false;
+ } // Need to make this better, eg position of fuel selector switch.
+
+ // Check if we are turning the starter motor
+ if(cranking != starter) {
+ // This check saves .../cranking from getting updated every loop - they only update when changed.
+ cranking = starter;
+ crank_counter = 0;
+ }
+ // Note that although /engines/engine[0]/starter and /engines/engine[0]/cranking might appear to be duplication it is
+ // not since the starter may be engaged with the battery voltage too low for cranking to occur (or perhaps the master
+ // switch just left off) and the sound manager will read .../cranking to determine wether to play a cranking sound.
+ // For now though none of that is implemented so cranking can be set equal to .../starter without further checks.
+
+// int Alternate_Air_Pos =0; // Off = 0. Reduces power by 3 % for same throttle setting
+ // DCL - don't know what this Alternate_Air_Pos is - this is a leftover from the Schubert code.
+
+ //Check mode of engine operation
+ if(cranking) {
+ crank_counter++;
+ if(RPM <= 480) {
+ RPM += 100;
+ if(RPM > 480)
+ RPM = 480;
+ } else {
+ // consider making a horrible noise if the starter is engaged with the engine running
+ }
+ }
+ if((!running) && (spark) && (fuel) && (crank_counter > 120)) {
+ // start the engine if revs high enough
+ if(RPM > 450) {
+ // For now just instantaneously start but later we should maybe crank for a bit
+ running = true;
+// RPM = 600;
+ }
+ }
+ if( (running) && ((!spark)||(!fuel)) ) {
+ // Cut the engine
+ // note that we only cut the power - the engine may continue to spin if the prop is in a moving airstream
+ running = false;
+ }
+
+ // Now we've ascertained whether the engine is running or not we can start to do the engine calculations 'proper'
+
+ // Calculate Sea Level Manifold Pressure
+ Manifold_Pressure = Calc_Manifold_Pressure( Throttle_Lever_Pos, Max_Manifold_Pressure, Min_Manifold_Pressure );
+ // cout << "manifold pressure = " << Manifold_Pressure << endl;
+
+ //Then find the actual manifold pressure (the calculated one is the sea level pressure)
+ True_Manifold_Pressure = Manifold_Pressure * p_amb / p_amb_sea_level;
+
+ //Do the fuel flow calculations
+ Calc_Fuel_Flow_Gals_Hr();
+
+ //Calculate engine power
+ Calc_Percentage_Power(Magneto_Left, Magneto_Right);
+ HP = Percentage_Power * MaxHP / 100.0;
+ Power_SI = HP * CONVERT_HP_TO_WATTS;
+
+ // FMEP calculation. For now we will just use this during motored operation.
+ // Eventually we will calculate IMEP and use the FMEP all the time to give BMEP (maybe!)
+ if(!running) {
+ // This FMEP data is from the Patton paper, assumes fully warm conditions.
+ FMEP = 1e-12*pow(RPM,4) - 1e-8*pow(RPM,3) + 5e-5*pow(RPM,2) - 0.0722*RPM + 154.85;
+ // Gives FMEP in kPa - now convert to Pa
+ FMEP *= 1000;
+ } else {
+ FMEP = 0.0;
+ }
+ // Is this total FMEP or friction FMEP ???
+
+ Torque_FMEP = (FMEP * displacement_SI) / (2.0 * LS_PI * n_R);
+
+ // Calculate Engine Torque. Check for div by zero since percentage power correlation does not guarantee zero power at zero rpm.
+ // However this is problematical since there is a resistance to movement even at rest
+ // Ie this is a dynamics equation not a statics one. This can be solved by going over to MEP based torque calculations.
+ if(RPM == 0) {
+ Torque_SI = 0 - Torque_FMEP;
+ }
+ else {
+ Torque_SI = ((Power_SI * 60.0) / (2.0 * LS_PI * RPM)) - Torque_FMEP; //Torque = power / angular velocity
+ // cout << Torque << " Nm\n";
+ }
+
+ //Calculate Exhaust gas temperature
+ if(running)
+ Calc_EGT();
+ else
+ EGT = 298.0;
+
+ // Calculate Cylinder Head Temperature
+ Calc_CHT();
+
+ // Calculate oil temperature
+ current_oil_temp = Calc_Oil_Temp(current_oil_temp);
+
+ // Calculate Oil Pressure
+ Oil_Pressure = Calc_Oil_Press( Oil_Temp, RPM );
+
+ // Now do the Propeller Calculations
+ Do_Prop_Calcs();
+
+// Now do the engine - prop torque balance to calculate final RPM
+
+ //Calculate new RPM from torque balance and inertia.
+ Torque_Imbalance = Torque_SI - prop_torque; //This gives a +ve value when the engine torque exeeds the prop torque
+ // (Engine torque is +ve when it acts in the direction of engine revolution, prop torque is +ve when it opposes the direction of engine revolution)
+
+ angular_acceleration = Torque_Imbalance / (engine_inertia + prop_inertia);
+ angular_velocity_SI += (angular_acceleration * time_step);
+ // Don't let the engine go into reverse
+ if(angular_velocity_SI < 0)
+ angular_velocity_SI = 0;
+ RPM = (angular_velocity_SI * 60) / (2.0 * LS_PI);
+
+ // And finally a last check on the engine state after doing the torque balance with the prop - have we stalled?
+ if(running) {
+ //Check if we have stalled the engine
+ if (RPM == 0) {
+ running = false;
+ } else if((RPM <= 480) && (cranking)) {
+ //Make sure the engine noise dosn't play if the engine won't start due to eg mixture lever pulled out.
+ running = false;
+ EGT = 298.0;
+ }
+ }
+
+ // And finally, do any unit conversions from internal units to output units
+ EGT_degF = (EGT * 1.8) - 459.67;
+ CHT_degF = (CHT * 1.8) - 459.67;
+}
+
+//*****************************************************************************************************
+
+// FGNewEngine member functions
+
+////////////////////////////////////////////////////////////////////////////////////////////
+// Return the combustion efficiency as a function of equivalence ratio
+//
+// Combustion efficiency values from Heywood,
+// "Internal Combustion Engine Fundamentals", ISBN 0-07-100499-8
+////////////////////////////////////////////////////////////////////////////////////////////
+float FGNewEngine::Lookup_Combustion_Efficiency(float thi_actual)
+{
+ const int NUM_ELEMENTS = 11;
+ float thi[NUM_ELEMENTS] = {0.0, 0.9, 1.0, 1.05, 1.1, 1.15, 1.2, 1.3, 1.4, 1.5, 1.6}; //array of equivalence ratio values
+ float neta_comb[NUM_ELEMENTS] = {0.98, 0.98, 0.97, 0.95, 0.9, 0.85, 0.79, 0.7, 0.63, 0.57, 0.525}; //corresponding array of combustion efficiency values
+ float neta_comb_actual = 0.0f;
+ float factor;
+
+ int i;
+ int j = NUM_ELEMENTS; //This must be equal to the number of elements in the lookup table arrays
+
+ for(i=0;i<j;i++)
+ {
+ if(i == (j-1)) {
+ // Assume linear extrapolation of the slope between the last two points beyond the last point
+ float dydx = (neta_comb[i] - neta_comb[i-1]) / (thi[i] - thi[i-1]);
+ neta_comb_actual = neta_comb[i] + dydx * (thi_actual - thi[i]);
+ return neta_comb_actual;
+ }
+ if(thi_actual == thi[i]) {
+ neta_comb_actual = neta_comb[i];
+ return neta_comb_actual;
+ }
+ if((thi_actual > thi[i]) && (thi_actual < thi[i + 1])) {
+ //do linear interpolation between the two points
+ factor = (thi_actual - thi[i]) / (thi[i+1] - thi[i]);
+ neta_comb_actual = (factor * (neta_comb[i+1] - neta_comb[i])) + neta_comb[i];
+ return neta_comb_actual;
+ }
+ }
+
+ //if we get here something has gone badly wrong
+ SG_LOG(SG_AIRCRAFT, SG_ALERT, "error in FGNewEngine::Lookup_Combustion_Efficiency");
+ return neta_comb_actual;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+// Return the percentage of best mixture power available at a given mixture strength
+//
+// Based on data from "Technical Considerations for Catalysts for the European Market"
+// by H S Gandi, published 1988 by IMechE
+//
+// Note that currently no attempt is made to set a lean limit on stable combustion
+////////////////////////////////////////////////////////////////////////////////////////////
+float FGNewEngine::Power_Mixture_Correlation(float thi_actual)
+{
+ float AFR_actual = 14.7 / thi_actual;
+ // thi and thi_actual are equivalence ratio
+ const int NUM_ELEMENTS = 13;
+ // The lookup table is in AFR because the source data was. I added the two end elements to make sure we are almost always in it.
+ float AFR[NUM_ELEMENTS] = {(14.7/1.6), 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, (14.7/0.6)}; //array of equivalence ratio values
+ float mixPerPow[NUM_ELEMENTS] = {78, 86, 93.5, 98, 100, 99, 96.4, 92.5, 88, 83, 78.5, 74, 58}; //corresponding array of combustion efficiency values
+ float mixPerPow_actual = 0.0f;
+ float factor;
+ float dydx;
+
+ int i;
+ int j = NUM_ELEMENTS; //This must be equal to the number of elements in the lookup table arrays
+
+ for(i=0;i<j;i++)
+ {
+ if(i == (j-1)) {
+ // Assume linear extrapolation of the slope between the last two points beyond the last point
+ dydx = (mixPerPow[i] - mixPerPow[i-1]) / (AFR[i] - AFR[i-1]);
+ mixPerPow_actual = mixPerPow[i] + dydx * (AFR_actual - AFR[i]);
+ return mixPerPow_actual;
+ }
+ if((i == 0) && (AFR_actual < AFR[i])) {
+ // Assume linear extrapolation of the slope between the first two points for points before the first point
+ dydx = (mixPerPow[i] - mixPerPow[i-1]) / (AFR[i] - AFR[i-1]);
+ mixPerPow_actual = mixPerPow[i] + dydx * (AFR_actual - AFR[i]);
+ return mixPerPow_actual;
+ }
+ if(AFR_actual == AFR[i]) {
+ mixPerPow_actual = mixPerPow[i];
+ return mixPerPow_actual;
+ }
+ if((AFR_actual > AFR[i]) && (AFR_actual < AFR[i + 1])) {
+ //do linear interpolation between the two points
+ factor = (AFR_actual - AFR[i]) / (AFR[i+1] - AFR[i]);
+ mixPerPow_actual = (factor * (mixPerPow[i+1] - mixPerPow[i])) + mixPerPow[i];
+ return mixPerPow_actual;
+ }
+ }
+
+ //if we get here something has gone badly wrong
+ SG_LOG(SG_AIRCRAFT, SG_ALERT, "error in FGNewEngine::Power_Mixture_Correlation");
+ return mixPerPow_actual;
+}
+
+// Calculate Cylinder Head Temperature
+// Crudely models the cylinder head as an arbitary lump of arbitary size and area with one third of combustion energy
+// as heat input and heat output as a function of airspeed and temperature. Could be improved!!!
+void FGNewEngine::Calc_CHT()
+{
+ float h1 = -95.0; //co-efficient for free convection
+ float h2 = -3.95; //co-efficient for forced convection
+ float h3 = -0.05; //co-efficient for forced convection due to prop backwash
+ float v_apparent; //air velocity over cylinder head in m/s
+ float v_dot_cooling_air;
+ float m_dot_cooling_air;
+ float temperature_difference;
+ float arbitary_area = 1.0;
+ float dqdt_from_combustion;
+ float dqdt_forced; //Rate of energy transfer to/from cylinder head due to forced convection (Joules) (sign convention: to cylinder head is +ve)
+ float dqdt_free; //Rate of energy transfer to/from cylinder head due to free convection (Joules) (sign convention: to cylinder head is +ve)
+ float dqdt_cylinder_head; //Overall energy change in cylinder head
+ float CpCylinderHead = 800.0; //FIXME - this is a guess - I need to look up the correct value
+ float MassCylinderHead = 8.0; //Kg - this is a guess - it dosn't have to be absolutely accurate but can be varied to alter the warm-up rate
+ float HeatCapacityCylinderHead;
+ float dCHTdt;
+
+ // The above values are hardwired to give reasonable results for an IO360 (C172 engine)
+ // Now adjust to try to compensate for arbitary sized engines - this is currently untested
+ arbitary_area *= (MaxHP / 180.0);
+ MassCylinderHead *= (MaxHP / 180.0);
+
+ temperature_difference = CHT - T_amb;
+
+ v_apparent = IAS * 0.5144444; //convert from knots to m/s
+ v_dot_cooling_air = arbitary_area * v_apparent;
+ m_dot_cooling_air = v_dot_cooling_air * rho_air;
+
+ //Calculate rate of energy transfer to cylinder head from combustion
+ dqdt_from_combustion = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
+
+ //Calculate rate of energy transfer from cylinder head due to cooling NOTE is calculated as rate to but negative
+ dqdt_forced = (h2 * m_dot_cooling_air * temperature_difference) + (h3 * RPM * temperature_difference);
+ dqdt_free = h1 * temperature_difference;
+
+ //Calculate net rate of energy transfer to or from cylinder head
+ dqdt_cylinder_head = dqdt_from_combustion + dqdt_forced + dqdt_free;
+
+ HeatCapacityCylinderHead = CpCylinderHead * MassCylinderHead;
+
+ dCHTdt = dqdt_cylinder_head / HeatCapacityCylinderHead;
+
+ CHT += (dCHTdt * time_step);
+}
+
+// Calculate exhaust gas temperature
+void FGNewEngine::Calc_EGT()
+{
+ combustion_efficiency = Lookup_Combustion_Efficiency(equivalence_ratio); //The combustion efficiency basically tells us what proportion of the fuels calorific value is released
+
+ //now calculate energy release to exhaust
+ //We will assume a three way split of fuel energy between useful work, the coolant system and the exhaust system
+ //This is a reasonable first suck of the thumb estimate for a water cooled automotive engine - whether it holds for an air cooled aero engine is probably open to question
+ //Regardless - it won't affect the variation of EGT with mixture, and we can always put a multiplier on EGT to get a reasonable peak value.
+ enthalpy_exhaust = m_dot_fuel * calorific_value_fuel * combustion_efficiency * 0.33;
+ heat_capacity_exhaust = (Cp_air * m_dot_air) + (Cp_fuel * m_dot_fuel);
+ delta_T_exhaust = enthalpy_exhaust / heat_capacity_exhaust;
+
+ EGT = T_amb + delta_T_exhaust;
+
+ //The above gives the exhaust temperature immediately prior to leaving the combustion chamber
+ //Now derate to give a more realistic figure as measured downstream
+ //For now we will aim for a peak of around 400 degC (750 degF)
+
+ EGT *= 0.444 + ((0.544 - 0.444) * Percentage_Power / 100.0);
+}
+
+// Calculate Manifold Pressure based on Throttle lever Position
+float FGNewEngine::Calc_Manifold_Pressure ( float LeverPosn, float MaxMan, float MinMan)
+{
+ float Inches;
+
+ //Note that setting the manifold pressure as a function of lever position only is not strictly accurate
+ //MAP is also a function of engine speed. (and ambient pressure if we are going for an actual MAP model)
+ Inches = MinMan + (LeverPosn * (MaxMan - MinMan) / 100);
+
+ //allow for idle bypass valve or slightly open throttle stop
+ if(Inches < MinMan)
+ Inches = MinMan;
+
+ //Check for stopped engine (crudest way of compensating for engine speed)
+ if(RPM == 0)
+ Inches = 29.92;
+
+ return Inches;
+}
+
+// Calculate fuel flow in gals/hr
+void FGNewEngine::Calc_Fuel_Flow_Gals_Hr()
+{
+ //DCL - calculate mass air flow into engine based on speed and load - separate this out into a function eventually
+ //t_amb is actual temperature calculated from altitude
+ //calculate density from ideal gas equation
+ rho_air = p_amb / ( R_air * T_amb );
+ rho_air_manifold = rho_air * Manifold_Pressure / 29.6; //This is a bit of a roundabout way of calculating this but it works !! If we put manifold pressure into SI units we could do it simpler.
+ //calculate ideal engine volume inducted per second
+ swept_volume = (displacement_SI * (RPM / 60)) / 2; //This equation is only valid for a four stroke engine
+ //calculate volumetric efficiency - for now we will just use 0.8, but actually it is a function of engine speed and the exhaust to manifold pressure ratio
+ //Note that this is cylinder vol eff - the throttle drop is already accounted for in the MAP calculation
+ volumetric_efficiency = 0.8;
+ //Now use volumetric efficiency to calculate actual air volume inducted per second
+ v_dot_air = swept_volume * volumetric_efficiency;
+ //Now calculate mass flow rate of air into engine
+ m_dot_air = v_dot_air * rho_air_manifold;
+
+//**************
+
+ //DCL - now calculate fuel flow into engine based on air flow and mixture lever position
+ //assume lever runs from no flow at fully out to thi = 1.3 at fully in at sea level
+ //also assume that the injector linkage is ideal - hence the set mixture is maintained at a given altitude throughout the speed and load range
+ thi_sea_level = 1.3 * ( Mixture_Lever_Pos / 100.0 );
+ equivalence_ratio = thi_sea_level * p_amb_sea_level / p_amb; //ie as we go higher the mixture gets richer for a given lever position
+ m_dot_fuel = m_dot_air / 14.7 * equivalence_ratio;
+ Fuel_Flow_gals_hr = (m_dot_fuel / rho_fuel) * 264.172 * 3600.0; // Note this assumes US gallons
+}
+
+// Calculate the percentage of maximum rated power delivered as a function of Manifold pressure multiplied by engine speed (rpm)
+// This is not necessarilly the best approach but seems to work for now.
+// May well need tweaking at the bottom end if the prop model is changed.
+void FGNewEngine::Calc_Percentage_Power(bool mag_left, bool mag_right)
+{
+ // For a given Manifold Pressure and RPM calculate the % Power
+ // Multiply Manifold Pressure by RPM
+ float ManXRPM = True_Manifold_Pressure * RPM;
+
+/*
+// Phil's %power correlation
+ // Calculate % Power
+ Percentage_Power = (+ 7E-09 * ManXRPM * ManXRPM) + ( + 7E-04 * ManXRPM) - 0.1218;
+ // cout << Percentage_Power << "%" << "\t";
+*/
+
+// DCL %power correlation - basically Phil's correlation modified to give slighty less power at the low end
+// might need some adjustment as the prop model is adjusted
+// My aim is to match the prop model and engine model at the low end to give the manufacturer's recommended idle speed with the throttle closed - 600rpm for the Continental IO520
+ // Calculate % Power for Nev's prop model
+ //Percentage_Power = (+ 6E-09 * ManXRPM * ManXRPM) + ( + 8E-04 * ManXRPM) - 1.8524;
+
+ // Calculate %power for DCL prop model
+ Percentage_Power = (7e-9 * ManXRPM * ManXRPM) + (7e-4 * ManXRPM) - 1.0;
+
+ // Power de-rating for altitude has been removed now that we are basing the power
+ // on the actual manifold pressure, which takes air pressure into account. However - this fails to
+ // take the temperature into account - this is TODO.
+
+ // Adjust power for temperature - this is temporary until the power is done as a function of mass flow rate induced
+ // Adjust for Temperature - Temperature above Standard decrease
+ // power by 7/120 % per degree F increase, and incease power for
+ // temps below at the same ratio
+ float T_amb_degF = (T_amb * 1.8) - 459.67;
+ float T_amb_sea_lev_degF = (288 * 1.8) - 459.67;
+ Percentage_Power = Percentage_Power + ((T_amb_sea_lev_degF - T_amb_degF) * 7 /120);
+
+ //DCL - now adjust power to compensate for mixture
+ Percentage_of_best_power_mixture_power = Power_Mixture_Correlation(equivalence_ratio);
+ Percentage_Power = Percentage_Power * Percentage_of_best_power_mixture_power / 100.0;
+
+ // Now Derate engine for the effects of Bad/Switched off magnetos
+ //if (Magneto_Left == 0 && Magneto_Right == 0) {
+ if(!running) {
+ // cout << "Both OFF\n";
+ Percentage_Power = 0;
+ } else if (mag_left && mag_right) {
+ // cout << "Both On ";
+ } else if (mag_left == 0 || mag_right== 0) {
+ // cout << "1 Magneto Failed ";
+ Percentage_Power = Percentage_Power * ((100.0 - Mag_Derate_Percent)/100.0);
+ // cout << FGEng1_Percentage_Power << "%" << "\t";
+ }
+/*
+ //DCL - stall the engine if RPM drops below 450 - this is possible if mixture lever is pulled right out
+ //This is a kludge that I should eliminate by adding a total fmep estimation.
+ if(RPM < 450)
+ Percentage_Power = 0;
+*/
+ if(Percentage_Power < 0)
+ Percentage_Power = 0;
+}
+
+// Calculate Oil Temperature in degrees Kelvin
+float FGNewEngine::Calc_Oil_Temp (float oil_temp)
+{
+ float idle_percentage_power = 2.3; // approximately
+ float target_oil_temp; // Steady state oil temp at the current engine conditions
+ float time_constant; // The time constant for the differential equation
+ if(running) {
+ target_oil_temp = 363;
+ time_constant = 500; // Time constant for engine-on idling.
+ if(Percentage_Power > idle_percentage_power) {
+ time_constant /= ((Percentage_Power / idle_percentage_power) / 10.0); // adjust for power
+ }
+ } else {
+ target_oil_temp = 298;
+ time_constant = 1000; // Time constant for engine-off; reflects the fact that oil is no longer getting circulated
+ }
+
+ float dOilTempdt = (target_oil_temp - oil_temp) / time_constant;
+
+ oil_temp += (dOilTempdt * time_step);
+
+ return (oil_temp);
+}
+
+// Calculate Oil Pressure
+float FGNewEngine::Calc_Oil_Press (float Oil_Temp, float Engine_RPM)
+{
+ float Oil_Pressure = 0; //PSI
+ float Oil_Press_Relief_Valve = 60; //PSI
+ float Oil_Press_RPM_Max = 1800;
+ float Design_Oil_Temp = 85; //Celsius
+ float Oil_Viscosity_Index = 0.25; // PSI/Deg C
+// float Temp_Deviation = 0; // Deg C
+
+ Oil_Pressure = (Oil_Press_Relief_Valve / Oil_Press_RPM_Max) * Engine_RPM;
+
+ // Pressure relief valve opens at Oil_Press_Relief_Valve PSI setting
+ if (Oil_Pressure >= Oil_Press_Relief_Valve) {
+ Oil_Pressure = Oil_Press_Relief_Valve;
+ }
+
+ // Now adjust pressure according to Temp which affects the viscosity
+
+ Oil_Pressure += (Design_Oil_Temp - Oil_Temp) * Oil_Viscosity_Index;
+
+ return Oil_Pressure;
+}
+
+
+// Propeller calculations.
+void FGNewEngine::Do_Prop_Calcs()
+{
+ float Gear_Ratio = 1.0;
+ float forward_velocity; // m/s
+ float prop_power_consumed_SI; // Watts
+ double J; // advance ratio - dimensionless
+ double Cp_20; // coefficient of power for 20 degree blade angle
+ double Cp_25; // coefficient of power for 25 degree blade angle
+ double Cp; // Our actual coefficient of power
+ double neta_prop_20;
+ double neta_prop_25;
+ double neta_prop; // prop efficiency
+
+ FGProp1_RPS = RPM * Gear_Ratio / 60.0;
+ angular_velocity_SI = 2.0 * LS_PI * RPM / 60.0;
+ forward_velocity = IAS * 0.514444444444; // Convert to m/s
+
+ if(FGProp1_RPS == 0)
+ J = 0;
+ else
+ J = forward_velocity / (FGProp1_RPS * prop_diameter);
+ //cout << "advance_ratio = " << J << '\n';
+
+ //Cp correlations based on data from McCormick
+ Cp_20 = 0.0342*J*J*J*J - 0.1102*J*J*J + 0.0365*J*J - 0.0133*J + 0.064;
+ Cp_25 = 0.0119*J*J*J*J - 0.0652*J*J*J + 0.018*J*J - 0.0077*J + 0.0921;
+
+ //cout << "Cp_20 = " << Cp_20 << '\n';
+ //cout << "Cp_25 = " << Cp_25 << '\n';
+
+ //Assume that the blade angle is between 20 and 25 deg - reasonable for fixed pitch prop but won't hold for variable one !!!
+ Cp = Cp_20 + ( (Cp_25 - Cp_20) * ((blade_angle - 20)/(25 - 20)) );
+ //cout << "Cp = " << Cp << '\n';
+ //cout << "RPM = " << RPM << '\n';
+ //cout << "angular_velocity_SI = " << angular_velocity_SI << '\n';
+
+ prop_power_consumed_SI = Cp * rho_air * pow(FGProp1_RPS,3.0f) * pow(float(prop_diameter),5.0f);
+ //cout << "prop HP consumed = " << prop_power_consumed_SI / 745.699 << '\n';
+ if(angular_velocity_SI == 0)
+ prop_torque = 0;
+ // However this can give problems - if rpm == 0 but forward velocity increases the prop should be able to generate a torque to start the engine spinning
+ // Unlikely to happen in practice - but I suppose someone could move the lever of a stopped large piston engine from feathered to windmilling.
+ // This *does* give problems - if the plane is put into a steep climb whilst windmilling the engine friction will eventually stop it spinning.
+ // When put back into a dive it never starts re-spinning again. Although it is unlikely that anyone would do this in real life, they might well do it in a sim!!!
+ else
+ prop_torque = prop_power_consumed_SI / angular_velocity_SI;
+
+ // calculate neta_prop here
+ neta_prop_20 = 0.1328*J*J*J*J - 1.3073*J*J*J + 0.3525*J*J + 1.5591*J + 0.0007;
+ neta_prop_25 = -0.3121*J*J*J*J + 0.4234*J*J*J - 0.7686*J*J + 1.5237*J - 0.0004;
+ neta_prop = neta_prop_20 + ( (neta_prop_25 - neta_prop_20) * ((blade_angle - 20)/(25 - 20)) );
+
+ // Check for zero forward velocity to avoid divide by zero
+ if(forward_velocity < 0.0001)
+ prop_thrust = 0.0;
+ // I don't see how this works - how can the plane possibly start from rest ???
+ // Hmmmm - it works because the forward_velocity at present never drops below about 0.03 even at rest
+ // We can't really rely on this in the future - needs fixing !!!!
+ // The problem is that we're doing this calculation backwards - we're working out the thrust from the power consumed and the velocity, which becomes invalid as velocity goes to zero.
+ // It would be far more natural to work out the power consumed from the thrust - FIXME!!!!!.
+ else
+ prop_thrust = neta_prop * prop_power_consumed_SI / forward_velocity; //TODO - rename forward_velocity to IAS_SI
+ //cout << "prop_thrust = " << prop_thrust << '\n';
+}
--- /dev/null
+// IO360.hxx - a piston engine model currently for the IO360 engine fitted to the C172
+// but with the potential to model other naturally aspirated piston engines
+// given appropriate config input.
+//
+// Written by David Luff, started 2000.
+// Based on code by Phil Schubert, started 1999.
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of the
+// License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+
+
+#ifndef _IO360_HXX_
+#define _IO360_HXX_
+
+#include <simgear/compiler.h>
+
+#include <math.h>
+
+#include STL_IOSTREAM
+#include STL_FSTREAM
+
+SG_USING_STD(ofstream);
+
+class FGNewEngine {
+
+private:
+
+ // These constants should probably be moved eventually
+ float CONVERT_CUBIC_INCHES_TO_METERS_CUBED;
+ float CONVERT_HP_TO_WATTS;
+
+ // Properties of working fluids
+ float Cp_air; // J/KgK
+ float Cp_fuel; // J/KgK
+ float calorific_value_fuel; // W/Kg
+ float rho_fuel; // kg/m^3
+ float rho_air; // kg/m^3
+
+ // environment inputs
+ float p_amb_sea_level; // Sea level ambient pressure in Pascals
+ float p_amb; // Ambient pressure at current altitude in Pascals
+ float T_amb; // ditto deg Kelvin
+
+ // Control inputs
+ float Throttle_Lever_Pos; // 0 = Closed, 100 = Fully Open
+ float Propeller_Lever_Pos; // 0 = Full Course 100 = Full Fine
+ float Mixture_Lever_Pos; // 0 = Idle Cut Off 100 = Full Rich
+ int mag_pos; // 0=off, 1=left, 2=right, 3=both.
+ bool starter;
+
+ //misc
+ float IAS;
+ double time_step;
+
+ // Engine Specific Variables that should be read in from a config file
+ float MaxHP; // Horsepower
+ float displacement; // Cubic inches
+ float displacement_SI; //m^3 (derived from above rather than read in)
+ float engine_inertia; //kg.m^2
+ float prop_inertia; //kg.m^2
+ float Max_Fuel_Flow; // Units??? Do we need this variable any more??
+
+ // Engine specific variables that maybe should be read in from config but are pretty generic and won't vary much for a naturally aspirated piston engine.
+ float Max_Manifold_Pressure; // inches Hg - typical manifold pressure at full throttle and typical max rpm
+ float Min_Manifold_Pressure; // inches Hg - typical manifold pressure at idle (closed throttle)
+ float Max_RPM; // rpm - this is really a bit of a hack and could be make redundant if the prop model works properly and takes tips at speed of sound into account.
+ float Min_RPM; // rpm - possibly redundant ???
+ float Mag_Derate_Percent; // Percentage reduction in power when mags are switched from 'both' to either 'L' or 'R'
+ float Gear_Ratio; // Gearing between engine and propellor
+ float n_R; // Number of cycles per power stroke - 2 for a 4 stroke engine.
+
+ // Engine Variables not read in from config file
+ float RPM; // rpm
+ float Percentage_Power; // Power output as percentage of maximum power output
+ float Manifold_Pressure; // Inches Hg
+ float Fuel_Flow_gals_hr; // USgals/hour
+ float Torque_lbft; // lb-ft
+ float Torque_SI; // Nm
+ float CHT; // Cylinder head temperature deg K
+ float CHT_degF; // Ditto in deg Fahrenheit
+ float Mixture;
+ float Oil_Pressure; // PSI
+ float Oil_Temp; // Deg C
+ float current_oil_temp; // deg K
+ /**** one of these is superfluous !!!!***/
+ float HP; // Current power output in HP
+ float Power_SI; // Current power output in Watts
+ float RPS;
+ float angular_velocity_SI; // rad/s
+ float Torque_FMEP; // The component of Engine torque due to FMEP (Nm)
+ float Torque_Imbalance; // difference between engine and prop torque
+ float EGT; // Exhaust gas temperature deg K
+ float EGT_degF; // Exhaust gas temperature deg Fahrenheit
+ float volumetric_efficiency;
+ float combustion_efficiency;
+ float equivalence_ratio; // ratio of stoichiometric AFR over actual AFR
+ float thi_sea_level; // the equivalence ratio we would be running at assuming sea level air denisity in the manifold
+ float v_dot_air; // volume flow rate of air into engine - m^3/s
+ float m_dot_air; // mass flow rate of air into engine - kg/s
+ float m_dot_fuel; // mass flow rate of fuel into engine - kg/s
+ float swept_volume; // total engine swept volume - m^3
+ /********* swept volume or the geometry used to calculate it should be in the config read section surely ??? ******/
+ float True_Manifold_Pressure; //in Hg - actual manifold gauge pressure
+ float rho_air_manifold; // denisty of air in the manifold - kg/m^3
+ float R_air; // Gas constant of air (287) UNITS???
+ float delta_T_exhaust; // Temperature change of exhaust this time step - degK
+ float heat_capacity_exhaust; // exhaust gas specific heat capacity - J/kgK
+ float enthalpy_exhaust; // Enthalpy at current exhaust gas conditions - UNITS???
+ float Percentage_of_best_power_mixture_power; // Current power as a percentage of what power we would have at the same conditions but at best power mixture
+ float abstract_mixture; //temporary hack
+ float angular_acceleration; //rad/s^2
+ float FMEP; //Friction Mean Effective Pressure (Pa)
+
+ // Various bits of housekeeping describing the engines state.
+ bool running; // flag to indicate the engine is running self sustaining
+ bool cranking; // flag to indicate the engine is being cranked
+ int crank_counter; // Number of iterations that the engine has been cranked non-stop
+ bool spark; // flag to indicate a spark is available
+ bool fuel; // flag to indicate fuel is available
+
+ // Propellor Variables
+ float FGProp1_RPS; // rps
+ float prop_torque; // Nm
+ float prop_thrust; // Newtons
+ double prop_diameter; // meters
+ double blade_angle; // degrees
+
+
+// MEMBER FUNCTIONS
+
+ // Calculate Engine RPM based on Propellor Lever Position
+ float Calc_Engine_RPM(float Position);
+
+ // Calculate Manifold Pressure based on throttle lever position
+ // Note that this is simplistic and needs altering to include engine speed effects
+ float Calc_Manifold_Pressure( float LeverPosn, float MaxMan, float MinMan);
+
+ // Calculate combustion efficiency based on equivalence ratio
+ float Lookup_Combustion_Efficiency(float thi_actual);
+
+ // Calculate percentage of best power mixture power based on equivalence ratio
+ float Power_Mixture_Correlation(float thi_actual);
+
+ // Calculate exhaust gas temperature change
+ float Calculate_Delta_T_Exhaust(void);
+
+ // Calculate cylinder head temperature
+ void Calc_CHT(void);
+
+ // Calculate exhaust gas temperature
+ void Calc_EGT(void);
+
+ // Calculate fuel flow in gals/hr
+ void Calc_Fuel_Flow_Gals_Hr(void);
+
+ // Calculate current percentage power
+ void Calc_Percentage_Power(bool mag_left, bool mag_right);
+
+ // Calculate Oil Temperature
+ float Calc_Oil_Temp (float oil_temp);
+
+ // Calculate Oil Pressure
+ float Calc_Oil_Press (float Oil_Temp, float Engine_RPM);
+
+ // Propeller calculations.
+ void Do_Prop_Calcs(void);
+
+public:
+
+// ofstream outfile;
+
+ //constructor
+ FGNewEngine() {
+// outfile.open("FGNewEngine.dat", ios::out|ios::trunc);
+ }
+
+ //destructor
+ ~FGNewEngine() {
+// outfile.close();
+ }
+
+ // set initial default values
+ void init(double dt);
+
+ // update the engine model based on current control positions
+ void update();
+
+ inline void set_IAS( float value ) { IAS = value; }
+ inline void set_Throttle_Lever_Pos( float value ) {
+ Throttle_Lever_Pos = value;
+ }
+ inline void set_Propeller_Lever_Pos( float value ) {
+ Propeller_Lever_Pos = value;
+ }
+ inline void set_Mixture_Lever_Pos( float value ) {
+ Mixture_Lever_Pos = value;
+ }
+ // set the magneto switch position
+ inline void set_Magneto_Switch_Pos( int value ) {
+ mag_pos = value;
+ }
+ inline void setStarterFlag( bool flag ) {
+ starter = flag;
+ }
+ // set ambient pressure - takes pounds per square foot
+ inline void set_p_amb( float value ) {
+ p_amb = value * 47.88026;
+ // Convert to Pascals
+ }
+ // set ambient temperature - takes degrees Rankine
+ inline void set_T_amb( float value ) {
+ T_amb = value * 0.555555555556;
+ // Convert to degrees Kelvin
+ }
+
+ // accessors
+ inline float get_RPM() const { return RPM; }
+ inline float get_Manifold_Pressure() const { return True_Manifold_Pressure; }
+ // inline float get_Rho() const { return Rho; }
+ inline float get_MaxHP() const { return MaxHP; }
+ inline float get_Percentage_Power() const { return Percentage_Power; }
+ inline float get_EGT() const { return EGT_degF; } // Returns EGT in Fahrenheit
+ inline float get_CHT() const { return CHT_degF; } // Note this returns CHT in Fahrenheit
+ inline float get_prop_thrust_SI() const { return prop_thrust; }
+ inline float get_prop_thrust_lbs() const { return (prop_thrust * 0.2248); }
+ inline float get_fuel_flow_gals_hr() const { return (Fuel_Flow_gals_hr); }
+ inline float get_oil_temp() const { return ((current_oil_temp * 1.8) - 459.67); }
+ inline bool getRunningFlag() const { return running; }
+ inline bool getCrankingFlag() const { return cranking; }
+ inline bool getStarterFlag() const { return starter; }
+};
+
+
+#endif // _IO360_HXX_
#include <Aircraft/aircraft.hxx>
#include <Controls/controls.hxx>
#include <FDM/flight.hxx>
-#include <FDM/IO360.hxx>
-#include <FDM/LaRCsimIC.hxx>
#include <FDM/UIUCModel/uiuc_aircraft.h>
#include <Main/fg_props.hxx>
#include <Model/acmodel.hxx>
#include "ls_interface.h"
#include "ls_constants.h"
+#include "IO360.hxx"
+#include "LaRCsimIC.hxx"
#include "LaRCsim.hxx"
#define _LARCSIM_HXX
-#include <FDM/IO360.hxx>
#include <FDM/flight.hxx>
-#include <FDM/LaRCsimIC.hxx>
+
+#include "IO360.hxx"
+#include "LaRCsimIC.hxx"
class FGLaRCsim: public FGInterface {
--- /dev/null
+/*******************************************************************************
+
+ Header: LaRCsimIC.cxx
+ Author: Tony Peden
+ Date started: 10/9/99
+
+ ------------- Copyright (C) 1999 Anthony K. Peden (apeden@earthlink.net) -------------
+
+ This program is free software; you can redistribute it and/or modify it under
+ the terms of the GNU General Public License as published by the Free Software
+ Foundation; either version 2 of the License, or (at your option) any later
+ version.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ details.
+
+ You should have received a copy of the GNU General Public License along with
+ this program; if not, write to the Free Software Foundation, Inc., 59 Temple
+ Place - Suite 330, Boston, MA 02111-1307, USA.
+
+ Further information about the GNU General Public License can also be found on
+ the world wide web at http://www.gnu.org.
+*/
+
+
+#include <simgear/compiler.h>
+
+#include <math.h>
+#include STL_IOSTREAM
+
+#include "ls_cockpit.h"
+#include "ls_generic.h"
+#include "ls_interface.h"
+#include "atmos_62.h"
+#include "ls_constants.h"
+#include "ls_geodesy.h"
+
+#include "LaRCsimIC.hxx"
+
+SG_USING_STD(cout);
+SG_USING_STD(endl);
+
+
+LaRCsimIC::LaRCsimIC(void) {
+ vt=vtg=vw=vc=ve=0;
+ mach=0;
+ alpha=beta=gamma=0;
+ theta=phi=psi=0;
+ altitude=runway_altitude=hdot=alt_agl=0;
+ latitude_gd=latitude_gc=longitude=0;
+ u=v=w=0;
+ vnorth=veast=vdown=0;
+ vnorthwind=veastwind=vdownwind=0;
+ lastSpeedSet=lssetvt;
+ lastAltSet=lssetasl;
+ get_atmosphere();
+ ls_geod_to_geoc( latitude_gd,altitude,&sea_level_radius,&latitude_gc);
+
+}
+
+
+LaRCsimIC::~LaRCsimIC(void) {}
+
+void LaRCsimIC::get_atmosphere(void) {
+ Altitude=altitude; //set LaRCsim variable
+ ls_atmos(Altitude, &density_ratio, &soundspeed, &T, &p);
+ rho=density_ratio*SEA_LEVEL_DENSITY;
+}
+
+void LaRCsimIC::SetVcalibratedKtsIC( SCALAR tt) {
+ vc=tt*KTS_TO_FPS;
+ lastSpeedSet=lssetvc;
+ vt=sqrt(1/density_ratio*vc*vc);
+}
+
+void LaRCsimIC::SetVtrueFpsIC( SCALAR tt) {
+ vt=tt;
+ lastSpeedSet=lssetvt;
+}
+
+void LaRCsimIC::SetMachIC( SCALAR tt) {
+ mach=tt;
+ vt=mach*soundspeed;
+ lastSpeedSet=lssetmach;
+}
+
+void LaRCsimIC::SetVequivalentKtsIC(SCALAR tt) {
+ ve=tt*KTS_TO_FPS;
+ lastSpeedSet=lssetve;
+ vt=sqrt(SEA_LEVEL_DENSITY/rho)*ve;
+}
+
+void LaRCsimIC::SetClimbRateFpmIC( SCALAR tt) {
+ SetClimbRateFpsIC(tt/60.0);
+}
+
+void LaRCsimIC::SetClimbRateFpsIC( SCALAR tt) {
+ vtg=vt+vw;
+ if(vtg > 0.1) {
+ hdot = tt - vdownwind;
+ gamma=asin(hdot/vtg);
+ getTheta();
+ vdown=-1*hdot;
+ }
+}
+
+void LaRCsimIC::SetFlightPathAngleRadIC( SCALAR tt) {
+ gamma=tt;
+ getTheta();
+ vtg=vt+vw;
+ hdot=vtg*sin(tt);
+ vdown=-1*hdot;
+}
+
+void LaRCsimIC::SetPitchAngleRadIC(SCALAR tt) {
+ if(alt_agl < (DEFAULT_AGL_ALT + 0.1) || vt < 10 )
+ theta=DEFAULT_PITCH_ON_GROUND;
+ else
+ theta=tt;
+ getAlpha();
+}
+
+void LaRCsimIC::SetUVWFpsIC(SCALAR vu, SCALAR vv, SCALAR vw) {
+ u=vu; v=vv; w=vw;
+ vt=sqrt(u*u+v*v+w*w);
+ lastSpeedSet=lssetuvw;
+}
+
+
+void LaRCsimIC::SetVNEDAirmassFpsIC(SCALAR north, SCALAR east, SCALAR down ) {
+ vnorthwind=north; veastwind=east; vdownwind=down;
+ vw=sqrt(vnorthwind*vnorthwind + veastwind*veastwind + vdownwind*vdownwind);
+ vtg=vt+vw;
+ SetClimbRateFpsIC(-1*(vdown+vdownwind));
+}
+
+void LaRCsimIC::SetAltitudeFtIC( SCALAR tt) {
+ if(tt > (runway_altitude + DEFAULT_AGL_ALT)) {
+ altitude=tt;
+ } else {
+ altitude=runway_altitude + DEFAULT_AGL_ALT;
+ alt_agl=DEFAULT_AGL_ALT;
+ theta=DEFAULT_PITCH_ON_GROUND;
+ }
+ lastAltSet=lssetasl;
+ get_atmosphere();
+ //lets try to make sure the user gets what they intended
+
+ switch(lastSpeedSet) {
+ case lssetned:
+ calcVtfromNED();
+ break;
+ case lssetuvw:
+ case lssetvt:
+ SetVtrueFpsIC(vt);
+ break;
+ case lssetvc:
+ SetVcalibratedKtsIC(vc*V_TO_KNOTS);
+ break;
+ case lssetve:
+ SetVequivalentKtsIC(ve*V_TO_KNOTS);
+ break;
+ case lssetmach:
+ SetMachIC(mach);
+ break;
+ }
+}
+
+void LaRCsimIC::SetAltitudeAGLFtIC( SCALAR tt) {
+ alt_agl=tt;
+ lastAltSet=lssetagl;
+ altitude=runway_altitude + alt_agl;
+}
+
+void LaRCsimIC::SetRunwayAltitudeFtIC( SCALAR tt) {
+ runway_altitude=tt;
+ if(lastAltSet == lssetagl)
+ altitude = runway_altitude + alt_agl;
+}
+
+void LaRCsimIC::calcVtfromNED(void) {
+ //take the earth's rotation out of veast first
+ //float veastner = veast- OMEGA_EARTH*sea_level_radius*cos( latitude_gd );
+ float veastner=veast;
+ u = (vnorth - vnorthwind)*cos(theta)*cos(psi) +
+ (veastner - veastwind)*cos(theta)*sin(psi) -
+ (vdown - vdownwind)*sin(theta);
+ v = (vnorth - vnorthwind)*(sin(phi)*sin(theta)*cos(psi) - cos(phi)*sin(psi)) +
+ (veastner - veastwind)*(sin(phi)*sin(theta)*sin(psi) + cos(phi)*cos(psi)) +
+ (vdown - vdownwind)*sin(phi)*cos(theta);
+ w = (vnorth - vnorthwind)*(cos(phi)*sin(theta)*cos(psi) + sin(phi)*sin(psi)) +
+ (veastner - veastwind)*(cos(phi)*sin(theta)*sin(psi) - sin(phi)*cos(psi)) +
+ (vdown - vdownwind)*cos(phi)*cos(theta);
+ vt = sqrt(u*u + v*v + w*w);
+}
+
+void LaRCsimIC::calcNEDfromVt(void) {
+ float veastner;
+
+ //get the body components of vt
+ u = GetUBodyFpsIC();
+ v = GetVBodyFpsIC();
+ w = GetWBodyFpsIC();
+
+ //transform them to local axes and add on the wind components
+ vnorth = u*cos(theta)*cos(psi) +
+ v*(sin(phi)*sin(theta)*cos(psi) - cos(phi)*sin(psi)) +
+ w*(cos(phi)*sin(theta)*cos(psi) + sin(phi)*sin(psi)) +
+ vnorthwind;
+
+ //need to account for the earth's rotation here
+ veastner = u*cos(theta)*sin(psi) +
+ v*(sin(phi)*sin(theta)*sin(psi) + cos(phi)*cos(psi)) +
+ w*(cos(phi)*sin(theta)*sin(psi) - sin(phi)*cos(psi)) +
+ veastwind;
+ veast = veastner;
+ //veast = veastner + OMEGA_EARTH*sea_level_radius*cos( latitude_gd );
+
+ vdown = u*sin(theta) +
+ v*sin(phi)*cos(theta) +
+ w*cos(phi)*cos(theta) +
+ vdownwind;
+}
+
+void LaRCsimIC::SetVNEDFpsIC( SCALAR north, SCALAR east, SCALAR down) {
+ vnorth=north;
+ veast=east;
+ vdown=down;
+ SetClimbRateFpsIC(-1*vdown);
+ lastSpeedSet=lssetned;
+ calcVtfromNED();
+}
+
+void LaRCsimIC::SetLatitudeGDRadIC(SCALAR tt) {
+ latitude_gd=tt;
+ ls_geod_to_geoc(latitude_gd,altitude,&sea_level_radius, &latitude_gc);
+}
+
+bool LaRCsimIC::getAlpha(void) {
+ bool result=false;
+ float guess=theta-gamma;
+ xlo=xhi=0;
+ xmin=-89;
+ xmax=89;
+ sfunc=&LaRCsimIC::GammaEqOfAlpha;
+ if(findInterval(0,guess)){
+ if(solve(&alpha,0)){
+ result=true;
+ }
+ }
+ return result;
+}
+
+
+bool LaRCsimIC::getTheta(void) {
+ bool result=false;
+ float guess=alpha+gamma;
+ xlo=xhi=0;
+ xmin=-89;xmax=89;
+ sfunc=&LaRCsimIC::GammaEqOfTheta;
+ if(findInterval(0,guess)){
+ if(solve(&theta,0)){
+ result=true;
+ }
+ }
+ return result;
+}
+
+
+
+SCALAR LaRCsimIC::GammaEqOfTheta( SCALAR theta_arg) {
+ SCALAR a,b,c;
+
+ a=cos(alpha)*cos(beta)*sin(theta_arg);
+ b=sin(beta)*sin(phi);
+ c=sin(alpha)*cos(beta)*cos(phi);
+ return sin(gamma)-a+(b+c)*cos(theta_arg);
+}
+
+SCALAR LaRCsimIC::GammaEqOfAlpha( SCALAR alpha_arg) {
+ float a,b,c;
+
+ a=cos(alpha_arg)*cos(beta)*sin(theta);
+ b=sin(beta)*sin(phi);
+ c=sin(alpha_arg)*cos(beta)*cos(phi);
+ return sin(gamma)-a+(b+c)*cos(theta);
+}
+
+
+
+
+
+
+bool LaRCsimIC::findInterval( SCALAR x,SCALAR guess) {
+ //void find_interval(inter_params &ip,eqfunc f,float y,float constant, int &flag){
+
+ int i=0;
+ bool found=false;
+ float flo,fhi,fguess;
+ float lo,hi,step;
+ step=0.1;
+ fguess=(this->*sfunc)(guess)-x;
+ lo=hi=guess;
+ do {
+ step=2*step;
+ lo-=step;
+ hi+=step;
+ if(lo < xmin) lo=xmin;
+ if(hi > xmax) hi=xmax;
+ i++;
+ flo=(this->*sfunc)(lo)-x;
+ fhi=(this->*sfunc)(hi)-x;
+ if(flo*fhi <=0) { //found interval with root
+ found=true;
+ if(flo*fguess <= 0) { //narrow interval down a bit
+ hi=lo+step; //to pass solver interval that is as
+ //small as possible
+ }
+ else if(fhi*fguess <= 0) {
+ lo=hi-step;
+ }
+ }
+ //cout << "findInterval: i=" << i << " Lo= " << lo << " Hi= " << hi << endl;
+ }
+ while((found == 0) && (i <= 100));
+ xlo=lo;
+ xhi=hi;
+ return found;
+}
+
+
+
+
+bool LaRCsimIC::solve( SCALAR *y,SCALAR x) {
+ float x1,x2,x3,f1,f2,f3,d,d0;
+ float eps=1E-5;
+ float const relax =0.9;
+ int i;
+ bool success=false;
+
+ //initializations
+ d=1;
+
+ x1=xlo;x3=xhi;
+ f1=(this->*sfunc)(x1)-x;
+ f3=(this->*sfunc)(x3)-x;
+ d0=fabs(x3-x1);
+
+ //iterations
+ i=0;
+ while ((fabs(d) > eps) && (i < 100)) {
+ d=(x3-x1)/d0;
+ x2=x1-d*d0*f1/(f3-f1);
+
+ f2=(this->*sfunc)(x2)-x;
+ //cout << "solve x1,x2,x3: " << x1 << "," << x2 << "," << x3 << endl;
+ //cout << " " << f1 << "," << f2 << "," << f3 << endl;
+
+ if(fabs(f2) <= 0.001) {
+ x1=x3=x2;
+ } else if(f1*f2 <= 0.0) {
+ x3=x2;
+ f3=f2;
+ f1=relax*f1;
+ } else if(f2*f3 <= 0) {
+ x1=x2;
+ f1=f2;
+ f3=relax*f3;
+ }
+ //cout << i << endl;
+ i++;
+ }//end while
+ if(i < 100) {
+ success=true;
+ *y=x2;
+ }
+
+ //cout << "Success= " << success << " Vcas: " << vcas*V_TO_KNOTS << " Mach: " << x2 << endl;
+ return success;
+}
--- /dev/null
+/*******************************************************************************
+
+ Header: LaRCsimIC.hxx
+ Author: Tony Peden
+ Date started: 10/9/00
+
+ ------------- Copyright (C) 2000 Anthony K. Peden (apeden@earthlink.net) -------------
+
+ This program is free software; you can redistribute it and/or modify it under
+ the terms of the GNU General Public License as published by the Free Software
+ Foundation; either version 2 of the License, or (at your option) any later
+ version.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ details.
+
+ You should have received a copy of the GNU General Public License along with
+ this program; if not, write to the Free Software Foundation, Inc., 59 Temple
+ Place - Suite 330, Boston, MA 02111-1307, USA.
+
+ Further information about the GNU General Public License can also be found on
+ the world wide web at http://www.gnu.org.
+*/
+#ifndef _LARCSIMIC_HXX
+#define _LARCSIMIC_HXX
+
+/*******************************************************************************
+INCLUDES
+*******************************************************************************/
+
+#include <plib/sg.h>
+
+#include "ls_constants.h"
+#include "ls_types.h"
+
+/*******************************************************************************
+CLASS DECLARATION
+*******************************************************************************/
+
+#define KTS_TO_FPS 1.6889
+#define M_TO_FT 3.2808399
+#define DEFAULT_AGL_ALT 3.758099
+#define DEFAULT_PITCH_ON_GROUND 0.0074002
+
+enum lsspeedset { lssetvt, lssetvc, lssetve, lssetmach, lssetuvw, lssetned };
+enum lsaltset { lssetasl, lssetagl };
+
+
+class LaRCsimIC {
+public:
+
+ LaRCsimIC(void);
+ ~LaRCsimIC(void);
+
+ void SetVcalibratedKtsIC(SCALAR tt);
+ void SetMachIC(SCALAR tt);
+
+ void SetVtrueFpsIC(SCALAR tt);
+
+ void SetVequivalentKtsIC(SCALAR tt);
+
+ void SetUVWFpsIC(SCALAR vu, SCALAR vv, SCALAR vw);
+
+ void SetVNEDFpsIC(SCALAR north, SCALAR east, SCALAR down);
+
+ void SetVNEDAirmassFpsIC(SCALAR north, SCALAR east, SCALAR down );
+
+ void SetAltitudeFtIC(SCALAR tt);
+ void SetAltitudeAGLFtIC(SCALAR tt);
+
+ //"vertical" flight path, recalculate theta
+ inline void SetFlightPathAngleDegIC(SCALAR tt) { SetFlightPathAngleRadIC(tt*SGD_DEGREES_TO_RADIANS); }
+ void SetFlightPathAngleRadIC(SCALAR tt);
+
+ //set speed first
+ void SetClimbRateFpmIC(SCALAR tt);
+ void SetClimbRateFpsIC(SCALAR tt);
+
+ //use currently stored gamma, recalcualte theta
+ inline void SetAlphaDegIC(SCALAR tt) { alpha=tt*SGD_DEGREES_TO_RADIANS; getTheta(); }
+ inline void SetAlphaRadIC(SCALAR tt) { alpha=tt; getTheta(); }
+
+ //use currently stored gamma, recalcualte alpha
+ inline void SetPitchAngleDegIC(SCALAR tt) { SetPitchAngleRadIC(tt*SGD_DEGREES_TO_RADIANS); }
+ void SetPitchAngleRadIC(SCALAR tt);
+
+ inline void SetBetaDegIC(SCALAR tt) { beta=tt*SGD_DEGREES_TO_RADIANS; getTheta();}
+ inline void SetBetaRadIC(SCALAR tt) { beta=tt; getTheta(); }
+
+ inline void SetRollAngleDegIC(SCALAR tt) { phi=tt*SGD_DEGREES_TO_RADIANS; getTheta(); }
+ inline void SetRollAngleRadIC(SCALAR tt) { phi=tt; getTheta(); }
+
+ inline void SetHeadingDegIC(SCALAR tt) { psi=tt*SGD_DEGREES_TO_RADIANS; }
+ inline void SetHeadingRadIC(SCALAR tt) { psi=tt; }
+
+ inline void SetLatitudeGDDegIC(SCALAR tt) { SetLatitudeGDRadIC(tt*SGD_DEGREES_TO_RADIANS); }
+ void SetLatitudeGDRadIC(SCALAR tt);
+
+ inline void SetLongitudeDegIC(SCALAR tt) { longitude=tt*SGD_DEGREES_TO_RADIANS; }
+ inline void SetLongitudeRadIC(SCALAR tt) { longitude=tt; }
+
+ void SetRunwayAltitudeFtIC(SCALAR tt);
+
+ inline SCALAR GetVcalibratedKtsIC(void) { return sqrt(density_ratio*vt*vt)*V_TO_KNOTS; }
+ inline SCALAR GetVequivalentKtsIC(void) { return sqrt(density_ratio)*vt*V_TO_KNOTS; }
+ inline SCALAR GetVtrueKtsIC(void) { return vt*V_TO_KNOTS; }
+ inline SCALAR GetVtrueFpsIC(void) { return vt; }
+ inline SCALAR GetMachIC(void) { return vt/soundspeed; }
+
+ inline SCALAR GetAltitudeFtIC(void) { return altitude; }
+ inline SCALAR GetAltitudeAGLFtIC(void) { return alt_agl; }
+
+ inline SCALAR GetRunwayAltitudeFtIC(void) { return runway_altitude; }
+
+ inline SCALAR GetFlightPathAngleDegIC(void) { return gamma*SGD_RADIANS_TO_DEGREES; }
+ inline SCALAR GetFlightPathAngleRadIC(void) { return gamma; }
+
+ inline SCALAR GetClimbRateFpmIC(void) { return hdot*60; }
+ inline SCALAR GetClimbRateFpsIC(void) { return hdot; }
+
+ inline SCALAR GetAlphaDegIC(void) { return alpha*SGD_RADIANS_TO_DEGREES; }
+ inline SCALAR GetAlphaRadIC(void) { return alpha; }
+
+ inline SCALAR GetPitchAngleDegIC(void) { return theta*SGD_RADIANS_TO_DEGREES; }
+ inline SCALAR GetPitchAngleRadIC(void) { return theta; }
+
+
+ inline SCALAR GetBetaDegIC(void) { return beta*SGD_RADIANS_TO_DEGREES; }
+ inline SCALAR GetBetaRadIC(void) { return beta*SGD_RADIANS_TO_DEGREES; }
+
+ inline SCALAR GetRollAngleDegIC(void) { return phi*SGD_RADIANS_TO_DEGREES; }
+ inline SCALAR GetRollAngleRadIC(void) { return phi; }
+
+ inline SCALAR GetHeadingDegIC(void) { return psi*SGD_RADIANS_TO_DEGREES; }
+ inline SCALAR GetHeadingRadIC(void) { return psi; }
+
+ inline SCALAR GetLatitudeGDDegIC(void) { return latitude_gd*SGD_RADIANS_TO_DEGREES; }
+ inline SCALAR GetLatitudeGDRadIC(void) { return latitude_gd; }
+
+ inline SCALAR GetLongitudeDegIC(void) { return longitude*SGD_RADIANS_TO_DEGREES; }
+ inline SCALAR GetLongitudeRadIC(void) { return longitude; }
+
+ inline SCALAR GetUBodyFpsIC(void) { return vt*cos(alpha)*cos(beta); }
+ inline SCALAR GetVBodyFpsIC(void) { return vt*sin(beta); }
+ inline SCALAR GetWBodyFpsIC(void) { return vt*sin(alpha)*cos(beta); }
+
+ inline SCALAR GetVnorthFpsIC(void) { calcNEDfromVt();return vnorth; }
+ inline SCALAR GetVeastFpsIC(void) { calcNEDfromVt();return veast; }
+ inline SCALAR GetVdownFpsIC(void) { calcNEDfromVt();return vdown; }
+
+ inline SCALAR GetVnorthAirmassFpsIC(void) { return vnorthwind; }
+ inline SCALAR GetVeastAirmassFpsIC(void) { return veastwind; }
+ inline SCALAR GetVdownAirmassFpsIC(void) { return vdownwind; }
+
+ inline SCALAR GetThetaRadIC(void) { return theta; }
+ inline SCALAR GetPhiRadIC(void) { return phi; }
+ inline SCALAR GetPsiRadIC(void) { return psi; }
+
+
+
+private:
+ SCALAR vt,vtg,vw,vc,ve;
+ SCALAR alpha,beta,gamma,theta,phi,psi;
+ SCALAR mach;
+ SCALAR altitude,runway_altitude,hdot,alt_agl,sea_level_radius;
+ SCALAR latitude_gd,latitude_gc,longitude;
+ SCALAR u,v,w;
+ SCALAR vnorth,veast,vdown;
+ SCALAR vnorthwind, veastwind, vdownwind;
+ SCALAR p,T,soundspeed,density_ratio,rho;
+
+ SCALAR xlo, xhi,xmin,xmax;
+
+ typedef SCALAR (LaRCsimIC::*fp)(SCALAR x);
+ fp sfunc;
+
+ lsspeedset lastSpeedSet;
+ lsaltset lastAltSet;
+
+
+ void calcVtfromNED(void);
+ void calcNEDfromVt(void);
+ void calcSpeeds(void);
+
+
+ bool getAlpha(void);
+ bool getTheta(void);
+ SCALAR GammaEqOfTheta(SCALAR tt);
+ SCALAR GammaEqOfAlpha(SCALAR tt);
+ void get_atmosphere(void);
+
+
+ bool findInterval(SCALAR x,SCALAR guess);
+ bool solve(SCALAR *y, SCALAR x);
+};
+
+#endif
libLaRCsim_a_SOURCES = \
LaRCsim.cxx LaRCsim.hxx \
+ LaRCsimIC.cxx LaRCsimIC.hxx \
+ IO360.cxx IO360.hxx \
atmos_62.c atmos_62.h \
default_model_routines.c default_model_routines.h \
ls_accel.c ls_accel.h \
+++ /dev/null
-/*******************************************************************************
-
- Header: LaRCsimIC.cxx
- Author: Tony Peden
- Date started: 10/9/99
-
- ------------- Copyright (C) 1999 Anthony K. Peden (apeden@earthlink.net) -------------
-
- This program is free software; you can redistribute it and/or modify it under
- the terms of the GNU General Public License as published by the Free Software
- Foundation; either version 2 of the License, or (at your option) any later
- version.
-
- This program is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
- FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
- details.
-
- You should have received a copy of the GNU General Public License along with
- this program; if not, write to the Free Software Foundation, Inc., 59 Temple
- Place - Suite 330, Boston, MA 02111-1307, USA.
-
- Further information about the GNU General Public License can also be found on
- the world wide web at http://www.gnu.org.
-*/
-
-
-#include <simgear/compiler.h>
-
-#include <math.h>
-#include STL_IOSTREAM
-
-#include "FDM/LaRCsimIC.hxx"
-#include <FDM/LaRCsim/ls_cockpit.h>
-#include <FDM/LaRCsim/ls_generic.h>
-#include <FDM/LaRCsim/ls_interface.h>
-#include <FDM/LaRCsim/atmos_62.h>
-#include <FDM/LaRCsim/ls_constants.h>
-#include <FDM/LaRCsim/ls_geodesy.h>
-
-SG_USING_STD(cout);
-SG_USING_STD(endl);
-
-
-LaRCsimIC::LaRCsimIC(void) {
- vt=vtg=vw=vc=ve=0;
- mach=0;
- alpha=beta=gamma=0;
- theta=phi=psi=0;
- altitude=runway_altitude=hdot=alt_agl=0;
- latitude_gd=latitude_gc=longitude=0;
- u=v=w=0;
- vnorth=veast=vdown=0;
- vnorthwind=veastwind=vdownwind=0;
- lastSpeedSet=lssetvt;
- lastAltSet=lssetasl;
- get_atmosphere();
- ls_geod_to_geoc( latitude_gd,altitude,&sea_level_radius,&latitude_gc);
-
-}
-
-
-LaRCsimIC::~LaRCsimIC(void) {}
-
-void LaRCsimIC::get_atmosphere(void) {
- Altitude=altitude; //set LaRCsim variable
- ls_atmos(Altitude, &density_ratio, &soundspeed, &T, &p);
- rho=density_ratio*SEA_LEVEL_DENSITY;
-}
-
-void LaRCsimIC::SetVcalibratedKtsIC( SCALAR tt) {
- vc=tt*KTS_TO_FPS;
- lastSpeedSet=lssetvc;
- vt=sqrt(1/density_ratio*vc*vc);
-}
-
-void LaRCsimIC::SetVtrueFpsIC( SCALAR tt) {
- vt=tt;
- lastSpeedSet=lssetvt;
-}
-
-void LaRCsimIC::SetMachIC( SCALAR tt) {
- mach=tt;
- vt=mach*soundspeed;
- lastSpeedSet=lssetmach;
-}
-
-void LaRCsimIC::SetVequivalentKtsIC(SCALAR tt) {
- ve=tt*KTS_TO_FPS;
- lastSpeedSet=lssetve;
- vt=sqrt(SEA_LEVEL_DENSITY/rho)*ve;
-}
-
-void LaRCsimIC::SetClimbRateFpmIC( SCALAR tt) {
- SetClimbRateFpsIC(tt/60.0);
-}
-
-void LaRCsimIC::SetClimbRateFpsIC( SCALAR tt) {
- vtg=vt+vw;
- if(vtg > 0.1) {
- hdot = tt - vdownwind;
- gamma=asin(hdot/vtg);
- getTheta();
- vdown=-1*hdot;
- }
-}
-
-void LaRCsimIC::SetFlightPathAngleRadIC( SCALAR tt) {
- gamma=tt;
- getTheta();
- vtg=vt+vw;
- hdot=vtg*sin(tt);
- vdown=-1*hdot;
-}
-
-void LaRCsimIC::SetPitchAngleRadIC(SCALAR tt) {
- if(alt_agl < (DEFAULT_AGL_ALT + 0.1) || vt < 10 )
- theta=DEFAULT_PITCH_ON_GROUND;
- else
- theta=tt;
- getAlpha();
-}
-
-void LaRCsimIC::SetUVWFpsIC(SCALAR vu, SCALAR vv, SCALAR vw) {
- u=vu; v=vv; w=vw;
- vt=sqrt(u*u+v*v+w*w);
- lastSpeedSet=lssetuvw;
-}
-
-
-void LaRCsimIC::SetVNEDAirmassFpsIC(SCALAR north, SCALAR east, SCALAR down ) {
- vnorthwind=north; veastwind=east; vdownwind=down;
- vw=sqrt(vnorthwind*vnorthwind + veastwind*veastwind + vdownwind*vdownwind);
- vtg=vt+vw;
- SetClimbRateFpsIC(-1*(vdown+vdownwind));
-}
-
-void LaRCsimIC::SetAltitudeFtIC( SCALAR tt) {
- if(tt > (runway_altitude + DEFAULT_AGL_ALT)) {
- altitude=tt;
- } else {
- altitude=runway_altitude + DEFAULT_AGL_ALT;
- alt_agl=DEFAULT_AGL_ALT;
- theta=DEFAULT_PITCH_ON_GROUND;
- }
- lastAltSet=lssetasl;
- get_atmosphere();
- //lets try to make sure the user gets what they intended
-
- switch(lastSpeedSet) {
- case lssetned:
- calcVtfromNED();
- break;
- case lssetuvw:
- case lssetvt:
- SetVtrueFpsIC(vt);
- break;
- case lssetvc:
- SetVcalibratedKtsIC(vc*V_TO_KNOTS);
- break;
- case lssetve:
- SetVequivalentKtsIC(ve*V_TO_KNOTS);
- break;
- case lssetmach:
- SetMachIC(mach);
- break;
- }
-}
-
-void LaRCsimIC::SetAltitudeAGLFtIC( SCALAR tt) {
- alt_agl=tt;
- lastAltSet=lssetagl;
- altitude=runway_altitude + alt_agl;
-}
-
-void LaRCsimIC::SetRunwayAltitudeFtIC( SCALAR tt) {
- runway_altitude=tt;
- if(lastAltSet == lssetagl)
- altitude = runway_altitude + alt_agl;
-}
-
-void LaRCsimIC::calcVtfromNED(void) {
- //take the earth's rotation out of veast first
- //float veastner = veast- OMEGA_EARTH*sea_level_radius*cos( latitude_gd );
- float veastner=veast;
- u = (vnorth - vnorthwind)*cos(theta)*cos(psi) +
- (veastner - veastwind)*cos(theta)*sin(psi) -
- (vdown - vdownwind)*sin(theta);
- v = (vnorth - vnorthwind)*(sin(phi)*sin(theta)*cos(psi) - cos(phi)*sin(psi)) +
- (veastner - veastwind)*(sin(phi)*sin(theta)*sin(psi) + cos(phi)*cos(psi)) +
- (vdown - vdownwind)*sin(phi)*cos(theta);
- w = (vnorth - vnorthwind)*(cos(phi)*sin(theta)*cos(psi) + sin(phi)*sin(psi)) +
- (veastner - veastwind)*(cos(phi)*sin(theta)*sin(psi) - sin(phi)*cos(psi)) +
- (vdown - vdownwind)*cos(phi)*cos(theta);
- vt = sqrt(u*u + v*v + w*w);
-}
-
-void LaRCsimIC::calcNEDfromVt(void) {
- float veastner;
-
- //get the body components of vt
- u = GetUBodyFpsIC();
- v = GetVBodyFpsIC();
- w = GetWBodyFpsIC();
-
- //transform them to local axes and add on the wind components
- vnorth = u*cos(theta)*cos(psi) +
- v*(sin(phi)*sin(theta)*cos(psi) - cos(phi)*sin(psi)) +
- w*(cos(phi)*sin(theta)*cos(psi) + sin(phi)*sin(psi)) +
- vnorthwind;
-
- //need to account for the earth's rotation here
- veastner = u*cos(theta)*sin(psi) +
- v*(sin(phi)*sin(theta)*sin(psi) + cos(phi)*cos(psi)) +
- w*(cos(phi)*sin(theta)*sin(psi) - sin(phi)*cos(psi)) +
- veastwind;
- veast = veastner;
- //veast = veastner + OMEGA_EARTH*sea_level_radius*cos( latitude_gd );
-
- vdown = u*sin(theta) +
- v*sin(phi)*cos(theta) +
- w*cos(phi)*cos(theta) +
- vdownwind;
-}
-
-void LaRCsimIC::SetVNEDFpsIC( SCALAR north, SCALAR east, SCALAR down) {
- vnorth=north;
- veast=east;
- vdown=down;
- SetClimbRateFpsIC(-1*vdown);
- lastSpeedSet=lssetned;
- calcVtfromNED();
-}
-
-void LaRCsimIC::SetLatitudeGDRadIC(SCALAR tt) {
- latitude_gd=tt;
- ls_geod_to_geoc(latitude_gd,altitude,&sea_level_radius, &latitude_gc);
-}
-
-bool LaRCsimIC::getAlpha(void) {
- bool result=false;
- float guess=theta-gamma;
- xlo=xhi=0;
- xmin=-89;
- xmax=89;
- sfunc=&LaRCsimIC::GammaEqOfAlpha;
- if(findInterval(0,guess)){
- if(solve(&alpha,0)){
- result=true;
- }
- }
- return result;
-}
-
-
-bool LaRCsimIC::getTheta(void) {
- bool result=false;
- float guess=alpha+gamma;
- xlo=xhi=0;
- xmin=-89;xmax=89;
- sfunc=&LaRCsimIC::GammaEqOfTheta;
- if(findInterval(0,guess)){
- if(solve(&theta,0)){
- result=true;
- }
- }
- return result;
-}
-
-
-
-SCALAR LaRCsimIC::GammaEqOfTheta( SCALAR theta_arg) {
- SCALAR a,b,c;
-
- a=cos(alpha)*cos(beta)*sin(theta_arg);
- b=sin(beta)*sin(phi);
- c=sin(alpha)*cos(beta)*cos(phi);
- return sin(gamma)-a+(b+c)*cos(theta_arg);
-}
-
-SCALAR LaRCsimIC::GammaEqOfAlpha( SCALAR alpha_arg) {
- float a,b,c;
-
- a=cos(alpha_arg)*cos(beta)*sin(theta);
- b=sin(beta)*sin(phi);
- c=sin(alpha_arg)*cos(beta)*cos(phi);
- return sin(gamma)-a+(b+c)*cos(theta);
-}
-
-
-
-
-
-
-bool LaRCsimIC::findInterval( SCALAR x,SCALAR guess) {
- //void find_interval(inter_params &ip,eqfunc f,float y,float constant, int &flag){
-
- int i=0;
- bool found=false;
- float flo,fhi,fguess;
- float lo,hi,step;
- step=0.1;
- fguess=(this->*sfunc)(guess)-x;
- lo=hi=guess;
- do {
- step=2*step;
- lo-=step;
- hi+=step;
- if(lo < xmin) lo=xmin;
- if(hi > xmax) hi=xmax;
- i++;
- flo=(this->*sfunc)(lo)-x;
- fhi=(this->*sfunc)(hi)-x;
- if(flo*fhi <=0) { //found interval with root
- found=true;
- if(flo*fguess <= 0) { //narrow interval down a bit
- hi=lo+step; //to pass solver interval that is as
- //small as possible
- }
- else if(fhi*fguess <= 0) {
- lo=hi-step;
- }
- }
- //cout << "findInterval: i=" << i << " Lo= " << lo << " Hi= " << hi << endl;
- }
- while((found == 0) && (i <= 100));
- xlo=lo;
- xhi=hi;
- return found;
-}
-
-
-
-
-bool LaRCsimIC::solve( SCALAR *y,SCALAR x) {
- float x1,x2,x3,f1,f2,f3,d,d0;
- float eps=1E-5;
- float const relax =0.9;
- int i;
- bool success=false;
-
- //initializations
- d=1;
-
- x1=xlo;x3=xhi;
- f1=(this->*sfunc)(x1)-x;
- f3=(this->*sfunc)(x3)-x;
- d0=fabs(x3-x1);
-
- //iterations
- i=0;
- while ((fabs(d) > eps) && (i < 100)) {
- d=(x3-x1)/d0;
- x2=x1-d*d0*f1/(f3-f1);
-
- f2=(this->*sfunc)(x2)-x;
- //cout << "solve x1,x2,x3: " << x1 << "," << x2 << "," << x3 << endl;
- //cout << " " << f1 << "," << f2 << "," << f3 << endl;
-
- if(fabs(f2) <= 0.001) {
- x1=x3=x2;
- } else if(f1*f2 <= 0.0) {
- x3=x2;
- f3=f2;
- f1=relax*f1;
- } else if(f2*f3 <= 0) {
- x1=x2;
- f1=f2;
- f3=relax*f3;
- }
- //cout << i << endl;
- i++;
- }//end while
- if(i < 100) {
- success=true;
- *y=x2;
- }
-
- //cout << "Success= " << success << " Vcas: " << vcas*V_TO_KNOTS << " Mach: " << x2 << endl;
- return success;
-}
+++ /dev/null
-/*******************************************************************************
-
- Header: LaRCsimIC.hxx
- Author: Tony Peden
- Date started: 10/9/00
-
- ------------- Copyright (C) 2000 Anthony K. Peden (apeden@earthlink.net) -------------
-
- This program is free software; you can redistribute it and/or modify it under
- the terms of the GNU General Public License as published by the Free Software
- Foundation; either version 2 of the License, or (at your option) any later
- version.
-
- This program is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
- FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
- details.
-
- You should have received a copy of the GNU General Public License along with
- this program; if not, write to the Free Software Foundation, Inc., 59 Temple
- Place - Suite 330, Boston, MA 02111-1307, USA.
-
- Further information about the GNU General Public License can also be found on
- the world wide web at http://www.gnu.org.
-*/
-#ifndef _LARCSIMIC_HXX
-#define _LARCSIMIC_HXX
-
-/*******************************************************************************
-INCLUDES
-*******************************************************************************/
-
-#include <plib/sg.h>
-
-#include <FDM/LaRCsim/ls_constants.h>
-#include <FDM/LaRCsim/ls_types.h>
-
-/*******************************************************************************
-CLASS DECLARATION
-*******************************************************************************/
-
-#define KTS_TO_FPS 1.6889
-#define M_TO_FT 3.2808399
-#define DEFAULT_AGL_ALT 3.758099
-#define DEFAULT_PITCH_ON_GROUND 0.0074002
-
-enum lsspeedset { lssetvt, lssetvc, lssetve, lssetmach, lssetuvw, lssetned };
-enum lsaltset { lssetasl, lssetagl };
-
-
-class LaRCsimIC {
-public:
-
- LaRCsimIC(void);
- ~LaRCsimIC(void);
-
- void SetVcalibratedKtsIC(SCALAR tt);
- void SetMachIC(SCALAR tt);
-
- void SetVtrueFpsIC(SCALAR tt);
-
- void SetVequivalentKtsIC(SCALAR tt);
-
- void SetUVWFpsIC(SCALAR vu, SCALAR vv, SCALAR vw);
-
- void SetVNEDFpsIC(SCALAR north, SCALAR east, SCALAR down);
-
- void SetVNEDAirmassFpsIC(SCALAR north, SCALAR east, SCALAR down );
-
- void SetAltitudeFtIC(SCALAR tt);
- void SetAltitudeAGLFtIC(SCALAR tt);
-
- //"vertical" flight path, recalculate theta
- inline void SetFlightPathAngleDegIC(SCALAR tt) { SetFlightPathAngleRadIC(tt*SGD_DEGREES_TO_RADIANS); }
- void SetFlightPathAngleRadIC(SCALAR tt);
-
- //set speed first
- void SetClimbRateFpmIC(SCALAR tt);
- void SetClimbRateFpsIC(SCALAR tt);
-
- //use currently stored gamma, recalcualte theta
- inline void SetAlphaDegIC(SCALAR tt) { alpha=tt*SGD_DEGREES_TO_RADIANS; getTheta(); }
- inline void SetAlphaRadIC(SCALAR tt) { alpha=tt; getTheta(); }
-
- //use currently stored gamma, recalcualte alpha
- inline void SetPitchAngleDegIC(SCALAR tt) { SetPitchAngleRadIC(tt*SGD_DEGREES_TO_RADIANS); }
- void SetPitchAngleRadIC(SCALAR tt);
-
- inline void SetBetaDegIC(SCALAR tt) { beta=tt*SGD_DEGREES_TO_RADIANS; getTheta();}
- inline void SetBetaRadIC(SCALAR tt) { beta=tt; getTheta(); }
-
- inline void SetRollAngleDegIC(SCALAR tt) { phi=tt*SGD_DEGREES_TO_RADIANS; getTheta(); }
- inline void SetRollAngleRadIC(SCALAR tt) { phi=tt; getTheta(); }
-
- inline void SetHeadingDegIC(SCALAR tt) { psi=tt*SGD_DEGREES_TO_RADIANS; }
- inline void SetHeadingRadIC(SCALAR tt) { psi=tt; }
-
- inline void SetLatitudeGDDegIC(SCALAR tt) { SetLatitudeGDRadIC(tt*SGD_DEGREES_TO_RADIANS); }
- void SetLatitudeGDRadIC(SCALAR tt);
-
- inline void SetLongitudeDegIC(SCALAR tt) { longitude=tt*SGD_DEGREES_TO_RADIANS; }
- inline void SetLongitudeRadIC(SCALAR tt) { longitude=tt; }
-
- void SetRunwayAltitudeFtIC(SCALAR tt);
-
- inline SCALAR GetVcalibratedKtsIC(void) { return sqrt(density_ratio*vt*vt)*V_TO_KNOTS; }
- inline SCALAR GetVequivalentKtsIC(void) { return sqrt(density_ratio)*vt*V_TO_KNOTS; }
- inline SCALAR GetVtrueKtsIC(void) { return vt*V_TO_KNOTS; }
- inline SCALAR GetVtrueFpsIC(void) { return vt; }
- inline SCALAR GetMachIC(void) { return vt/soundspeed; }
-
- inline SCALAR GetAltitudeFtIC(void) { return altitude; }
- inline SCALAR GetAltitudeAGLFtIC(void) { return alt_agl; }
-
- inline SCALAR GetRunwayAltitudeFtIC(void) { return runway_altitude; }
-
- inline SCALAR GetFlightPathAngleDegIC(void) { return gamma*SGD_RADIANS_TO_DEGREES; }
- inline SCALAR GetFlightPathAngleRadIC(void) { return gamma; }
-
- inline SCALAR GetClimbRateFpmIC(void) { return hdot*60; }
- inline SCALAR GetClimbRateFpsIC(void) { return hdot; }
-
- inline SCALAR GetAlphaDegIC(void) { return alpha*SGD_RADIANS_TO_DEGREES; }
- inline SCALAR GetAlphaRadIC(void) { return alpha; }
-
- inline SCALAR GetPitchAngleDegIC(void) { return theta*SGD_RADIANS_TO_DEGREES; }
- inline SCALAR GetPitchAngleRadIC(void) { return theta; }
-
-
- inline SCALAR GetBetaDegIC(void) { return beta*SGD_RADIANS_TO_DEGREES; }
- inline SCALAR GetBetaRadIC(void) { return beta*SGD_RADIANS_TO_DEGREES; }
-
- inline SCALAR GetRollAngleDegIC(void) { return phi*SGD_RADIANS_TO_DEGREES; }
- inline SCALAR GetRollAngleRadIC(void) { return phi; }
-
- inline SCALAR GetHeadingDegIC(void) { return psi*SGD_RADIANS_TO_DEGREES; }
- inline SCALAR GetHeadingRadIC(void) { return psi; }
-
- inline SCALAR GetLatitudeGDDegIC(void) { return latitude_gd*SGD_RADIANS_TO_DEGREES; }
- inline SCALAR GetLatitudeGDRadIC(void) { return latitude_gd; }
-
- inline SCALAR GetLongitudeDegIC(void) { return longitude*SGD_RADIANS_TO_DEGREES; }
- inline SCALAR GetLongitudeRadIC(void) { return longitude; }
-
- inline SCALAR GetUBodyFpsIC(void) { return vt*cos(alpha)*cos(beta); }
- inline SCALAR GetVBodyFpsIC(void) { return vt*sin(beta); }
- inline SCALAR GetWBodyFpsIC(void) { return vt*sin(alpha)*cos(beta); }
-
- inline SCALAR GetVnorthFpsIC(void) { calcNEDfromVt();return vnorth; }
- inline SCALAR GetVeastFpsIC(void) { calcNEDfromVt();return veast; }
- inline SCALAR GetVdownFpsIC(void) { calcNEDfromVt();return vdown; }
-
- inline SCALAR GetVnorthAirmassFpsIC(void) { return vnorthwind; }
- inline SCALAR GetVeastAirmassFpsIC(void) { return veastwind; }
- inline SCALAR GetVdownAirmassFpsIC(void) { return vdownwind; }
-
- inline SCALAR GetThetaRadIC(void) { return theta; }
- inline SCALAR GetPhiRadIC(void) { return phi; }
- inline SCALAR GetPsiRadIC(void) { return psi; }
-
-
-
-private:
- SCALAR vt,vtg,vw,vc,ve;
- SCALAR alpha,beta,gamma,theta,phi,psi;
- SCALAR mach;
- SCALAR altitude,runway_altitude,hdot,alt_agl,sea_level_radius;
- SCALAR latitude_gd,latitude_gc,longitude;
- SCALAR u,v,w;
- SCALAR vnorth,veast,vdown;
- SCALAR vnorthwind, veastwind, vdownwind;
- SCALAR p,T,soundspeed,density_ratio,rho;
-
- SCALAR xlo, xhi,xmin,xmax;
-
- typedef SCALAR (LaRCsimIC::*fp)(SCALAR x);
- fp sfunc;
-
- lsspeedset lastSpeedSet;
- lsaltset lastAltSet;
-
-
- void calcVtfromNED(void);
- void calcNEDfromVt(void);
- void calcSpeeds(void);
-
-
- bool getAlpha(void);
- bool getTheta(void);
- SCALAR GammaEqOfTheta(SCALAR tt);
- SCALAR GammaEqOfAlpha(SCALAR tt);
- void get_atmosphere(void);
-
-
- bool findInterval(SCALAR x,SCALAR guess);
- bool solve(SCALAR *y, SCALAR x);
-};
-
-#endif
ADA.cxx ADA.hxx \
Balloon.cxx Balloon.h \
flight.cxx flight.hxx \
- IO360.cxx IO360.hxx \
- LaRCsimIC.cxx LaRCsimIC.hxx \
MagicCarpet.cxx MagicCarpet.hxx \
UFO.cxx UFO.hxx \
NullFDM.cxx NullFDM.hxx
projects / flightgear.git / commitdiff