[Jsbsim-cvslogs] CVS: JSBSim/src/initialization FGInitialCondition.h, 1.44, 1.45 FGInitialCondition.cpp, 1.107, 1.108

[Bertrand <bco...@us...>]

Update of /cvsroot/jsbsim/JSBSim/src/initialization
In directory sfp-cvs-1.v30.ch3.sourceforge.com:/tmp/cvs-serv17298/src/initialization

Modified Files:
	FGInitialCondition.h FGInitialCondition.cpp 
Log Message:
Fixed the geodetic computations: now the altitude (ASL and AGL) is correctly set when a geodetic latitude is specified.

Index: FGInitialCondition.h
===================================================================
RCS file: /cvsroot/jsbsim/JSBSim/src/initialization/FGInitialCondition.h,v
retrieving revision 1.44
retrieving revision 1.45
diff -C2 -r1.44 -r1.45
*** FGInitialCondition.h	10 Jan 2016 16:35:28 -0000	1.44
--- FGInitialCondition.h	26 Jun 2016 20:33:04 -0000	1.45
***************
*** 678,681 ****
--- 678,682 ----
    FGMatrix33 Tw2b, Tb2w;
    double  alpha, beta;
+   double a, e2;
  
    speedset lastSpeedSet;
***************
*** 700,703 ****
--- 701,705 ----
    void calcAeroAngles(const FGColumnVector3& _vt_BODY);
    void calcThetaBeta(double alfa, const FGColumnVector3& _vt_NED);
+   bool LoadLatitude(Element* position_el);
    void Debug(int from);
  };

