1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 --------- Copyright (C) 1999 Anthony K. Peden (apeden@earthlink.net) ---------
9 This program is free software; you can redistribute it and/or modify it under
10 the terms of the GNU General Public License as published by the Free Software
11 Foundation; either version 2 of the License, or (at your option) any later
14 This program is distributed in the hope that it will be useful, but WITHOUT
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
19 You should have received a copy of the GNU General Public License along with
20 this program; if not, write to the Free Software Foundation, Inc., 59 Temple
21 Place - Suite 330, Boston, MA 02111-1307, USA.
23 Further information about the GNU General Public License can also be found on
24 the world wide web at http://www.gnu.org.
28 --------------------------------------------------------------------------------
31 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
38 #include "FGFDMExec.h"
39 #include "FGAtmosphere.h"
40 #include "FGInitialCondition.h"
41 #include "FGTrimAxis.h"
42 #include "FGAircraft.h"
43 #include "FGPropulsion.h"
45 static const char *IdSrc = "$Id$";
46 static const char *IdHdr = ID_TRIMAXIS;
48 extern short debug_lvl;
50 /*****************************************************************************/
52 FGTrimAxis::FGTrimAxis(FGFDMExec* fdex, FGInitialCondition* ic, State st,
62 its_to_stable_value=0;
64 total_stability_iterations=0;
75 control_min=-30*DEGTORAD;
76 control_max=30*DEGTORAD;
77 control_convert=RADTODEG;
80 control_min=fdmex->GetAircraft()->GetAlphaCLMin();
81 control_max=fdmex->GetAircraft()->GetAlphaCLMax();
82 if(control_max <= control_min) {
83 control_max=20*DEGTORAD;
84 control_min=-5*DEGTORAD;
86 control_value= (control_min+control_max)/2;
87 control_convert=RADTODEG;
88 solver_eps=tolerance/100;
98 state_convert=RADTODEG;
99 solver_eps=tolerance/100;
104 control_value=fdmex->GetPosition()->GetDistanceAGL();
105 solver_eps=tolerance/100;
108 control_min=fdmex->GetRotation()->Gettht() - 5*DEGTORAD;
109 control_max=fdmex->GetRotation()->Gettht() + 5*DEGTORAD;
110 state_convert=RADTODEG;
113 control_min=fdmex->GetRotation()->Getphi() - 30*DEGTORAD;
114 control_max=fdmex->GetRotation()->Getphi() + 30*DEGTORAD;
115 state_convert=RADTODEG;
116 control_convert=RADTODEG;
119 solver_eps=tolerance/100;
120 control_min=-80*DEGTORAD;
121 control_max=80*DEGTORAD;
122 control_convert=RADTODEG;
125 control_min=fdmex->GetRotation()->Getpsi() - 30*DEGTORAD;
126 control_max=fdmex->GetRotation()->Getpsi() + 30*DEGTORAD;
127 state_convert=RADTODEG;
132 case tUdot: tolerance = DEFAULT_TOLERANCE; break;
133 case tVdot: tolerance = DEFAULT_TOLERANCE; break;
134 case tWdot: tolerance = DEFAULT_TOLERANCE; break;
135 case tQdot: tolerance = DEFAULT_TOLERANCE / 10; break;
136 case tPdot: tolerance = DEFAULT_TOLERANCE / 10; break;
137 case tRdot: tolerance = DEFAULT_TOLERANCE / 10; break;
138 case tHmgt: tolerance = 0.01; break;
141 if (debug_lvl & 2) cout << "Instantiated: FGTrimAxis" << endl;
144 /*****************************************************************************/
146 FGTrimAxis::~FGTrimAxis()
148 if (debug_lvl & 2) cout << "Destroyed: FGTrimAxis" << endl;
151 /*****************************************************************************/
153 void FGTrimAxis::getState(void) {
155 case tUdot: state_value=fdmex->GetTranslation()->GetUVWdot(1); break;
156 case tVdot: state_value=fdmex->GetTranslation()->GetUVWdot(2); break;
157 case tWdot: state_value=fdmex->GetTranslation()->GetUVWdot(3); break;
158 case tQdot: state_value=fdmex->GetRotation()->GetPQRdot(2);break;
159 case tPdot: state_value=fdmex->GetRotation()->GetPQRdot(1); break;
160 case tRdot: state_value=fdmex->GetRotation()->GetPQRdot(3); break;
161 case tHmgt: state_value=computeHmgt(); break;
165 /*****************************************************************************/
167 //States are not settable
169 void FGTrimAxis::getControl(void) {
171 case tThrottle: control_value=fdmex->GetFCS()->GetThrottleCmd(0); break;
172 case tBeta: control_value=fdmex->GetTranslation()->Getalpha(); break;
173 case tAlpha: control_value=fdmex->GetTranslation()->Getbeta(); break;
174 case tPitchTrim: control_value=fdmex->GetFCS() -> GetPitchTrimCmd(); break;
175 case tElevator: control_value=fdmex->GetFCS() -> GetDeCmd(); break;
177 case tAileron: control_value=fdmex->GetFCS() -> GetDaCmd(); break;
179 case tRudder: control_value=fdmex->GetFCS() -> GetDrCmd(); break;
180 case tAltAGL: control_value=fdmex->GetPosition()->GetDistanceAGL();break;
181 case tTheta: control_value=fdmex->GetRotation()->Gettht(); break;
182 case tPhi: control_value=fdmex->GetRotation()->Getphi(); break;
183 case tGamma: control_value=fdmex->GetPosition()->GetGamma();break;
184 case tHeading: control_value=fdmex->GetRotation()->Getpsi(); break;
188 /*****************************************************************************/
190 float FGTrimAxis::computeHmgt(void) {
193 diff = fdmex->GetRotation()->Getpsi() -
194 fdmex->GetPosition()->GetGroundTrack();
197 return (diff + 2*M_PI);
198 } else if( diff > M_PI ) {
199 return (diff - 2*M_PI);
206 /*****************************************************************************/
209 void FGTrimAxis::setControl(void) {
211 case tThrottle: setThrottlesPct(); break;
212 case tBeta: fgic->SetBetaRadIC(control_value); break;
213 case tAlpha: fgic->SetAlphaRadIC(control_value); break;
214 case tPitchTrim: fdmex->GetFCS()->SetPitchTrimCmd(control_value); break;
215 case tElevator: fdmex->GetFCS()->SetDeCmd(control_value); break;
217 case tAileron: fdmex->GetFCS()->SetDaCmd(control_value); break;
219 case tRudder: fdmex->GetFCS()->SetDrCmd(control_value); break;
220 case tAltAGL: fgic->SetAltitudeAGLFtIC(control_value); break;
221 case tTheta: fgic->SetPitchAngleRadIC(control_value); break;
222 case tPhi: fgic->SetRollAngleRadIC(control_value); break;
223 case tGamma: fgic->SetFlightPathAngleRadIC(control_value); break;
224 case tHeading: fgic->SetTrueHeadingRadIC(control_value); break;
232 /*****************************************************************************/
234 // the aircraft center of rotation is no longer the cg once the gear
235 // contact the ground so the altitude needs to be changed when pitch
236 // and roll angle are adjusted. Instead of attempting to calculate the
237 // new center of rotation, pick a gear unit as a reference and use its
238 // location vector to calculate the new height change. i.e. new altitude =
239 // earth z component of that vector (which is in body axes )
240 void FGTrimAxis::SetThetaOnGround(float ff) {
243 // favor an off-center unit so that the same one can be used for both
244 // pitch and roll. An on-center unit is used (for pitch)if that's all
245 // that's in contact with the ground.
246 i=0; ref=-1; center=-1;
247 while( (ref < 0) && (i < fdmex->GetAircraft()->GetNumGearUnits()) ) {
248 if(fdmex->GetAircraft()->GetGearUnit(i)->GetWOW()) {
249 if(fabs(fdmex->GetAircraft()->GetGearUnit(i)->GetBodyLocation(2)) > 0.01)
256 if((ref < 0) && (center >= 0)) {
259 cout << "SetThetaOnGround ref gear: " << ref << endl;
261 float sp=fdmex->GetRotation()->GetSinphi();
262 float cp=fdmex->GetRotation()->GetCosphi();
263 float lx=fdmex->GetAircraft()->GetGearUnit(ref)->GetBodyLocation(1);
264 float ly=fdmex->GetAircraft()->GetGearUnit(ref)->GetBodyLocation(2);
265 float lz=fdmex->GetAircraft()->GetGearUnit(ref)->GetBodyLocation(3);
266 float hagl = -1*lx*sin(ff) +
270 fgic->SetAltitudeAGLFtIC(hagl);
271 cout << "SetThetaOnGround new alt: " << hagl << endl;
273 fgic->SetPitchAngleRadIC(ff);
274 cout << "SetThetaOnGround new theta: " << ff << endl;
277 /*****************************************************************************/
279 bool FGTrimAxis::initTheta(void) {
280 int i,N,iAft, iForward;
281 float zAft,zForward,zDiff,theta;
285 saveAlt=fgic->GetAltitudeAGLFtIC();
286 fgic->SetAltitudeAGLFtIC(100);
289 N=fdmex->GetAircraft()->GetNumGearUnits();
291 //find the first wheel unit forward of the cg
292 //the list is short so a simple linear search is fine
293 for( i=0; i<N; i++ ) {
294 if(fdmex->GetAircraft()->GetGearUnit(i)->GetBodyLocation(1) > 0 ) {
299 //now find the first wheel unit aft of the cg
300 for( i=0; i<N; i++ ) {
301 if(fdmex->GetAircraft()->GetGearUnit(i)->GetBodyLocation(1) < 0 ) {
307 // now adjust theta till the wheels are the same distance from the ground
308 zAft=fdmex->GetAircraft()->GetGearUnit(1)->GetLocalGear(3);
309 zForward=fdmex->GetAircraft()->GetGearUnit(0)->GetLocalGear(3);
310 zDiff = zForward - zAft;
312 theta=fgic->GetPitchAngleDegIC();
313 while(!level && (i < 100)) {
315 fgic->SetPitchAngleDegIC(theta);
317 zAft=fdmex->GetAircraft()->GetGearUnit(1)->GetLocalGear(3);
318 zForward=fdmex->GetAircraft()->GetGearUnit(0)->GetLocalGear(3);
319 zDiff = zForward - zAft;
320 //cout << endl << theta << " " << zDiff << endl;
321 //cout << "0: " << fdmex->GetAircraft()->GetGearUnit(0)->GetLocalGear() << endl;
322 //cout << "1: " << fdmex->GetAircraft()->GetGearUnit(1)->GetLocalGear() << endl;
324 if(fabs(zDiff ) < 0.1)
329 cout << " Initial Theta: " << fdmex->GetRotation()->Gettht()*RADTODEG << endl;
330 control_min=(theta+5)*DEGTORAD;
331 control_max=(theta-5)*DEGTORAD;
332 fgic->SetAltitudeAGLFtIC(saveAlt);
339 /*****************************************************************************/
341 void FGTrimAxis::SetPhiOnGround(float ff) {
345 //must have an off-center unit here
346 while( (ref < 0) && (i < fdmex->GetAircraft()->GetNumGearUnits()) ) {
347 if( (fdmex->GetAircraft()->GetGearUnit(i)->GetWOW()) &&
348 (fabs(fdmex->GetAircraft()->GetGearUnit(i)->GetBodyLocation(2)) > 0.01))
353 float st=fdmex->GetRotation()->GetSintht();
354 float ct=fdmex->GetRotation()->GetCostht();
355 float lx=fdmex->GetAircraft()->GetGearUnit(ref)->GetBodyLocation(1);
356 float ly=fdmex->GetAircraft()->GetGearUnit(ref)->GetBodyLocation(2);
357 float lz=fdmex->GetAircraft()->GetGearUnit(ref)->GetBodyLocation(3);
358 float hagl = -1*lx*st +
362 fgic->SetAltitudeAGLFtIC(hagl);
364 fgic->SetRollAngleRadIC(ff);
368 /*****************************************************************************/
370 void FGTrimAxis::Run(void) {
372 float last_state_value;
375 //cout << "FGTrimAxis::Run: " << control_value << endl;
380 last_state_value=state_value;
384 if((fabs(last_state_value - state_value) < tolerance) || (i >= 100) )
389 its_to_stable_value=i;
390 total_stability_iterations+=its_to_stable_value;
394 /*****************************************************************************/
396 void FGTrimAxis::setThrottlesPct(void) {
398 for(unsigned i=0;i<fdmex->GetPropulsion()->GetNumEngines();i++) {
399 tMin=fdmex->GetPropulsion()->GetEngine(i)->GetThrottleMin();
400 tMax=fdmex->GetPropulsion()->GetEngine(i)->GetThrottleMax();
401 //cout << "setThrottlespct: " << i << ", " << control_min << ", " << control_max << ", " << control_value;
402 fdmex -> GetFCS() -> SetThrottleCmd(i,tMin+control_value*(tMax-tMin));
406 /*****************************************************************************/
408 void FGTrimAxis::AxisReport(void) {
411 sprintf(out," %20s: %6.2f %5s: %9.2e Tolerance: %3.0e\n",
412 GetControlName().c_str(), GetControl()*control_convert,
413 GetStateName().c_str(), GetState(), GetTolerance());
418 /*****************************************************************************/
420 float FGTrimAxis::GetAvgStability( void ) {
421 if(total_iterations > 0) {
422 return float(total_stability_iterations)/float(total_iterations);
427 /*****************************************************************************/
429 void FGTrimAxis::Debug(void)