From: curt Date: Tue, 20 May 2003 11:29:06 +0000 (+0000) Subject: Oops missed a couple things when I moved LaRCsim.cxx into src/FDM/LaRCsim/ X-Git-Url: https://git.mxchange.org/?a=commitdiff_plain;h=676e4c8846bee1dd92bf4b2b34797c4194995dce;p=flightgear.git Oops missed a couple things when I moved LaRCsim.cxx into src/FDM/LaRCsim/ This was masked because I didn't wipe src/FDM/libFlight.a and recreate it. --- diff --git a/src/FDM/IO360.cxx b/src/FDM/IO360.cxx deleted file mode 100644 index 2c2d0db23..000000000 --- a/src/FDM/IO360.cxx +++ /dev/null @@ -1,696 +0,0 @@ -// 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 - -#include - -#include STL_FSTREAM -#include STL_IOSTREAM - -#include "IO360.hxx" -#include "LaRCsim/ls_constants.h" - -#include
- -//************************************************************************************* -// 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 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 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'; -} diff --git a/src/FDM/IO360.hxx b/src/FDM/IO360.hxx deleted file mode 100644 index 1f0df8d7c..000000000 --- a/src/FDM/IO360.hxx +++ /dev/null @@ -1,245 +0,0 @@ -// 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 - -#include - -#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_ diff --git a/src/FDM/LaRCsim/IO360.cxx b/src/FDM/LaRCsim/IO360.cxx new file mode 100644 index 000000000..16483848f --- /dev/null +++ b/src/FDM/LaRCsim/IO360.cxx @@ -0,0 +1,697 @@ +// 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 + +#include + +#include STL_FSTREAM +#include STL_IOSTREAM + +#include
+ +#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 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 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'; +} diff --git a/src/FDM/LaRCsim/IO360.hxx b/src/FDM/LaRCsim/IO360.hxx new file mode 100644 index 000000000..1f0df8d7c --- /dev/null +++ b/src/FDM/LaRCsim/IO360.hxx @@ -0,0 +1,245 @@ +// 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 + +#include + +#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_ diff --git a/src/FDM/LaRCsim/LaRCsim.cxx b/src/FDM/LaRCsim/LaRCsim.cxx index 4235db3a0..32da763b5 100644 --- a/src/FDM/LaRCsim/LaRCsim.cxx +++ b/src/FDM/LaRCsim/LaRCsim.cxx @@ -31,8 +31,6 @@ #include #include #include -#include -#include #include #include
#include @@ -42,6 +40,8 @@ #include "ls_interface.h" #include "ls_constants.h" +#include "IO360.hxx" +#include "LaRCsimIC.hxx" #include "LaRCsim.hxx" diff --git a/src/FDM/LaRCsim/LaRCsim.hxx b/src/FDM/LaRCsim/LaRCsim.hxx index 1e0e6d948..f729e68fb 100644 --- a/src/FDM/LaRCsim/LaRCsim.hxx +++ b/src/FDM/LaRCsim/LaRCsim.hxx @@ -27,9 +27,10 @@ #define _LARCSIM_HXX -#include #include -#include + +#include "IO360.hxx" +#include "LaRCsimIC.hxx" class FGLaRCsim: public FGInterface { diff --git a/src/FDM/LaRCsim/LaRCsimIC.cxx b/src/FDM/LaRCsim/LaRCsimIC.cxx new file mode 100644 index 000000000..562d0c5fa --- /dev/null +++ b/src/FDM/LaRCsim/LaRCsimIC.cxx @@ -0,0 +1,382 @@ +/******************************************************************************* + + 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 + +#include +#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; +} diff --git a/src/FDM/LaRCsim/LaRCsimIC.hxx b/src/FDM/LaRCsim/LaRCsimIC.hxx new file mode 100644 index 000000000..e927fd90f --- /dev/null +++ b/src/FDM/LaRCsim/LaRCsimIC.hxx @@ -0,0 +1,199 @@ +/******************************************************************************* + + 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 + +#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 diff --git a/src/FDM/LaRCsim/Makefile.am b/src/FDM/LaRCsim/Makefile.am index b8f7b68de..f98b7f31e 100644 --- a/src/FDM/LaRCsim/Makefile.am +++ b/src/FDM/LaRCsim/Makefile.am @@ -17,6 +17,8 @@ noinst_LIBRARIES = libLaRCsim.a 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 \ diff --git a/src/FDM/LaRCsimIC.cxx b/src/FDM/LaRCsimIC.cxx deleted file mode 100644 index 2bc6caa05..000000000 --- a/src/FDM/LaRCsimIC.cxx +++ /dev/null @@ -1,381 +0,0 @@ -/******************************************************************************* - - 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 - -#include -#include STL_IOSTREAM - -#include "FDM/LaRCsimIC.hxx" -#include -#include -#include -#include -#include -#include - -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; -} diff --git a/src/FDM/LaRCsimIC.hxx b/src/FDM/LaRCsimIC.hxx deleted file mode 100644 index 7f2bb0d85..000000000 --- a/src/FDM/LaRCsimIC.hxx +++ /dev/null @@ -1,199 +0,0 @@ -/******************************************************************************* - - 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 - -#include -#include - -/******************************************************************************* -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 diff --git a/src/FDM/Makefile.am b/src/FDM/Makefile.am index 3458725c4..6d87e9372 100644 --- a/src/FDM/Makefile.am +++ b/src/FDM/Makefile.am @@ -6,8 +6,6 @@ libFlight_a_SOURCES = \ 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