Index: FGInitialCondition.cpp
===================================================================
RCS file: /cvsroot/jsbsim/JSBSim/src/initialization/FGInitialCondition.cpp,v
retrieving revision 1.107
retrieving revision 1.108
diff -C2 -r1.107 -r1.108
*** FGInitialCondition.cpp	24 Jan 2016 18:18:38 -0000	1.107
--- FGInitialCondition.cpp	26 Jun 2016 20:33:04 -0000	1.108
***************
*** 126,131 ****
  {
    alpha=beta=0;
  
!   position.SetEllipse(fdmex->GetInertial()->GetSemimajor(), fdmex->GetInertial()->GetSemiminor());
  
    position.SetPositionGeodetic(0.0, 0.0, 0.0);
--- 126,135 ----
  {
    alpha=beta=0;
+   a = fdmex->GetInertial()->GetSemimajor();
+   double b = fdmex->GetInertial()->GetSemiminor();
+   double ec = b/a;
+   e2 = 1.0 - ec*ec;
  
!   position.SetEllipse(a, b);
  
    position.SetPositionGeodetic(0.0, 0.0, 0.0);
***************
*** 905,934 ****
  
  //******************************************************************************
  
! bool FGInitialCondition::Load_v1(Element* document)
  {
!   bool result = true;
! 
!   if (document->FindElement("longitude"))
!     SetLongitudeRadIC(document->FindElementValueAsNumberConvertTo("longitude", "RAD"));
!   if (document->FindElement("elevation"))
!     SetTerrainElevationFtIC(document->FindElementValueAsNumberConvertTo("elevation", "FT"));
! 
!   if (document->FindElement("altitude")) // This is feet above ground level
!     SetAltitudeAGLFtIC(document->FindElementValueAsNumberConvertTo("altitude", "FT"));
!   else if (document->FindElement("altitudeAGL")) // This is feet above ground level
!     SetAltitudeAGLFtIC(document->FindElementValueAsNumberConvertTo("altitudeAGL", "FT"));
!   else if (document->FindElement("altitudeMSL")) // This is feet above sea level
!     SetAltitudeASLFtIC(document->FindElementValueAsNumberConvertTo("altitudeMSL", "FT"));
  
-   double altitude = GetAltitudeASLFtIC();
-   double longitude = GetLongitudeRadIC();
- 
-   Element* latitude_el = document->FindElement("latitude");
    if (latitude_el) {
!     double latitude = document->FindElementValueAsNumberConvertTo("latitude", "RAD");
      if (fabs(latitude) > 0.5*M_PI) {
        string unit_type = latitude_el->GetAttributeValue("unit");
        if (unit_type.empty()) unit_type="RAD";
        cerr << latitude_el->ReadFrom() << "The latitude value "
             << latitude_el->GetDataAsNumber() << " " << unit_type
--- 909,933 ----
  
  //******************************************************************************
+ // Load the latitude from the XML file. The computations below assume that the
+ // terrain is a sphere and that the elevation is uniform all over the Earth.
+ // Would that assumption fail, the computation below would need to be adapted
+ // since the position radius would depend on the terrain elevation which depends
+ // itself on the latitude.
+ //
+ // This is an acceptable trade off because this routine is only used by
+ // standalone JSBSim which uses FGDefaultGroundCallback which assumes that the
+ // Earth is a sphere.
  
! bool FGInitialCondition::LoadLatitude(Element* position_el)
  {
!   Element* latitude_el = position_el->FindElement("latitude");
  
    if (latitude_el) {
!     double latitude = position_el->FindElementValueAsNumberConvertTo("latitude", "RAD");
! 
      if (fabs(latitude) > 0.5*M_PI) {
        string unit_type = latitude_el->GetAttributeValue("unit");
        if (unit_type.empty()) unit_type="RAD";
+ 
        cerr << latitude_el->ReadFrom() << "The latitude value "
             << latitude_el->GetDataAsNumber() << " " << unit_type
***************
*** 938,951 ****
        else
          cerr << "-PI/2 RAD; +PI/2 RAD]" << endl;
!       result = false;
      }
  
      string lat_type = latitude_el->GetAttributeValue("type");
!     if (lat_type == "geod" || lat_type == "geodetic")
!       position.SetPositionGeodetic(longitude, latitude, altitude); // Longitude and altitude will be set later on
      else
        position.SetLatitude(latitude);
    }
  
    FGColumnVector3 vOrient = orientation.GetEuler();
  
--- 937,982 ----
        else
          cerr << "-PI/2 RAD; +PI/2 RAD]" << endl;
! 
!       return false;
      }
  
      string lat_type = latitude_el->GetAttributeValue("type");
! 
!     if (lat_type == "geod" || lat_type == "geodetic") {
!       double R = position.GetRadius();
!       double slat = sin(latitude);
!       double RN = a / sqrt(1.0 - e2*slat*slat);
!       double p1 = e2*RN*slat*slat;
!       double p2 = e2*e2*RN*RN*slat*slat-R*R;
!       double h = p1 + sqrt(p1*p1-p2) - RN;
!       position.SetPositionGeodetic(position.GetLongitude(), latitude, h);
!     }
      else
        position.SetLatitude(latitude);
    }
  
+   return true;
+ }
+ 
+ //******************************************************************************
+ 
+ bool FGInitialCondition::Load_v1(Element* document)
+ {
+   bool result = true;
+ 
+   if (document->FindElement("longitude"))
+     SetLongitudeRadIC(document->FindElementValueAsNumberConvertTo("longitude", "RAD"));
+   if (document->FindElement("elevation"))
+     SetTerrainElevationFtIC(document->FindElementValueAsNumberConvertTo("elevation", "FT"));
+ 
+   if (document->FindElement("altitude")) // This is feet above ground level
+     SetAltitudeAGLFtIC(document->FindElementValueAsNumberConvertTo("altitude", "FT"));
+   else if (document->FindElement("altitudeAGL")) // This is feet above ground level
+     SetAltitudeAGLFtIC(document->FindElementValueAsNumberConvertTo("altitudeAGL", "FT"));
+   else if (document->FindElement("altitudeMSL")) // This is feet above sea level
+     SetAltitudeASLFtIC(document->FindElementValueAsNumberConvertTo("altitudeMSL", "FT"));
+ 
+   result = LoadLatitude(document);
+ 
    FGColumnVector3 vOrient = orientation.GetEuler();
  
***************
*** 996,1002 ****
      SetCrossWindKtsIC(document->FindElementValueAsNumberConvertTo("xwind", "KTS"));
    if (document->FindElement("targetNlf"))
-   {
      SetTargetNlfIC(document->FindElementValueAsNumber("targetNlf"));
-   }
    if (document->FindElement("trim"))
      needTrim = document->FindElementValueAsNumber("trim");
--- 1027,1031 ----
***************
*** 1007,1013 ****
    double radInv = 1.0 / position.GetRadius();
    FGColumnVector3 vOmegaLocal = FGColumnVector3(
!    radInv*vUVW_NED(eEast),
!   -radInv*vUVW_NED(eNorth),
!   -radInv*vUVW_NED(eEast)*position.GetTanLatitude() );
  
    vPQR_body = Tl2b * vOmegaLocal;
--- 1036,1042 ----
    double radInv = 1.0 / position.GetRadius();
    FGColumnVector3 vOmegaLocal = FGColumnVector3(
!                                                 radInv*vUVW_NED(eEast),
!                                                 -radInv*vUVW_NED(eNorth),
!                                                 -radInv*vUVW_NED(eEast)*position.GetTanLatitude() );
  
    vPQR_body = Tl2b * vOmegaLocal;
***************
*** 1036,1039 ****
--- 1065,1071 ----
    FGColumnVector3 vOmegaEarth = fdmex->GetInertial()->GetOmegaPlanet();
  
+   if (document->FindElement("elevation"))
+     fdmex->GetGroundCallback()->SetTerrainGeoCentRadius(document->FindElementValueAsNumberConvertTo("elevation", "FT")+position.GetSeaLevelRadius());
+ 
    // Initialize vehicle position
    //
***************
*** 1050,1054 ****
      } else if (frame == "ecef") {
        if (!position_el->FindElement("x") && !position_el->FindElement("y") && !position_el->FindElement("z")) {
-         Element* latitude_el = position_el->FindElement("latitude");
          if (position_el->FindElement("longitude"))
            position.SetLongitude(position_el->FindElementValueAsNumberConvertTo("longitude", "RAD"));
--- 1082,1085 ----
***************
*** 1065,1079 ****
          }
  
!         double altitude = position.GetAltitudeASL();
!         double longitude = position.GetLongitude();
! 
!         if (latitude_el) {
!           string lat_type = latitude_el->GetAttributeValue("type");
!           double latitude = position_el->FindElementValueAsNumberConvertTo("latitude", "RAD");
!           if (lat_type == "geod" || lat_type == "geodetic")
!             position.SetPositionGeodetic(longitude, latitude, altitude);
!           else
!             position.SetLatitude(latitude);
!         }
  
        } else {
--- 1096,1100 ----
          }
  
!         result = LoadLatitude(position_el);
  
        } else {
***************
*** 1089,1095 ****
    }
  
-   if (document->FindElement("elevation"))
-     fdmex->GetGroundCallback()->SetTerrainGeoCentRadius(document->FindElementValueAsNumberConvertTo("elevation", "FT")+position.GetSeaLevelRadius());
- 
    // End of position initialization
  
--- 1110,1113 ----