]> git.mxchange.org Git - flightgear.git/blob - src/FDM/JSBSim/FGTrimAxis.cpp
Added FGTurbine.[ch]*
[flightgear.git] / src / FDM / JSBSim / FGTrimAxis.cpp
1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2  
3  Header:       FGTrimAxis.cpp
4  Author:       Tony Peden
5  Date started: 7/3/00
6  
7  --------- Copyright (C) 1999  Anthony K. Peden (apeden@earthlink.net) ---------
8  
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
12  version.
13  
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
17  details.
18  
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.
22  
23  Further information about the GNU General Public License can also be found on
24  the world wide web at http://www.gnu.org.
25  
26  
27  HISTORY
28 --------------------------------------------------------------------------------
29 7/3/00   TP   Created
30  
31 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32 INCLUDES
33 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
34
35 #ifdef _MSC_VER
36 #  pragma warning (disable : 4786)
37 #endif
38
39 #include <string>
40 #include <stdlib.h>
41
42 #include "FGFDMExec.h"
43 #include "FGAtmosphere.h"
44 #include "FGInitialCondition.h"
45 #include "FGTrimAxis.h"
46 #include "FGAircraft.h"
47 #include "FGPropulsion.h"
48
49 static const char *IdSrc = "$Id$";
50 static const char *IdHdr = ID_TRIMAXIS;
51
52 /*****************************************************************************/
53
54 FGTrimAxis::FGTrimAxis(FGFDMExec* fdex, FGInitialCondition* ic, State st,
55                        Control ctrl) {
56
57   fdmex=fdex;
58   fgic=ic;
59   state=st;
60   control=ctrl;
61   max_iterations=10;
62   control_value=0;
63   its_to_stable_value=0;
64   total_iterations=0;
65   total_stability_iterations=0;
66   state_convert=1.0;
67   control_convert=1.0;
68   state_value=0;
69   state_target=0;
70   switch(state) {
71     case tUdot: tolerance = DEFAULT_TOLERANCE; break;
72     case tVdot: tolerance = DEFAULT_TOLERANCE; break;
73     case tWdot: tolerance = DEFAULT_TOLERANCE; break;
74     case tQdot: tolerance = DEFAULT_TOLERANCE / 10; break;
75     case tPdot: tolerance = DEFAULT_TOLERANCE / 10; break;
76     case tRdot: tolerance = DEFAULT_TOLERANCE / 10; break;
77     case tHmgt: tolerance = 0.01; break;
78     case  tNlf: state_target=1.0; tolerance = 1E-5; break;
79     case tAll: break;
80   }  
81   
82   solver_eps=tolerance;
83   switch(control) {
84   case tThrottle:
85     control_min=0;
86     control_max=1;
87     control_value=0.5;
88     break;
89   case tBeta:
90     control_min=-30*degtorad;
91     control_max=30*degtorad;
92     control_convert=radtodeg;
93     break;
94   case tAlpha:
95     control_min=fdmex->GetAerodynamics()->GetAlphaCLMin();
96     control_max=fdmex->GetAerodynamics()->GetAlphaCLMax();
97     if(control_max <= control_min) {
98       control_max=20*degtorad;
99       control_min=-5*degtorad;
100     }
101     control_value= (control_min+control_max)/2;
102     control_convert=radtodeg;
103     solver_eps=tolerance/100;
104     break;
105   case tPitchTrim:
106   case tElevator:
107   case tRollTrim:
108   case tAileron:
109   case tYawTrim:
110   case tRudder:
111     control_min=-1;
112     control_max=1;
113     state_convert=radtodeg;
114     solver_eps=tolerance/100;
115     break;
116   case tAltAGL:
117     control_min=0;
118     control_max=30;
119     control_value=fdmex->GetPosition()->GetDistanceAGL();
120     solver_eps=tolerance/100;
121     break;
122   case tTheta:
123     control_min=fdmex->GetRotation()->Gettht() - 5*degtorad;
124     control_max=fdmex->GetRotation()->Gettht() + 5*degtorad;
125     state_convert=radtodeg;
126     break;
127   case tPhi:
128     control_min=fdmex->GetRotation()->Getphi() - 30*degtorad;
129     control_max=fdmex->GetRotation()->Getphi() + 30*degtorad;
130     state_convert=radtodeg;
131     control_convert=radtodeg;
132     break;
133   case tGamma:
134     solver_eps=tolerance/100;
135     control_min=-80*degtorad;
136     control_max=80*degtorad;
137     control_convert=radtodeg;
138     break;
139   case tHeading:
140     control_min=fdmex->GetRotation()->Getpsi() - 30*degtorad;
141     control_max=fdmex->GetRotation()->Getpsi() + 30*degtorad;
142     state_convert=radtodeg;
143     break;
144   }
145   
146   
147   Debug(0);
148 }
149
150 /*****************************************************************************/
151
152 FGTrimAxis::~FGTrimAxis(void)
153 {
154   Debug(1);
155 }
156
157 /*****************************************************************************/
158
159 void FGTrimAxis::getState(void) {
160   switch(state) {
161   case tUdot: state_value=fdmex->GetTranslation()->GetUVWdot(1)-state_target; break;
162   case tVdot: state_value=fdmex->GetTranslation()->GetUVWdot(2)-state_target; break;
163   case tWdot: state_value=fdmex->GetTranslation()->GetUVWdot(3)-state_target; break;
164   case tQdot: state_value=fdmex->GetRotation()->GetPQRdot(2)-state_target;break;
165   case tPdot: state_value=fdmex->GetRotation()->GetPQRdot(1)-state_target; break;
166   case tRdot: state_value=fdmex->GetRotation()->GetPQRdot(3)-state_target; break;
167   case tHmgt: state_value=computeHmgt()-state_target; break;
168   case tNlf:  state_value=fdmex->GetAircraft()->GetNlf()-state_target; break;
169   case tAll: break;
170   }
171 }
172
173 /*****************************************************************************/
174
175 //States are not settable
176
177 void FGTrimAxis::getControl(void) {
178   switch(control) {
179   case tThrottle:  control_value=fdmex->GetFCS()->GetThrottleCmd(0); break;
180   case tBeta:      control_value=fdmex->GetTranslation()->Getalpha(); break;
181   case tAlpha:     control_value=fdmex->GetTranslation()->Getbeta();  break;
182   case tPitchTrim: control_value=fdmex->GetFCS() -> GetPitchTrimCmd(); break;
183   case tElevator:  control_value=fdmex->GetFCS() -> GetDeCmd(); break;
184   case tRollTrim:
185   case tAileron:   control_value=fdmex->GetFCS() -> GetDaCmd(); break;
186   case tYawTrim:
187   case tRudder:    control_value=fdmex->GetFCS() -> GetDrCmd(); break;
188   case tAltAGL:    control_value=fdmex->GetPosition()->GetDistanceAGL();break;
189   case tTheta:     control_value=fdmex->GetRotation()->Gettht(); break;
190   case tPhi:       control_value=fdmex->GetRotation()->Getphi(); break;
191   case tGamma:     control_value=fdmex->GetPosition()->GetGamma();break;
192   case tHeading:   control_value=fdmex->GetRotation()->Getpsi(); break;
193   }
194 }
195
196 /*****************************************************************************/
197
198 double FGTrimAxis::computeHmgt(void) {
199   double diff;
200   
201   diff   = fdmex->GetRotation()->Getpsi() - 
202              fdmex->GetPosition()->GetGroundTrack();
203   
204   if( diff < -M_PI ) {
205      return (diff + 2*M_PI);
206   } else if( diff > M_PI ) {
207      return (diff - 2*M_PI);
208   } else {
209      return diff;
210   }
211
212 }
213        
214 /*****************************************************************************/
215
216
217 void FGTrimAxis::setControl(void) {
218   switch(control) {
219   case tThrottle:  setThrottlesPct(); break;
220   case tBeta:      fgic->SetBetaRadIC(control_value); break;
221   case tAlpha:     fgic->SetAlphaRadIC(control_value);  break;
222   case tPitchTrim: fdmex->GetFCS()->SetPitchTrimCmd(control_value); break;
223   case tElevator:  fdmex->GetFCS()->SetDeCmd(control_value); break;
224   case tRollTrim:
225   case tAileron:   fdmex->GetFCS()->SetDaCmd(control_value); break;
226   case tYawTrim:
227   case tRudder:    fdmex->GetFCS()->SetDrCmd(control_value); break;
228   case tAltAGL:    fgic->SetAltitudeAGLFtIC(control_value); break;
229   case tTheta:     fgic->SetPitchAngleRadIC(control_value); break;
230   case tPhi:       fgic->SetRollAngleRadIC(control_value); break;
231   case tGamma:     fgic->SetFlightPathAngleRadIC(control_value); break;
232   case tHeading:   fgic->SetTrueHeadingRadIC(control_value); break;
233   }
234 }
235
236
237   
238
239
240 /*****************************************************************************/
241
242 // the aircraft center of rotation is no longer the cg once the gear
243 // contact the ground so the altitude needs to be changed when pitch 
244 // and roll angle are adjusted.  Instead of attempting to calculate the 
245 // new center of rotation, pick a gear unit as a reference and use its
246 // location vector to calculate the new height change. i.e. new altitude =
247 // earth z component of that vector (which is in body axes )  
248 void FGTrimAxis::SetThetaOnGround(double ff) {
249   int center,i,ref;
250
251   // favor an off-center unit so that the same one can be used for both
252   // pitch and roll.  An on-center unit is used (for pitch)if that's all 
253   // that's in contact with the ground.
254   i=0; ref=-1; center=-1;
255   while( (ref < 0) && (i < fdmex->GetGroundReactions()->GetNumGearUnits()) ) {
256     if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetWOW()) {
257       if(fabs(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(2)) > 0.01)
258         ref=i;
259       else
260         center=i;
261     } 
262     i++; 
263   }
264   if((ref < 0) && (center >= 0)) {
265     ref=center;
266   }
267   cout << "SetThetaOnGround ref gear: " << ref << endl;
268   if(ref >= 0) {
269     double sp=fdmex->GetRotation()->GetSinphi();
270     double cp=fdmex->GetRotation()->GetCosphi();
271     double lx=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(1);
272     double ly=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(2);
273     double lz=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(3);
274     double hagl = -1*lx*sin(ff) +
275                     ly*sp*cos(ff) +
276                     lz*cp*cos(ff);
277    
278     fgic->SetAltitudeAGLFtIC(hagl);
279     cout << "SetThetaOnGround new alt: " << hagl << endl;
280   }                   
281   fgic->SetPitchAngleRadIC(ff);  
282   cout << "SetThetaOnGround new theta: " << ff << endl;      
283 }      
284
285 /*****************************************************************************/
286
287 bool FGTrimAxis::initTheta(void) {
288   int i,N,iAft, iForward;
289   double zAft,zForward,zDiff,theta;
290   bool level;  
291   double saveAlt;
292   
293   saveAlt=fgic->GetAltitudeAGLFtIC();
294   fgic->SetAltitudeAGLFtIC(100);
295   
296   
297   N=fdmex->GetGroundReactions()->GetNumGearUnits();
298   
299   //find the first wheel unit forward of the cg
300   //the list is short so a simple linear search is fine
301   for( i=0; i<N; i++ ) {
302     if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(1) > 0 ) {
303         iForward=i;
304         break;
305     }
306   }
307   //now find the first wheel unit aft of the cg
308   for( i=0; i<N; i++ ) {
309     if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(1) < 0 ) {
310         iAft=i;
311         break;
312     }
313   }
314           
315   // now adjust theta till the wheels are the same distance from the ground
316   zAft=fdmex->GetGroundReactions()->GetGearUnit(1)->GetLocalGear(3);
317   zForward=fdmex->GetGroundReactions()->GetGearUnit(0)->GetLocalGear(3);
318   zDiff = zForward - zAft;
319   level=false;
320   theta=fgic->GetPitchAngleDegIC(); 
321   while(!level && (i < 100)) {
322         theta+=2.0*zDiff;
323         fgic->SetPitchAngleDegIC(theta);   
324         fdmex->RunIC(fgic);
325         zAft=fdmex->GetGroundReactions()->GetGearUnit(1)->GetLocalGear(3);
326         zForward=fdmex->GetGroundReactions()->GetGearUnit(0)->GetLocalGear(3);
327         zDiff = zForward - zAft;
328         //cout << endl << theta << "  " << zDiff << endl;
329         //cout << "0: " << fdmex->GetGroundReactions()->GetGearUnit(0)->GetLocalGear() << endl;
330         //cout << "1: " << fdmex->GetGroundReactions()->GetGearUnit(1)->GetLocalGear() << endl;
331
332         if(fabs(zDiff ) < 0.1) 
333             level=true;
334         i++;   
335   }                     
336   //cout << i << endl;
337   cout << "    Initial Theta: " << fdmex->GetRotation()->Gettht()*radtodeg << endl;
338   control_min=(theta+5)*degtorad;
339   control_max=(theta-5)*degtorad;
340   fgic->SetAltitudeAGLFtIC(saveAlt);
341   if(i < 100) 
342     return true;
343   else
344     return false;  
345
346
347 /*****************************************************************************/
348
349 void FGTrimAxis::SetPhiOnGround(double ff) {
350   int i,ref;
351
352   i=0; ref=-1;
353   //must have an off-center unit here 
354   while( (ref < 0) && (i < fdmex->GetGroundReactions()->GetNumGearUnits()) ) {
355     if( (fdmex->GetGroundReactions()->GetGearUnit(i)->GetWOW()) && 
356       (fabs(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(2)) > 0.01))
357         ref=i;
358     i++; 
359   }
360   if(ref >= 0) {
361     double st=fdmex->GetRotation()->GetSintht();
362     double ct=fdmex->GetRotation()->GetCostht();
363     double lx=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(1);
364     double ly=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(2);
365     double lz=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(3);
366     double hagl = -1*lx*st +
367                     ly*sin(ff)*ct +
368                     lz*cos(ff)*ct;
369    
370     fgic->SetAltitudeAGLFtIC(hagl);
371   }                   
372   fgic->SetRollAngleRadIC(ff);
373
374 }      
375
376 /*****************************************************************************/
377
378 void FGTrimAxis::Run(void) {
379
380   double last_state_value;
381   int i;
382   setControl();
383   //cout << "FGTrimAxis::Run: " << control_value << endl;
384   i=0;
385   bool stable=false;
386   while(!stable) {
387     i++;
388     last_state_value=state_value;
389     fdmex->RunIC(fgic);
390     getState();
391     if(i > 1) {
392       if((fabs(last_state_value - state_value) < tolerance) || (i >= 100) )
393         stable=true;
394     }
395   }
396
397   its_to_stable_value=i;
398   total_stability_iterations+=its_to_stable_value;
399   total_iterations++;
400 }
401
402 /*****************************************************************************/
403
404 void FGTrimAxis::setThrottlesPct(void) {
405   double tMin,tMax;
406   for(unsigned i=0;i<fdmex->GetPropulsion()->GetNumEngines();i++) {
407       tMin=fdmex->GetPropulsion()->GetEngine(i)->GetThrottleMin();
408       tMax=fdmex->GetPropulsion()->GetEngine(i)->GetThrottleMax();
409       //cout << "setThrottlespct: " << i << ", " << control_min << ", " << control_max << ", " << control_value;
410       fdmex->GetFCS()->SetThrottleCmd(i,tMin+control_value*(tMax-tMin));
411       //cout << "setThrottlespct: " << fdmex->GetFCS()->GetThrottleCmd(i) << endl;
412       fdmex->RunIC(fgic); //apply throttle change
413       fdmex->GetPropulsion()->GetSteadyState();
414   }
415 }
416
417 /*****************************************************************************/
418
419 void FGTrimAxis::AxisReport(void) {
420   
421   char out[80];
422   sprintf(out,"  %20s: %6.2f %5s: %9.2e Tolerance: %3.0e\n",
423            GetControlName().c_str(), GetControl()*control_convert,
424            GetStateName().c_str(), GetState()+state_target, GetTolerance()); 
425   cout << out;
426
427 }
428
429 /*****************************************************************************/
430
431 double FGTrimAxis::GetAvgStability( void ) {
432   if(total_iterations > 0) {
433     return double(total_stability_iterations)/double(total_iterations);
434   }
435   return 0;
436 }
437
438 /*****************************************************************************/
439 //    The bitmasked value choices are as follows:
440 //    unset: In this case (the default) JSBSim would only print
441 //       out the normally expected messages, essentially echoing
442 //       the config files as they are read. If the environment
443 //       variable is not set, debug_lvl is set to 1 internally
444 //    0: This requests JSBSim not to output any messages
445 //       whatsoever.
446 //    1: This value explicity requests the normal JSBSim
447 //       startup messages
448 //    2: This value asks for a message to be printed out when
449 //       a class is instantiated
450 //    4: When this value is set, a message is displayed when a
451 //       FGModel object executes its Run() method
452 //    8: When this value is set, various runtime state variables
453 //       are printed out periodically
454 //    16: When set various parameters are sanity checked and
455 //       a message is printed out when they go out of bounds
456
457 void FGTrimAxis::Debug(int from)
458 {
459
460   if (debug_lvl <= 0) return;
461   if (debug_lvl & 1 ) { // Standard console startup message output
462     if (from == 0) { // Constructor
463
464     }
465   }
466   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
467     if (from == 0) cout << "Instantiated: FGTrimAxis" << endl;
468     if (from == 1) cout << "Destroyed:    FGTrimAxis" << endl;
469   }
470   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
471   }
472   if (debug_lvl & 8 ) { // Runtime state variables
473   }
474   if (debug_lvl & 16) { // Sanity checking
475   }
476   if (debug_lvl & 64) {
477     if (from == 0) { // Constructor
478       cout << IdSrc << endl;
479       cout << IdHdr << endl;
480     }
481   }
482 }
483
484