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 /*****************************************************************************/
50 FGTrimAxis::FGTrimAxis(FGFDMExec* fdex, FGInitialCondition* ic, State st,
59 its_to_stable_value=0;
61 total_stability_iterations=0;
67 case tUdot: tolerance = DEFAULT_TOLERANCE; break;
68 case tVdot: tolerance = DEFAULT_TOLERANCE; break;
69 case tWdot: tolerance = DEFAULT_TOLERANCE; break;
70 case tQdot: tolerance = DEFAULT_TOLERANCE / 10; break;
71 case tPdot: tolerance = DEFAULT_TOLERANCE / 10; break;
72 case tRdot: tolerance = DEFAULT_TOLERANCE / 10; break;
73 case tHmgt: tolerance = 0.01; break;
74 case tNlf: state_target=1.0; tolerance = 1E-5; break;
86 control_min=-30*degtorad;
87 control_max=30*degtorad;
88 control_convert=radtodeg;
91 control_min=fdmex->GetAircraft()->GetAlphaCLMin();
92 control_max=fdmex->GetAircraft()->GetAlphaCLMax();
93 if(control_max <= control_min) {
94 control_max=20*degtorad;
95 control_min=-5*degtorad;
97 control_value= (control_min+control_max)/2;
98 control_convert=radtodeg;
99 solver_eps=tolerance/100;
109 state_convert=radtodeg;
110 solver_eps=tolerance/100;
115 control_value=fdmex->GetPosition()->GetDistanceAGL();
116 solver_eps=tolerance/100;
119 control_min=fdmex->GetRotation()->Gettht() - 5*degtorad;
120 control_max=fdmex->GetRotation()->Gettht() + 5*degtorad;
121 state_convert=radtodeg;
124 control_min=fdmex->GetRotation()->Getphi() - 30*degtorad;
125 control_max=fdmex->GetRotation()->Getphi() + 30*degtorad;
126 state_convert=radtodeg;
127 control_convert=radtodeg;
130 solver_eps=tolerance/100;
131 control_min=-80*degtorad;
132 control_max=80*degtorad;
133 control_convert=radtodeg;
136 control_min=fdmex->GetRotation()->Getpsi() - 30*degtorad;
137 control_max=fdmex->GetRotation()->Getpsi() + 30*degtorad;
138 state_convert=radtodeg;
146 /*****************************************************************************/
148 FGTrimAxis::~FGTrimAxis(void)
153 /*****************************************************************************/
155 void FGTrimAxis::getState(void) {
157 case tUdot: state_value=fdmex->GetTranslation()->GetUVWdot(1)-state_target; break;
158 case tVdot: state_value=fdmex->GetTranslation()->GetUVWdot(2)-state_target; break;
159 case tWdot: state_value=fdmex->GetTranslation()->GetUVWdot(3)-state_target; break;
160 case tQdot: state_value=fdmex->GetRotation()->GetPQRdot(2)-state_target;break;
161 case tPdot: state_value=fdmex->GetRotation()->GetPQRdot(1)-state_target; break;
162 case tRdot: state_value=fdmex->GetRotation()->GetPQRdot(3)-state_target; break;
163 case tHmgt: state_value=computeHmgt()-state_target; break;
164 case tNlf: state_value=fdmex->GetAircraft()->GetNlf()-state_target; break;
169 /*****************************************************************************/
171 //States are not settable
173 void FGTrimAxis::getControl(void) {
175 case tThrottle: control_value=fdmex->GetFCS()->GetThrottleCmd(0); break;
176 case tBeta: control_value=fdmex->GetTranslation()->Getalpha(); break;
177 case tAlpha: control_value=fdmex->GetTranslation()->Getbeta(); break;
178 case tPitchTrim: control_value=fdmex->GetFCS() -> GetPitchTrimCmd(); break;
179 case tElevator: control_value=fdmex->GetFCS() -> GetDeCmd(); break;
181 case tAileron: control_value=fdmex->GetFCS() -> GetDaCmd(); break;
183 case tRudder: control_value=fdmex->GetFCS() -> GetDrCmd(); break;
184 case tAltAGL: control_value=fdmex->GetPosition()->GetDistanceAGL();break;
185 case tTheta: control_value=fdmex->GetRotation()->Gettht(); break;
186 case tPhi: control_value=fdmex->GetRotation()->Getphi(); break;
187 case tGamma: control_value=fdmex->GetPosition()->GetGamma();break;
188 case tHeading: control_value=fdmex->GetRotation()->Getpsi(); break;
192 /*****************************************************************************/
194 double FGTrimAxis::computeHmgt(void) {
197 diff = fdmex->GetRotation()->Getpsi() -
198 fdmex->GetPosition()->GetGroundTrack();
201 return (diff + 2*M_PI);
202 } else if( diff > M_PI ) {
203 return (diff - 2*M_PI);
210 /*****************************************************************************/
213 void FGTrimAxis::setControl(void) {
215 case tThrottle: setThrottlesPct(); break;
216 case tBeta: fgic->SetBetaRadIC(control_value); break;
217 case tAlpha: fgic->SetAlphaRadIC(control_value); break;
218 case tPitchTrim: fdmex->GetFCS()->SetPitchTrimCmd(control_value); break;
219 case tElevator: fdmex->GetFCS()->SetDeCmd(control_value); break;
221 case tAileron: fdmex->GetFCS()->SetDaCmd(control_value); break;
223 case tRudder: fdmex->GetFCS()->SetDrCmd(control_value); break;
224 case tAltAGL: fgic->SetAltitudeAGLFtIC(control_value); break;
225 case tTheta: fgic->SetPitchAngleRadIC(control_value); break;
226 case tPhi: fgic->SetRollAngleRadIC(control_value); break;
227 case tGamma: fgic->SetFlightPathAngleRadIC(control_value); break;
228 case tHeading: fgic->SetTrueHeadingRadIC(control_value); break;
236 /*****************************************************************************/
238 // the aircraft center of rotation is no longer the cg once the gear
239 // contact the ground so the altitude needs to be changed when pitch
240 // and roll angle are adjusted. Instead of attempting to calculate the
241 // new center of rotation, pick a gear unit as a reference and use its
242 // location vector to calculate the new height change. i.e. new altitude =
243 // earth z component of that vector (which is in body axes )
244 void FGTrimAxis::SetThetaOnGround(double ff) {
247 // favor an off-center unit so that the same one can be used for both
248 // pitch and roll. An on-center unit is used (for pitch)if that's all
249 // that's in contact with the ground.
250 i=0; ref=-1; center=-1;
251 while( (ref < 0) && (i < fdmex->GetGroundReactions()->GetNumGearUnits()) ) {
252 if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetWOW()) {
253 if(fabs(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(2)) > 0.01)
260 if((ref < 0) && (center >= 0)) {
263 cout << "SetThetaOnGround ref gear: " << ref << endl;
265 double sp=fdmex->GetRotation()->GetSinphi();
266 double cp=fdmex->GetRotation()->GetCosphi();
267 double lx=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(1);
268 double ly=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(2);
269 double lz=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(3);
270 double hagl = -1*lx*sin(ff) +
274 fgic->SetAltitudeAGLFtIC(hagl);
275 cout << "SetThetaOnGround new alt: " << hagl << endl;
277 fgic->SetPitchAngleRadIC(ff);
278 cout << "SetThetaOnGround new theta: " << ff << endl;
281 /*****************************************************************************/
283 bool FGTrimAxis::initTheta(void) {
284 int i,N,iAft, iForward;
285 double zAft,zForward,zDiff,theta;
289 saveAlt=fgic->GetAltitudeAGLFtIC();
290 fgic->SetAltitudeAGLFtIC(100);
293 N=fdmex->GetGroundReactions()->GetNumGearUnits();
295 //find the first wheel unit forward of the cg
296 //the list is short so a simple linear search is fine
297 for( i=0; i<N; i++ ) {
298 if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(1) > 0 ) {
303 //now find the first wheel unit aft of the cg
304 for( i=0; i<N; i++ ) {
305 if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(1) < 0 ) {
311 // now adjust theta till the wheels are the same distance from the ground
312 zAft=fdmex->GetGroundReactions()->GetGearUnit(1)->GetLocalGear(3);
313 zForward=fdmex->GetGroundReactions()->GetGearUnit(0)->GetLocalGear(3);
314 zDiff = zForward - zAft;
316 theta=fgic->GetPitchAngleDegIC();
317 while(!level && (i < 100)) {
319 fgic->SetPitchAngleDegIC(theta);
321 zAft=fdmex->GetGroundReactions()->GetGearUnit(1)->GetLocalGear(3);
322 zForward=fdmex->GetGroundReactions()->GetGearUnit(0)->GetLocalGear(3);
323 zDiff = zForward - zAft;
324 //cout << endl << theta << " " << zDiff << endl;
325 //cout << "0: " << fdmex->GetGroundReactions()->GetGearUnit(0)->GetLocalGear() << endl;
326 //cout << "1: " << fdmex->GetGroundReactions()->GetGearUnit(1)->GetLocalGear() << endl;
328 if(fabs(zDiff ) < 0.1)
333 cout << " Initial Theta: " << fdmex->GetRotation()->Gettht()*radtodeg << endl;
334 control_min=(theta+5)*degtorad;
335 control_max=(theta-5)*degtorad;
336 fgic->SetAltitudeAGLFtIC(saveAlt);
343 /*****************************************************************************/
345 void FGTrimAxis::SetPhiOnGround(double ff) {
349 //must have an off-center unit here
350 while( (ref < 0) && (i < fdmex->GetGroundReactions()->GetNumGearUnits()) ) {
351 if( (fdmex->GetGroundReactions()->GetGearUnit(i)->GetWOW()) &&
352 (fabs(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(2)) > 0.01))
357 double st=fdmex->GetRotation()->GetSintht();
358 double ct=fdmex->GetRotation()->GetCostht();
359 double lx=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(1);
360 double ly=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(2);
361 double lz=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(3);
362 double hagl = -1*lx*st +
366 fgic->SetAltitudeAGLFtIC(hagl);
368 fgic->SetRollAngleRadIC(ff);
372 /*****************************************************************************/
374 void FGTrimAxis::Run(void) {
376 double last_state_value;
379 //cout << "FGTrimAxis::Run: " << control_value << endl;
384 last_state_value=state_value;
388 if((fabs(last_state_value - state_value) < tolerance) || (i >= 100) )
393 its_to_stable_value=i;
394 total_stability_iterations+=its_to_stable_value;
398 /*****************************************************************************/
400 void FGTrimAxis::setThrottlesPct(void) {
402 for(unsigned i=0;i<fdmex->GetPropulsion()->GetNumEngines();i++) {
403 tMin=fdmex->GetPropulsion()->GetEngine(i)->GetThrottleMin();
404 tMax=fdmex->GetPropulsion()->GetEngine(i)->GetThrottleMax();
405 //cout << "setThrottlespct: " << i << ", " << control_min << ", " << control_max << ", " << control_value;
406 fdmex->GetFCS()->SetThrottleCmd(i,tMin+control_value*(tMax-tMin));
407 //cout << "setThrottlespct: " << fdmex->GetFCS()->GetThrottleCmd(i) << endl;
408 fdmex->RunIC(fgic); //apply throttle change
409 fdmex->GetPropulsion()->GetSteadyState();
413 /*****************************************************************************/
415 void FGTrimAxis::AxisReport(void) {
418 sprintf(out," %20s: %6.2f %5s: %9.2e Tolerance: %3.0e\n",
419 GetControlName().c_str(), GetControl()*control_convert,
420 GetStateName().c_str(), GetState()+state_target, GetTolerance());
425 /*****************************************************************************/
427 double FGTrimAxis::GetAvgStability( void ) {
428 if(total_iterations > 0) {
429 return double(total_stability_iterations)/double(total_iterations);
434 /*****************************************************************************/
435 // The bitmasked value choices are as follows:
436 // unset: In this case (the default) JSBSim would only print
437 // out the normally expected messages, essentially echoing
438 // the config files as they are read. If the environment
439 // variable is not set, debug_lvl is set to 1 internally
440 // 0: This requests JSBSim not to output any messages
442 // 1: This value explicity requests the normal JSBSim
444 // 2: This value asks for a message to be printed out when
445 // a class is instantiated
446 // 4: When this value is set, a message is displayed when a
447 // FGModel object executes its Run() method
448 // 8: When this value is set, various runtime state variables
449 // are printed out periodically
450 // 16: When set various parameters are sanity checked and
451 // a message is printed out when they go out of bounds
453 void FGTrimAxis::Debug(int from)
456 if (debug_lvl <= 0) return;
457 if (debug_lvl & 1 ) { // Standard console startup message output
458 if (from == 0) { // Constructor
462 if (debug_lvl & 2 ) { // Instantiation/Destruction notification
463 if (from == 0) cout << "Instantiated: FGTrimAxis" << endl;
464 if (from == 1) cout << "Destroyed: FGTrimAxis" << endl;
466 if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
468 if (debug_lvl & 8 ) { // Runtime state variables
470 if (debug_lvl & 16) { // Sanity checking
472 if (debug_lvl & 64) {
473 if (from == 0) { // Constructor
474 cout << IdSrc << endl;
475 cout << IdHdr << endl;