Better fix for a compilation problem with MSVC 2012

[flightgear.git] / src / FDM / JSBSim / models / propulsion / FGPropeller.cpp

diff --git a/src/FDM/JSBSim/models/propulsion/FGPropeller.cpp b/src/FDM/JSBSim/models/propulsion/FGPropeller.cpp
index 7837d39ff3f32c0f23f469d15ac33f9f64893582..52635b7ef6b13447b128e62963a02d525d230fc9 100644 (file)
--- a/src/FDM/JSBSim/models/propulsion/FGPropeller.cpp
+++ b/src/FDM/JSBSim/models/propulsion/FGPropeller.cpp
@@ -45,7 +45,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGPropeller.cpp,v 1.41 2011/11/17 21:07:30 jentron Exp $";
+static const char *IdSrc = "$Id: FGPropeller.cpp,v 1.44 2012/04/29 13:10:46 bcoconni Exp $";
 static const char *IdHdr = ID_PROPELLER;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -194,52 +194,77 @@ FGPropeller::~FGPropeller()
 
 double FGPropeller::Calculate(double EnginePower)
 {
+  FGColumnVector3 localAeroVel = Transform().Transposed() * in.AeroUVW;
   double omega, PowerAvailable;
 
-  double Vel = in.AeroUVW(eU);
+  double Vel = localAeroVel(eU);
   double rho = in.Density;
   double RPS = RPM/60.0;
 
   // Calculate helical tip Mach
   double Area = 0.25*Diameter*Diameter*M_PI;
   double Vtip = RPS * Diameter * M_PI;
-  HelicalTipMach = sqrt(Vtip*Vtip + Vel*Vel) / in.Soundspeed; 
+  HelicalTipMach = sqrt(Vtip*Vtip + Vel*Vel) / in.Soundspeed;
 
   PowerAvailable = EnginePower - GetPowerRequired();
 
   if (RPS > 0.0) J = Vel / (Diameter * RPS); // Calculate J normally
-  else           J = Vel / Diameter;      
+  else           J = Vel / Diameter;
 
   if (MaxPitch == MinPitch) {    // Fixed pitch prop
-         ThrustCoeff = cThrust->GetValue(J);
+    ThrustCoeff = cThrust->GetValue(J);
   } else {                       // Variable pitch prop
-         ThrustCoeff = cThrust->GetValue(J, Pitch);
+    ThrustCoeff = cThrust->GetValue(J, Pitch);
   }
- 
+
   // Apply optional scaling factor to Ct (default value = 1)
   ThrustCoeff *= CtFactor;
 
   // Apply optional Mach effects from CT_MACH table
   if (CtMach) ThrustCoeff *= CtMach->GetValue(HelicalTipMach);
 
-  if (P_Factor > 0.0001) {
-//    alpha = sin(fdmex->GetAuxiliary()->Getalpha() + FGThruster::GetPitch());
-//    beta  = sin(fdmex->GetAuxiliary()->Getbeta() + FGThruster::GetYaw());
-//    SetActingLocationY( GetLocationY() + P_Factor*alpha*Sense);
-//    SetActingLocationZ( GetLocationZ() + P_Factor*beta*Sense);
-    SetActingLocationY( GetLocationY() + P_Factor*in.Alpha*Sense);
-    SetActingLocationZ( GetLocationZ() + P_Factor*in.Beta*Sense);
-  }
-
   Thrust = ThrustCoeff*RPS*RPS*D4*rho;
 
-  // From B. W. McCormick, "Aerodynamics, Aeronautics, and Flight Mechanics"
-  // first edition, eqn. 6.15 (propeller analysis chapter).
-  Vinduced = 0.5 * (-Vel + sqrt(Vel*Vel + 2.0*Thrust/(rho*Area)));
+  // Induced velocity in the propeller disk area. This formula is obtained
+  // from momentum theory - see B. W. McCormick, "Aerodynamics, Aeronautics,
+  // and Flight Mechanics" 1st edition, eqn. 6.15 (propeller analysis chapter).
+  // Since Thrust and Vel can both be negative we need to adjust this formula
+  // To handle sign (direction) separately from magnitude.
+  double Vel2sum = Vel*abs(Vel) + 2.0*Thrust/(rho*Area);
+  
+  if( Vel2sum > 0.0)
+    Vinduced = 0.5 * (-Vel + sqrt(Vel2sum));
+  else
+    Vinduced = 0.5 * (-Vel - sqrt(-Vel2sum));
+
+  // We need to drop the case where the downstream velocity is opposite in
+  // direction to the aircraft velocity. For example, in such a case, the
+  // direction of the airflow on the tail would be opposite to the airflow on
+  // the wing tips. When such complicated airflows occur, the momentum theory
+  // breaks down and the formulas above are no longer applicable
+  // (see H. Glauert, "The Elements of Airfoil and Airscrew Theory",
+  // 2nd edition, §16.3, pp. 219-221)
+
+  if ((Vel+2.0*Vinduced)*Vel < 0.0)
+    Vinduced = 0.0; // We cannot calculate the induced velocity so let's assume it is zero.
+    
+  // P-factor is simulated by a shift of the acting location of the thrust.
+  // The shift is a multiple of the angle between the propeller shaft axis
+  // and the relative wind that goes through the propeller disk.
+  if (P_Factor > 0.0001) {
+    double tangentialVel = localAeroVel.Magnitude(eV, eW);
+
+    if (tangentialVel > 0.0001) {
+      double angle = atan2(tangentialVel, localAeroVel(eU));
+      double factor = Sense * P_Factor * angle / tangentialVel;
+      SetActingLocationY( GetLocationY() + factor * localAeroVel(eW));
+      SetActingLocationZ( GetLocationZ() + factor * localAeroVel(eV));
+    }
+  }
 
   omega = RPS*2.0*M_PI;
 
-  vFn(1) = Thrust;
+  vFn(eX) = Thrust;
 
   // The Ixx value and rotation speed given below are for rotation about the
   // natural axis of the engine. The transform takes place in the base class
@@ -273,7 +298,7 @@ double FGPropeller::GetPowerRequired(void)
   double RPS = RPM / 60.0;
 
   if (RPS != 0.0) J = Vel / (Diameter * RPS);
-  else            J = Vel / Diameter; 
+  else            J = Vel / Diameter;
 
   if (MaxPitch == MinPitch) {   // Fixed pitch prop
     cPReq = cPower->GetValue(J);
@@ -320,7 +345,7 @@ double FGPropeller::GetPowerRequired(void)
     } else { // Manual Pitch Mode, pitch is controlled externally
 
     }
-  
+
     cPReq = cPower->GetValue(J, Pitch);
   }