]> git.mxchange.org Git - flightgear.git/blobdiff - src/FDM/JSBSim/FGTrim.cpp
Updated to match changes in radiostack.[ch]xx
[flightgear.git] / src / FDM / JSBSim / FGTrim.cpp
index 55ed69a0d9edf42223549e1afecefedf2371c74c..e4a5c2ea658f412a824859a11780eec0637b8157 100644 (file)
@@ -53,6 +53,7 @@ INCLUDES
 #include "FGAircraft.h"
 #include "FGMassBalance.h"
 #include "FGAerodynamics.h"
+#include "FGColumnVector3.h"
 #if _MSC_VER
 #pragma warning (disable : 4786 4788)
 #endif
@@ -70,7 +71,7 @@ FGTrim::FGTrim(FGFDMExec *FDMExec,FGInitialCondition *FGIC, TrimMode tt ) {
   Tolerance=1E-3;
   A_Tolerance = Tolerance / 10;
   
-  Debug=0;
+  Debug=0;DebugLevel=0;
   fdmex=FDMExec;
   fgic=FGIC;
   total_its=0;
@@ -79,6 +80,8 @@ FGTrim::FGTrim(FGFDMExec *FDMExec,FGInitialCondition *FGIC, TrimMode tt ) {
   axis_count=0;
   mode=tt;
   xlo=xhi=alo=ahi;
+  targetNlf=1.0;
+  debug_axis=tAll;
   switch(mode) {
   case tFull:
     cout << "  Full Trim" << endl;
@@ -102,10 +105,30 @@ FGTrim::FGTrim(FGFDMExec *FDMExec,FGInitialCondition *FGIC, TrimMode tt ) {
     TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tQdot,tTheta ));
     //TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tPdot,tPhi ));
     break;
-  }
+  case tPullup:
+    TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tNlf,tAlpha ));
+    TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tUdot,tThrottle ));
+    TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tQdot,tPitchTrim ));
+    TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tHmgt,tBeta ));
+    TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tVdot,tPhi ));
+    TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tPdot,tAileron ));
+    TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tRdot,tRudder ));
+    break;
+  case tTurn:
+    TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tWdot,tAlpha ));
+    TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tUdot,tThrottle ));
+    TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tQdot,tPitchTrim ));
+    TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tVdot,tBeta ));
+    TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tPdot,tAileron ));
+    TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tRdot,tRudder ));
+    break;
+  case tCustom:
+  case tNone:
+    break;
+}
   //cout << "TrimAxes.size(): " << TrimAxes.size() << endl;
-  sub_iterations=new float[TrimAxes.size()];
-  successful=new float[TrimAxes.size()];
+  sub_iterations=new double[TrimAxes.size()];
+  successful=new double[TrimAxes.size()];
   solution=new bool[TrimAxes.size()];
   current_axis=0;
   
