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