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 Lesser 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 Lesser General Public License for more
19 You should have received a copy of the GNU Lesser 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 Lesser General Public License can also be found on
24 the world wide web at http://www.gnu.org.
28 --------------------------------------------------------------------------------
31 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
36 # pragma warning (disable : 4786)
42 #include "FGFDMExec.h"
43 #include "models/FGAtmosphere.h"
44 #include "FGInitialCondition.h"
45 #include "FGTrimAxis.h"
46 #include "models/FGAircraft.h"
47 #include "models/FGPropulsion.h"
48 #include "models/FGAerodynamics.h"
49 #include "models/FGFCS.h"
50 #include "models/propulsion/FGEngine.h"
51 #include "models/FGAuxiliary.h"
52 #include "models/FGGroundReactions.h"
58 static const char *IdSrc = "$Id$";
59 static const char *IdHdr = ID_TRIMAXIS;
61 /*****************************************************************************/
63 FGTrimAxis::FGTrimAxis(FGFDMExec* fdex, FGInitialCondition* ic, State st,
72 its_to_stable_value=0;
74 total_stability_iterations=0;
80 case tUdot: tolerance = DEFAULT_TOLERANCE; break;
81 case tVdot: tolerance = DEFAULT_TOLERANCE; break;
82 case tWdot: tolerance = DEFAULT_TOLERANCE; break;
83 case tQdot: tolerance = DEFAULT_TOLERANCE / 10; break;
84 case tPdot: tolerance = DEFAULT_TOLERANCE / 10; break;
85 case tRdot: tolerance = DEFAULT_TOLERANCE / 10; break;
86 case tHmgt: tolerance = 0.01; break;
87 case tNlf: state_target=1.0; tolerance = 1E-5; break;
99 control_min=-30*degtorad;
100 control_max=30*degtorad;
101 control_convert=radtodeg;
104 control_min=fdmex->GetAerodynamics()->GetAlphaCLMin();
105 control_max=fdmex->GetAerodynamics()->GetAlphaCLMax();
106 if(control_max <= control_min) {
107 control_max=20*degtorad;
108 control_min=-5*degtorad;
110 control_value= (control_min+control_max)/2;
111 control_convert=radtodeg;
112 solver_eps=tolerance/100;
122 state_convert=radtodeg;
123 solver_eps=tolerance/100;
128 control_value=fdmex->GetPropagate()->GetDistanceAGL();
129 solver_eps=tolerance/100;
132 control_min=fdmex->GetPropagate()->GetEuler(eTht) - 5*degtorad;
133 control_max=fdmex->GetPropagate()->GetEuler(eTht) + 5*degtorad;
134 state_convert=radtodeg;
137 control_min=fdmex->GetPropagate()->GetEuler(ePhi) - 30*degtorad;
138 control_max=fdmex->GetPropagate()->GetEuler(ePhi) + 30*degtorad;
139 state_convert=radtodeg;
140 control_convert=radtodeg;
143 solver_eps=tolerance/100;
144 control_min=-80*degtorad;
145 control_max=80*degtorad;
146 control_convert=radtodeg;
149 control_min=fdmex->GetPropagate()->GetEuler(ePsi) - 30*degtorad;
150 control_max=fdmex->GetPropagate()->GetEuler(ePsi) + 30*degtorad;
151 state_convert=radtodeg;
159 /*****************************************************************************/
161 FGTrimAxis::~FGTrimAxis(void)
166 /*****************************************************************************/
168 void FGTrimAxis::getState(void) {
170 case tUdot: state_value=fdmex->GetPropagate()->GetUVWdot(1)-state_target; break;
171 case tVdot: state_value=fdmex->GetPropagate()->GetUVWdot(2)-state_target; break;
172 case tWdot: state_value=fdmex->GetPropagate()->GetUVWdot(3)-state_target; break;
173 case tQdot: state_value=fdmex->GetPropagate()->GetPQRdot(2)-state_target;break;
174 case tPdot: state_value=fdmex->GetPropagate()->GetPQRdot(1)-state_target; break;
175 case tRdot: state_value=fdmex->GetPropagate()->GetPQRdot(3)-state_target; break;
176 case tHmgt: state_value=computeHmgt()-state_target; break;
177 case tNlf: state_value=fdmex->GetAircraft()->GetNlf()-state_target; break;
182 /*****************************************************************************/
184 //States are not settable
186 void FGTrimAxis::getControl(void) {
188 case tThrottle: control_value=fdmex->GetFCS()->GetThrottleCmd(0); break;
189 case tBeta: control_value=fdmex->GetAuxiliary()->Getbeta(); break;
190 case tAlpha: control_value=fdmex->GetAuxiliary()->Getalpha(); break;
191 case tPitchTrim: control_value=fdmex->GetFCS() -> GetPitchTrimCmd(); break;
192 case tElevator: control_value=fdmex->GetFCS() -> GetDeCmd(); break;
194 case tAileron: control_value=fdmex->GetFCS() -> GetDaCmd(); break;
196 case tRudder: control_value=fdmex->GetFCS() -> GetDrCmd(); break;
197 case tAltAGL: control_value=fdmex->GetPropagate()->GetDistanceAGL();break;
198 case tTheta: control_value=fdmex->GetPropagate()->GetEuler(eTht); break;
199 case tPhi: control_value=fdmex->GetPropagate()->GetEuler(ePhi); break;
200 case tGamma: control_value=fdmex->GetAuxiliary()->GetGamma();break;
201 case tHeading: control_value=fdmex->GetPropagate()->GetEuler(ePsi); break;
205 /*****************************************************************************/
207 double FGTrimAxis::computeHmgt(void) {
210 diff = fdmex->GetPropagate()->GetEuler(ePsi) -
211 fdmex->GetAuxiliary()->GetGroundTrack();
214 return (diff + 2*M_PI);
215 } else if( diff > M_PI ) {
216 return (diff - 2*M_PI);
223 /*****************************************************************************/
226 void FGTrimAxis::setControl(void) {
228 case tThrottle: setThrottlesPct(); break;
229 case tBeta: fgic->SetBetaRadIC(control_value); break;
230 case tAlpha: fgic->SetAlphaRadIC(control_value); break;
231 case tPitchTrim: fdmex->GetFCS()->SetPitchTrimCmd(control_value); break;
232 case tElevator: fdmex->GetFCS()->SetDeCmd(control_value); break;
234 case tAileron: fdmex->GetFCS()->SetDaCmd(control_value); break;
236 case tRudder: fdmex->GetFCS()->SetDrCmd(control_value); break;
237 case tAltAGL: fgic->SetAltitudeAGLFtIC(control_value); break;
238 case tTheta: fgic->SetThetaRadIC(control_value); break;
239 case tPhi: fgic->SetPhiRadIC(control_value); break;
240 case tGamma: fgic->SetFlightPathAngleRadIC(control_value); break;
241 case tHeading: fgic->SetPsiRadIC(control_value); break;
249 /*****************************************************************************/
251 // the aircraft center of rotation is no longer the cg once the gear
252 // contact the ground so the altitude needs to be changed when pitch
253 // and roll angle are adjusted. Instead of attempting to calculate the
254 // new center of rotation, pick a gear unit as a reference and use its
255 // location vector to calculate the new height change. i.e. new altitude =
256 // earth z component of that vector (which is in body axes )
257 void FGTrimAxis::SetThetaOnGround(double ff) {
260 // favor an off-center unit so that the same one can be used for both
261 // pitch and roll. An on-center unit is used (for pitch)if that's all
262 // that's in contact with the ground.
263 i=0; ref=-1; center=-1;
264 while( (ref < 0) && (i < fdmex->GetGroundReactions()->GetNumGearUnits()) ) {
265 if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetWOW()) {
266 if(fabs(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(2)) > 0.01)
273 if((ref < 0) && (center >= 0)) {
276 cout << "SetThetaOnGround ref gear: " << ref << endl;
278 double sp = fdmex->GetPropagate()->GetSinEuler(ePhi);
279 double cp = fdmex->GetPropagate()->GetCosEuler(ePhi);
280 double lx = fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(1);
281 double ly = fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(2);
282 double lz = fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(3);
283 double hagl = -1*lx*sin(ff) +
287 fgic->SetAltitudeAGLFtIC(hagl);
288 cout << "SetThetaOnGround new alt: " << hagl << endl;
290 fgic->SetThetaRadIC(ff);
291 cout << "SetThetaOnGround new theta: " << ff << endl;
294 /*****************************************************************************/
296 bool FGTrimAxis::initTheta(void) {
300 double zAft,zForward,zDiff,theta;
301 double xAft,xForward,xDiff;
305 saveAlt=fgic->GetAltitudeAGLFtIC();
306 fgic->SetAltitudeAGLFtIC(100);
309 N=fdmex->GetGroundReactions()->GetNumGearUnits();
311 //find the first wheel unit forward of the cg
312 //the list is short so a simple linear search is fine
313 for( i=0; i<N; i++ ) {
314 if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(1) > 0 ) {
319 //now find the first wheel unit aft of the cg
320 for( i=0; i<N; i++ ) {
321 if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(1) < 0 ) {
327 // now adjust theta till the wheels are the same distance from the ground
328 xAft=fdmex->GetGroundReactions()->GetGearUnit(iAft)->GetBodyLocation(1);
329 xForward=fdmex->GetGroundReactions()->GetGearUnit(iForward)->GetBodyLocation(1);
330 xDiff = xForward - xAft;
331 zAft=fdmex->GetGroundReactions()->GetGearUnit(iAft)->GetLocalGear(3);
332 zForward=fdmex->GetGroundReactions()->GetGearUnit(iForward)->GetLocalGear(3);
333 zDiff = zForward - zAft;
335 theta=fgic->GetThetaDegIC();
336 while(!level && (i < 100)) {
337 theta+=radtodeg*atan(zDiff/xDiff);
338 fgic->SetThetaDegIC(theta);
340 zAft=fdmex->GetGroundReactions()->GetGearUnit(iAft)->GetLocalGear(3);
341 zForward=fdmex->GetGroundReactions()->GetGearUnit(iForward)->GetLocalGear(3);
342 zDiff = zForward - zAft;
343 //cout << endl << theta << " " << zDiff << endl;
344 //cout << "0: " << fdmex->GetGroundReactions()->GetGearUnit(0)->GetLocalGear() << endl;
345 //cout << "1: " << fdmex->GetGroundReactions()->GetGearUnit(1)->GetLocalGear() << endl;
346 if(fabs(zDiff ) < 0.1)
352 cout << " Initial Theta: " << fdmex->GetPropagate()->GetEuler(eTht)*radtodeg << endl;
353 cout << " Used gear unit " << iAft << " as aft and " << iForward << " as forward" << endl;
355 control_min=(theta+5)*degtorad;
356 control_max=(theta-5)*degtorad;
357 fgic->SetAltitudeAGLFtIC(saveAlt);
364 /*****************************************************************************/
366 void FGTrimAxis::SetPhiOnGround(double ff) {
370 //must have an off-center unit here
371 while ( (ref < 0) && (i < fdmex->GetGroundReactions()->GetNumGearUnits()) ) {
372 if ( (fdmex->GetGroundReactions()->GetGearUnit(i)->GetWOW()) &&
373 (fabs(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(2)) > 0.01))
378 double st = fdmex->GetPropagate()->GetSinEuler(eTht);
379 double ct = fdmex->GetPropagate()->GetCosEuler(eTht);
380 double lx = fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(1);
381 double ly = fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(2);
382 double lz = fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(3);
383 double hagl = -1*lx*st +
387 fgic->SetAltitudeAGLFtIC(hagl);
389 fgic->SetPhiRadIC(ff);
393 /*****************************************************************************/
395 void FGTrimAxis::Run(void) {
397 double last_state_value;
400 //cout << "FGTrimAxis::Run: " << control_value << endl;
405 last_state_value=state_value;
409 if((fabs(last_state_value - state_value) < tolerance) || (i >= 100) )
414 its_to_stable_value=i;
415 total_stability_iterations+=its_to_stable_value;
419 /*****************************************************************************/
421 void FGTrimAxis::setThrottlesPct(void) {
423 for(unsigned i=0;i<fdmex->GetPropulsion()->GetNumEngines();i++) {
424 tMin=fdmex->GetPropulsion()->GetEngine(i)->GetThrottleMin();
425 tMax=fdmex->GetPropulsion()->GetEngine(i)->GetThrottleMax();
426 //cout << "setThrottlespct: " << i << ", " << control_min << ", " << control_max << ", " << control_value;
427 fdmex->GetFCS()->SetThrottleCmd(i,tMin+control_value*(tMax-tMin));
428 //cout << "setThrottlespct: " << fdmex->GetFCS()->GetThrottleCmd(i) << endl;
429 fdmex->RunIC(); //apply throttle change
430 fdmex->GetPropulsion()->GetSteadyState();
434 /*****************************************************************************/
436 void FGTrimAxis::AxisReport(void) {
437 cout << " " << setw(20) << GetControlName() << ": ";
438 cout << setw(6) << setprecision(2) << GetControl()*control_convert << ' ';
439 cout << setw(5) << GetStateName() << ": ";
440 cout << setw(9) << setprecision(2) << scientific << GetState()+state_target;
441 cout << " Tolerance: " << setw(3) << setprecision(0) << scientific << GetTolerance();
443 if( fabs(GetState()+state_target) < fabs(GetTolerance()) )
444 cout << " Passed" << endl;
446 cout << " Failed" << endl;
449 /*****************************************************************************/
451 double FGTrimAxis::GetAvgStability( void ) {
452 if(total_iterations > 0) {
453 return double(total_stability_iterations)/double(total_iterations);
458 /*****************************************************************************/
459 // The bitmasked value choices are as follows:
460 // unset: In this case (the default) JSBSim would only print
461 // out the normally expected messages, essentially echoing
462 // the config files as they are read. If the environment
463 // variable is not set, debug_lvl is set to 1 internally
464 // 0: This requests JSBSim not to output any messages
466 // 1: This value explicity requests the normal JSBSim
468 // 2: This value asks for a message to be printed out when
469 // a class is instantiated
470 // 4: When this value is set, a message is displayed when a
471 // FGModel object executes its Run() method
472 // 8: When this value is set, various runtime state variables
473 // are printed out periodically
474 // 16: When set various parameters are sanity checked and
475 // a message is printed out when they go out of bounds
477 void FGTrimAxis::Debug(int from)
480 if (debug_lvl <= 0) return;
481 if (debug_lvl & 1 ) { // Standard console startup message output
482 if (from == 0) { // Constructor
486 if (debug_lvl & 2 ) { // Instantiation/Destruction notification
487 if (from == 0) cout << "Instantiated: FGTrimAxis" << endl;
488 if (from == 1) cout << "Destroyed: FGTrimAxis" << endl;
490 if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
492 if (debug_lvl & 8 ) { // Runtime state variables
494 if (debug_lvl & 16) { // Sanity checking
496 if (debug_lvl & 64) {
497 if (from == 0) { // Constructor
498 cout << IdSrc << endl;
499 cout << IdHdr << endl;