@@ -138,7 +161,7 @@ void FGTrim::TrimStats() {
       snprintf(out,80,"   %5s: %3.0f average: %5.2f  successful: %3.0f  stability: %5.2f\n",
                   TrimAxes[current_axis]->GetStateName().c_str(),
                   sub_iterations[current_axis],
-                  sub_iterations[current_axis]/float(total_its),
+                  sub_iterations[current_axis]/double(total_its),
                   successful[current_axis],
                   TrimAxes[current_axis]->GetAvgStability() );
       cout << out;
@@ -158,76 +181,6 @@ void FGTrim::Report(void) {
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-void FGTrim::ReportState(void) {
-  char out[80], flap[10], gear[10];
-  
-  cout << endl << "  JSBSim State" << endl;
-  snprintf(out,80,"    Weight: %7.0f lbs.  CG: %5.1f, %5.1f, %5.1f inches\n",
-                   fdmex->GetMassBalance()->GetWeight(),
-                   fdmex->GetMassBalance()->GetXYZcg(1),
-                   fdmex->GetMassBalance()->GetXYZcg(2),
-                   fdmex->GetMassBalance()->GetXYZcg(3));
-  cout << out;             
-  if( fdmex->GetFCS()->GetDfPos() <= 0.01)
-    snprintf(flap,10,"Up");
-  else
-    snprintf(flap,10,"%2.0f",fdmex->GetFCS()->GetDfPos());
-  if(fdmex->GetAircraft()->GetGearUp() == true)
-    snprintf(gear,10,"Up");
-  else
-    snprintf(gear,10,"Down");
-  snprintf(out,80, "    Flaps: %3s  Gear: %4s\n",flap,gear);
-  cout << out;
-  snprintf(out,80, "    Speed: %4.0f KCAS  Mach: %5.2f\n",
-                    fdmex->GetAuxiliary()->GetVcalibratedKTS(),
-                    fdmex->GetState()->GetParameter(FG_MACH) );
-  cout << out;
-  snprintf(out,80, "    Altitude: %7.0f ft.  AGL Altitude: %7.0f ft.\n",
-                    fdmex->GetPosition()->Geth(),
-                    fdmex->GetPosition()->GetDistanceAGL() );
-  cout << out;
-  snprintf(out,80, "    Angle of Attack: %6.2f deg  Pitch Angle: %6.2f deg\n",
-                    fdmex->GetState()->GetParameter(FG_ALPHA)*RADTODEG,
-                    fdmex->GetRotation()->Gettht()*RADTODEG );
-  cout << out;
-  snprintf(out,80, "    Flight Path Angle: %6.2f deg  Climb Rate: %5.0f ft/min\n",
-                    fdmex->GetPosition()->GetGamma()*RADTODEG,
-                    fdmex->GetPosition()->Gethdot()*60 );
-  cout << out;                  
-  snprintf(out,80, "    Normal Load Factor: %4.2f g's  Pitch Rate: %5.2f deg/s\n",
-                    fdmex->GetAerodynamics()->GetNlf(),
-                    fdmex->GetState()->GetParameter(FG_PITCHRATE)*RADTODEG );
-  cout << out;
-  snprintf(out,80, "    Heading: %3.0f deg true  Sideslip: %5.2f deg\n",
-                    fdmex->GetRotation()->Getpsi()*RADTODEG,
-                    fdmex->GetState()->GetParameter(FG_BETA)*RADTODEG );                  
-  cout << out;
-  snprintf(out,80, "    Bank Angle: %5.2f deg\n",
-                    fdmex->GetRotation()->Getphi()*RADTODEG );
-  cout << out;
-  snprintf(out,80, "    Elevator: %5.2f deg  Left Aileron: %5.2f deg  Rudder: %5.2f deg\n",
-                    fdmex->GetState()->GetParameter(FG_ELEVATOR_POS)*RADTODEG,
-                    fdmex->GetState()->GetParameter(FG_AILERON_POS)*RADTODEG,
-                    fdmex->GetState()->GetParameter(FG_RUDDER_POS)*RADTODEG );
-  cout << out;                  
-  snprintf(out,80, "    Throttle: %5.2f%c\n",
-                    fdmex->GetFCS()->GetThrottlePos(0),'%' );
-  cout << out;
-  
-  snprintf(out,80, "    Wind Components: %5.2f kts head wind, %5.2f kts cross wind\n",
-                    fdmex->GetAuxiliary()->GetHeadWind()*jsbFPSTOKTS,
-                    fdmex->GetAuxiliary()->GetCrossWind()*jsbFPSTOKTS );
-  cout << out; 
-  
-  snprintf(out,80, "    Ground Speed: %4.0f knots , Ground Track: %3.0f deg true\n",
-                    fdmex->GetPosition()->GetVground()*jsbFPSTOKTS,
-                    fdmex->GetPosition()->GetGroundTrack()*RADTODEG );
-  cout << out;                                   
-
-}                  
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 void FGTrim::ClearStates(void) {
     FGTrimAxis* ta;
     
@@ -262,8 +215,8 @@ bool FGTrim::AddState( State state, Control control ) {
     delete[] sub_iterations;
     delete[] successful;
     delete[] solution;
-    sub_iterations=new float[TrimAxes.size()];
-    successful=new float[TrimAxes.size()];
+    sub_iterations=new double[TrimAxes.size()];
+    successful=new double[TrimAxes.size()];
     solution=new bool[TrimAxes.size()];
   }
   return result;
@@ -291,8 +244,8 @@ bool FGTrim::RemoveState( State state ) {
     delete[] sub_iterations;
     delete[] successful;
     delete[] solution;
-    sub_iterations=new float[TrimAxes.size()];
-    successful=new float[TrimAxes.size()];
+    sub_iterations=new double[TrimAxes.size()];
+    successful=new double[TrimAxes.size()];
     solution=new bool[TrimAxes.size()];
   }  
   return result;
@@ -317,7 +270,6 @@ bool FGTrim::EditState( State state, Control new_control ){
       }
       iAxes++;
   }
-  cout << "Exit FGTrim::EditState(...)" << endl;
   return result;
 }  
        
@@ -329,8 +281,8 @@ bool FGTrim::DoTrim(void) {
   trim_failed=false;
   int i;
 
-  for(i=0;i < fdmex->GetAircraft()->GetNumGearUnits();i++){
-    fdmex->GetAircraft()->GetGearUnit(i)->SetReport(false);
+  for(i=0;i < fdmex->GetGroundReactions()->GetNumGearUnits();i++){
+    fdmex->GetGroundReactions()->GetGearUnit(i)->SetReport(false);
   }
 
   fdmex->GetOutput()->Disable();
@@ -352,9 +304,24 @@ bool FGTrim::DoTrim(void) {
     successful[current_axis]=0;
     solution[current_axis]=false;
   }
+  
+  
+  if(mode == tPullup ) {
+    cout << "Setting pitch rate and nlf... " << endl;
+    setupPullup();
+    cout << "pitch rate done ... " << endl;
+    TrimAxes[0]->SetStateTarget(targetNlf);
+    cout << "nlf done" << endl;
+  } else if (mode == tTurn) {
+    setupTurn();
+    //TrimAxes[0]->SetStateTarget(targetNlf);
+  }  
+  
   do {
     axis_count=0;
     for(current_axis=0;current_axis<TrimAxes.size();current_axis++) {
+      setDebug();
+      updateRates();
       Nsub=0;
       if(!solution[current_axis]) {
         if(checkLimits()) { 
@@ -425,8 +392,8 @@ bool FGTrim::DoTrim(void) {
     total_its=N;
     cout << endl << "  Trim failed" << endl;
   }
-  for(i=0;i < fdmex->GetAircraft()->GetNumGearUnits();i++){
-    fdmex->GetAircraft()->GetGearUnit(i)->SetReport(true);
+  for(i=0;i < fdmex->GetGroundReactions()->GetNumGearUnits();i++){
+    fdmex->GetGroundReactions()->GetGearUnit(i)->SetReport(true);
   }
   fdmex->GetOutput()->Enable();
   return !trim_failed;
@@ -436,9 +403,9 @@ bool FGTrim::DoTrim(void) {
 
 bool FGTrim::solve(void) {
 
-  float x1,x2,x3,f1,f2,f3,d,d0;
-  const float relax =0.9;
-  float eps=TrimAxes[current_axis]->GetSolverEps();
+  double x1,x2,x3,f1,f2,f3,d,d0;
+  const double relax =0.9;
+  double eps=TrimAxes[current_axis]->GetSolverEps();
 
   x1=x2=x3=0;
   d=1;
@@ -452,12 +419,11 @@ bool FGTrim::solve(void) {
       x1=xhi;f1=ahi;
       x3=xlo;f3=alo;
     }   */
-
     d0=fabs(x3-x1);
     //iterations
     //max_sub_iterations=TrimAxes[current_axis]->GetIterationLimit();
-    while (!TrimAxes[current_axis]->InTolerance() && (fabs(d) > eps) 
-              && (Nsub < max_sub_iterations)) {
+    while ( (TrimAxes[current_axis]->InTolerance() == false )
+             && (fabs(d) > eps) && (Nsub < max_sub_iterations)) {
       Nsub++;
       d=(x3-x1)/d0;
       x2=x1-d*d0*f1/(f3-f1);
@@ -518,12 +484,12 @@ bool FGTrim::solve(void) {
 */
 bool FGTrim::findInterval(void) {
   bool found=false;
-  float step;
-  float current_control=TrimAxes[current_axis]->GetControl();
-  float current_accel=TrimAxes[current_axis]->GetState();;
-  float xmin=TrimAxes[current_axis]->GetControlMin();
-  float xmax=TrimAxes[current_axis]->GetControlMax();
-  float lastxlo,lastxhi,lastalo,lastahi;
+  double step;
+  double current_control=TrimAxes[current_axis]->GetControl();
+  double current_accel=TrimAxes[current_axis]->GetState();;
+  double xmin=TrimAxes[current_axis]->GetControlMin();
+  double xmax=TrimAxes[current_axis]->GetControlMax();
+  double lastxlo,lastxhi,lastalo,lastahi;
   
   step=0.025*fabs(xmax);
   xlo=xhi=current_control;
@@ -588,8 +554,8 @@ bool FGTrim::findInterval(void) {
 
 bool FGTrim::checkLimits(void) {
   bool solutionExists;
-  float current_control=TrimAxes[current_axis]->GetControl();
-  float current_accel=TrimAxes[current_axis]->GetState();
+  double current_control=TrimAxes[current_axis]->GetControl();
+  double current_accel=TrimAxes[current_axis]->GetState();
   xlo=TrimAxes[current_axis]->GetControlMin();
   xhi=TrimAxes[current_axis]->GetControlMax();
 
@@ -622,5 +588,68 @@ bool FGTrim::checkLimits(void) {
   return solutionExists;
 }
 
+void FGTrim::setupPullup() {
+  float g,q,cgamma;
+  FGColumnVector3 vPQR;
+  g=fdmex->GetInertial()->gravity();
+  cgamma=cos(fgic->GetFlightPathAngleRadIC());
+  cout << "setPitchRateInPullup():  " << g << ", " << cgamma << ", "
+       << fgic->GetVtrueFpsIC() << endl;
+  q=g*(targetNlf-cgamma)/fgic->GetVtrueFpsIC();
+  cout << targetNlf << ", " << q << endl;
+  fdmex->GetRotation()->SetPQR(0,q,0);
+  cout << "setPitchRateInPullup() complete" << endl;
+  
+}  
+  
+void FGTrim::setupTurn(void){
+  double g,phi;
+  phi = fgic->GetRollAngleRadIC();
+  if( fabs(phi) > 0.001 && fabs(phi) < 1.56 ) {
+    targetNlf = 1 / cos(phi);
+    g = fdmex->GetInertial()->gravity(); 
+    psidot = g*tan(phi) / fgic->GetUBodyFpsIC();
+    cout << targetNlf << ", " << psidot << endl;
+  }
+   
+}  
+
+void FGTrim::updateRates(void){
+  if( mode == tTurn ) {
+    double phi = fgic->GetRollAngleRadIC();
+    double g = fdmex->GetInertial()->gravity(); 
+    double p,q,r,theta;
+    if(fabs(phi) > 0.001 && fabs(phi) < 1.56 ) {
+      theta=fgic->GetPitchAngleRadIC();
+      phi=fgic->GetRollAngleRadIC();
+      psidot = g*tan(phi) / fgic->GetUBodyFpsIC();
+      p=-psidot*sin(theta);
+      q=psidot*cos(theta)*sin(phi);
+      r=psidot*cos(theta)*cos(phi);
+    } else {
+      p=q=r=0;
+    }      
+    fdmex->GetRotation()->SetPQR(p,q,r);
+  } else if( mode == tPullup && fabs(targetNlf-1) > 0.01) {
+      float g,q,cgamma;
+      FGColumnVector3 vPQR;
+      g=fdmex->GetInertial()->gravity();
+      cgamma=cos(fgic->GetFlightPathAngleRadIC());
+      q=g*(targetNlf-cgamma)/fgic->GetVtrueFpsIC();
+      fdmex->GetRotation()->SetPQR(0,q,0);
+  }  
+}  
+
+void FGTrim::setDebug(void) {
+  if(debug_axis == tAll ||
+      TrimAxes[current_axis]->GetStateType() == debug_axis ) {
+    Debug=DebugLevel; 
+    return;
+  } else {
+    Debug=0;
+    return;
+  }
+}      
+    
 //YOU WERE WARNED, BUT YOU DID IT ANYWAY.