src/FDM/JSBSim/initialization/FGInitialCondition.cpp

   1 /*******************************************************************************
   2
   3  Header:       FGInitialCondition.cpp
   4  Author:       Tony Peden
   5  Date started: 7/1/99
   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 Lesser 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 Lesser General Public License for more
  17  details.
  18
  19  You should have received a copy of the GNU Lesser 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 Lesser 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/1/99   TP   Created
  30
  31
  32 FUNCTIONAL DESCRIPTION
  33 --------------------------------------------------------------------------------
  34
  35 The purpose of this class is to take a set of initial conditions and provide
  36 a kinematically consistent set of body axis velocity components, euler
  37 angles, and altitude.  This class does not attempt to trim the model i.e.
  38 the sim will most likely start in a very dynamic state (unless, of course,
  39 you have chosen your IC's wisely) even after setting it up with this class.
  40
  41 ********************************************************************************
  42 INCLUDES
  43 *******************************************************************************/
  44
  45 #include "FGInitialCondition.h"
  46 #include "FGFDMExec.h"
  47 #include "models/FGInertial.h"
  48 #include "models/FGAtmosphere.h"
  49 #include "models/FGAerodynamics.h"
  50 #include "models/FGPropagate.h"
  51 #include "input_output/FGPropertyManager.h"
  52 #include "input_output/FGXMLElement.h"
  53 #include "models/FGPropulsion.h"
  54 #include "input_output/FGXMLParse.h"
  55 #include "math/FGQuaternion.h"
  56 #include <iostream>
  57 #include <fstream>
  58 #include <cstdlib>
  59 #include "input_output/string_utilities.h"
  60
  61 using namespace std;
  62
  63 namespace JSBSim {
  64
  65 static const char *IdSrc = "$Id: FGInitialCondition.cpp,v 1.40 2010/06/02 04:05:13 jberndt Exp $";
  66 static const char *IdHdr = ID_INITIALCONDITION;
  67
  68 //******************************************************************************
  69
  70 FGInitialCondition::FGInitialCondition(FGFDMExec *FDMExec) : fdmex(FDMExec)
  71 {
  72   InitializeIC();
  73
  74   if(FDMExec != NULL ) {
  75     fdmex->GetPropagate()->SetAltitudeASL(altitudeASL);
  76     fdmex->GetAtmosphere()->Run();
  77     PropertyManager=fdmex->GetPropertyManager();
  78     Constructing = true;
  79     bind();
  80     Constructing = false;
  81   } else {
  82     cout << "FGInitialCondition: This class requires a pointer to a valid FGFDMExec object" << endl;
  83   }
  84
  85   Debug(0);
  86 }
  87
  88 //******************************************************************************
  89
  90 FGInitialCondition::~FGInitialCondition()
  91 {
  92   Debug(1);
  93 }
  94
  95 //******************************************************************************
  96
  97 void FGInitialCondition::ResetIC(double u0, double v0, double w0,
  98                                  double p0, double q0, double r0,
  99                                  double alpha0, double beta0,
 100                                  double phi0, double theta0, double psi0,
 101                                  double latRad0, double lonRad0, double altAGLFt0,
 102                                  double gamma0)
 103 {
 104   InitializeIC();
 105
 106   u = u0;  v = v0;  w = w0;
 107   p = p0;  q = q0;  r = r0;
 108   alpha = alpha0;  beta = beta0;
 109   phi = phi0;  theta = theta0;  psi = psi0;
 110   gamma = gamma0;
 111
 112   latitude = latRad0;
 113   longitude = lonRad0;
 114   SetAltitudeAGLFtIC(altAGLFt0);
 115
 116   cphi   = cos(phi);   sphi   = sin(phi);   // phi, rad
 117   ctheta = cos(theta); stheta = sin(theta); // theta, rad
 118   cpsi   = cos(psi);   spsi   = sin(psi);   // psi, rad
 119
 120   FGQuaternion Quat( phi, theta, psi );
 121   Quat.Normalize();
 122
 123   const FGMatrix33& _Tl2b  = Quat.GetT();     // local to body frame
 124   const FGMatrix33& _Tb2l  = Quat.GetTInv();  // body to local
 125
 126   FGColumnVector3 _vUVW_BODY(u,v,w);
 127   FGColumnVector3 _vUVW_NED = _Tb2l * _vUVW_BODY;
 128   FGColumnVector3 _vWIND_NED(wnorth,weast,wdown);
 129 //  FGColumnVector3 _vUVWAero = _Tl2b * ( _vUVW_NED + _vWIND_NED );
 130
 131   uw=_vWIND_NED(1); vw=_vWIND_NED(2); ww=_vWIND_NED(3);
 132
 133 }
 134
 135 //******************************************************************************
 136
 137 void FGInitialCondition::InitializeIC(void)
 138 {
 139   vt=vc=ve=vg=0;
 140   mach=0;
 141   alpha=beta=gamma=0;
 142   theta=phi=psi=0;
 143   altitudeASL=hdot=0;
 144   latitude=longitude=0;
 145   u=v=w=0;
 146   p=q=r=0;
 147   uw=vw=ww=0;
 148   vnorth=veast=vdown=0;
 149   wnorth=weast=wdown=0;
 150   whead=wcross=0;
 151   wdir=wmag=0;
 152   lastSpeedSet=setvt;
 153   lastWindSet=setwned;
 154   radius_to_vehicle = sea_level_radius = fdmex->GetInertial()->GetRefRadius();
 155   terrain_elevation = 0;
 156
 157   targetNlfIC = 1.0;
 158
 159   salpha=sbeta=stheta=sphi=spsi=sgamma=0;
 160   calpha=cbeta=ctheta=cphi=cpsi=cgamma=1;
 161 }
 162
 163 //******************************************************************************
 164
 165 void FGInitialCondition::WriteStateFile(int num)
 166 {
 167   if (Constructing) return;
 168
 169   string filename = fdmex->GetFullAircraftPath();
 170
 171   if (filename.empty())
 172     filename = "initfile.xml";
 173   else
 174     filename.append("/initfile.xml");
 175
 176   ofstream outfile(filename.c_str());
 177   FGPropagate* Propagate = fdmex->GetPropagate();
 178
 179   if (outfile.is_open()) {
 180     outfile << "<?xml version=\"1.0\"?>" << endl;
 181     outfile << "<initialize name=\"reset00\">" << endl;
 182     outfile << "  <ubody unit=\"FT/SEC\"> " << Propagate->GetUVW(eX) << " </ubody> " << endl;
 183     outfile << "  <vbody unit=\"FT/SEC\"> " << Propagate->GetUVW(eY) << " </vbody> " << endl;
 184     outfile << "  <wbody unit=\"FT/SEC\"> " << Propagate->GetUVW(eZ) << " </wbody> " << endl;
 185     outfile << "  <phi unit=\"DEG\"> " << Propagate->GetEuler(ePhi) << " </phi>" << endl;
 186     outfile << "  <theta unit=\"DEG\"> " << Propagate->GetEuler(eTht) << " </theta>" << endl;
 187     outfile << "  <psi unit=\"DEG\"> " << Propagate->GetEuler(ePsi) << " </psi>" << endl;
 188     outfile << "  <longitude unit=\"DEG\"> " << Propagate->GetLongitudeDeg() << " </longitude>" << endl;
 189     outfile << "  <latitude unit=\"DEG\"> " << Propagate->GetLatitudeDeg() << " </latitude>" << endl;
 190     outfile << "  <altitude unit=\"FT\"> " << Propagate->GetAltitudeASL() << " </altitude>" << endl;
 191     outfile << "</initialize>" << endl;
 192     outfile.close();
 193   } else {
 194     cerr << "Could not open and/or write the state to the initial conditions file: " << filename << endl;
 195   }
 196 }
 197
 198 //******************************************************************************
 199
 200 void FGInitialCondition::SetVcalibratedKtsIC(double tt) {
 201
 202   if(getMachFromVcas(&mach,tt*ktstofps)) {
 203     //cout << "Mach: " << mach << endl;
 204     lastSpeedSet=setvc;
 205     vc=tt*ktstofps;
 206     vt=mach*fdmex->GetAtmosphere()->GetSoundSpeed();
 207     ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 208     //cout << "Vt: " << vt*fpstokts << " Vc: " << vc*fpstokts << endl;
 209   }
 210   else {
 211     cout << "Failed to get Mach number for given Vc and altitude, Vc unchanged." << endl;
 212     cout << "Please mail the set of initial conditions used to apeden@earthlink.net" << endl;
 213   }
 214 }
 215
 216 //******************************************************************************
 217
 218 void FGInitialCondition::SetVequivalentKtsIC(double tt) {
 219   ve=tt*ktstofps;
 220   lastSpeedSet=setve;
 221   vt=ve*1/sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 222   mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
 223   vc=calcVcas(mach);
 224 }
 225
 226 //******************************************************************************
 227
 228 void FGInitialCondition::SetVgroundFpsIC(double tt) {
 229   double ua,va,wa;
 230   double vxz;
 231
 232   vg=tt;
 233   lastSpeedSet=setvg;
 234   vnorth = vg*cos(psi); veast = vg*sin(psi); vdown = 0;
 235   calcUVWfromNED();
 236   ua = u + uw; va = v + vw; wa = w + ww;
 237   vt = sqrt( ua*ua + va*va + wa*wa );
 238   alpha = beta = 0;
 239   vxz = sqrt( u*u + w*w );
 240   if( w != 0 ) alpha = atan2( w, u );
 241   if( vxz != 0 ) beta = atan2( v, vxz );
 242   mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
 243   vc=calcVcas(mach);
 244   ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 245 }
 246
 247 //******************************************************************************
 248
 249 void FGInitialCondition::SetVtrueFpsIC(double tt) {
 250   vt=tt;
 251   lastSpeedSet=setvt;
 252   mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
 253   vc=calcVcas(mach);
 254   ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 255 }
 256
 257 //******************************************************************************
 258
 259 void FGInitialCondition::SetMachIC(double tt) {
 260   mach=tt;
 261   lastSpeedSet=setmach;
 262   vt=mach*fdmex->GetAtmosphere()->GetSoundSpeed();
 263   vc=calcVcas(mach);
 264   ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 265   //cout << "Vt: " << vt*fpstokts << " Vc: " << vc*fpstokts << endl;
 266 }
 267
 268 //******************************************************************************
 269
 270 void FGInitialCondition::SetClimbRateFpmIC(double tt) {
 271   SetClimbRateFpsIC(tt/60.0);
 272 }
 273
 274 //******************************************************************************
 275
 276 void FGInitialCondition::SetClimbRateFpsIC(double tt) {
 277
 278   if(vt > 0.1) {
 279     hdot=tt;
 280     gamma=asin(hdot/vt);
 281     sgamma=sin(gamma); cgamma=cos(gamma);
 282   }
 283 }
 284
 285 //******************************************************************************
 286
 287 void FGInitialCondition::SetFlightPathAngleRadIC(double tt) {
 288   gamma=tt;
 289   sgamma=sin(gamma); cgamma=cos(gamma);
 290   getTheta();
 291   hdot=vt*sgamma;
 292 }
 293
 294 //******************************************************************************
 295
 296 void FGInitialCondition::SetAlphaRadIC(double tt) {
 297   alpha=tt;
 298   salpha=sin(alpha); calpha=cos(alpha);
 299   getTheta();
 300 }
 301
 302 //******************************************************************************
 303
 304 void FGInitialCondition::SetThetaRadIC(double tt) {
 305   theta=tt;
 306   stheta=sin(theta); ctheta=cos(theta);
 307   getAlpha();
 308 }
 309
 310 //******************************************************************************
 311
 312 void FGInitialCondition::SetBetaRadIC(double tt) {
 313   beta=tt;
 314   sbeta=sin(beta); cbeta=cos(beta);
 315   getTheta();
 316
 317 }
 318
 319 //******************************************************************************
 320
 321 void FGInitialCondition::SetPhiRadIC(double tt) {
 322   phi=tt;
 323   sphi=sin(phi); cphi=cos(phi);
 324   getTheta();
 325 }
 326
 327 //******************************************************************************
 328
 329 void FGInitialCondition::SetPsiRadIC(double tt) {
 330     psi=tt;
 331     spsi=sin(psi); cpsi=cos(psi);
 332     calcWindUVW();
 333 }
 334
 335 //******************************************************************************
 336
 337 void FGInitialCondition::SetUBodyFpsIC(double tt) {
 338   u=tt;
 339   vt=sqrt(u*u + v*v + w*w);
 340   lastSpeedSet=setuvw;
 341 }
 342
 343 //******************************************************************************
 344
 345 void FGInitialCondition::SetVBodyFpsIC(double tt) {
 346   v=tt;
 347   vt=sqrt(u*u + v*v + w*w);
 348   lastSpeedSet=setuvw;
 349 }
 350
 351 //******************************************************************************
 352
 353 void FGInitialCondition::SetWBodyFpsIC(double tt) {
 354   w=tt;
 355   vt=sqrt( u*u + v*v + w*w );
 356   lastSpeedSet=setuvw;
 357 }
 358
 359
 360 //******************************************************************************
 361
 362 void FGInitialCondition::SetVNorthFpsIC(double tt) {
 363   double ua,va,wa;
 364   double vxz;
 365   vnorth = tt;
 366   calcUVWfromNED();
 367   ua = u + uw; va = v + vw; wa = w + ww;
 368   vt = sqrt( ua*ua + va*va + wa*wa );
 369   alpha = beta = 0;
 370   vxz = sqrt( u*u + w*w );
 371   if( w != 0 ) alpha = atan2( w, u );
 372   if( vxz != 0 ) beta = atan2( v, vxz );
 373   mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
 374   vc=calcVcas(mach);
 375   ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 376   lastSpeedSet=setned;
 377 }
 378
 379 //******************************************************************************
 380
 381 void FGInitialCondition::SetVEastFpsIC(double tt) {
 382   double ua,va,wa;
 383   double vxz;
 384   veast = tt;
 385   calcUVWfromNED();
 386   ua = u + uw; va = v + vw; wa = w + ww;
 387   vt = sqrt( ua*ua + va*va + wa*wa );
 388   alpha = beta = 0;
 389   vxz = sqrt( u*u + w*w );
 390   if( w != 0 ) alpha = atan2( w, u );
 391   if( vxz != 0 ) beta = atan2( v, vxz );
 392   mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
 393   vc=calcVcas(mach);
 394   ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 395   lastSpeedSet=setned;
 396 }
 397
 398 //******************************************************************************
 399
 400 void FGInitialCondition::SetVDownFpsIC(double tt) {
 401   double ua,va,wa;
 402   double vxz;
 403   vdown = tt;
 404   calcUVWfromNED();
 405   ua = u + uw; va = v + vw; wa = w + ww;
 406   vt = sqrt( ua*ua + va*va + wa*wa );
 407   alpha = beta = 0;
 408   vxz = sqrt( u*u + w*w );
 409   if( w != 0 ) alpha = atan2( w, u );
 410   if( vxz != 0 ) beta = atan2( v, vxz );
 411   mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
 412   vc=calcVcas(mach);
 413   ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 414   SetClimbRateFpsIC(-1*vdown);
 415   lastSpeedSet=setned;
 416 }
 417
 418 //******************************************************************************
 419
 420 double FGInitialCondition::GetUBodyFpsIC(void) const {
 421     if (lastSpeedSet == setvg || lastSpeedSet == setned)
 422       return u;
 423     else
 424       return vt*calpha*cbeta - uw;
 425 }
 426
 427 //******************************************************************************
 428
 429 double FGInitialCondition::GetVBodyFpsIC(void) const {
 430     if (lastSpeedSet == setvg || lastSpeedSet == setned)
 431       return v;
 432     else {
 433       return vt*sbeta - vw;
 434     }
 435 }
 436
 437 //******************************************************************************
 438
 439 double FGInitialCondition::GetWBodyFpsIC(void) const {
 440     if (lastSpeedSet == setvg || lastSpeedSet == setned)
 441       return w;
 442     else
 443       return vt*salpha*cbeta -ww;
 444 }
 445
 446 //******************************************************************************
 447
 448 void FGInitialCondition::SetWindNEDFpsIC(double wN, double wE, double wD ) {
 449   wnorth = wN; weast = wE; wdown = wD;
 450   lastWindSet = setwned;
 451   calcWindUVW();
 452   if(lastSpeedSet == setvg)
 453     SetVgroundFpsIC(vg);
 454 }
 455
 456 //******************************************************************************
 457
 458 void FGInitialCondition::SetCrossWindKtsIC(double cross){
 459     wcross=cross*ktstofps;
 460     lastWindSet=setwhc;
 461     calcWindUVW();
 462     if(lastSpeedSet == setvg)
 463       SetVgroundFpsIC(vg);
 464
 465 }
 466
 467 //******************************************************************************
 468
 469 // positive from left
 470 void FGInitialCondition::SetHeadWindKtsIC(double head){
 471     whead=head*ktstofps;
 472     lastWindSet=setwhc;
 473     calcWindUVW();
 474     if(lastSpeedSet == setvg)
 475       SetVgroundFpsIC(vg);
 476
 477 }
 478
 479 //******************************************************************************
 480
 481 void FGInitialCondition::SetWindDownKtsIC(double wD) {
 482     wdown=wD;
 483     calcWindUVW();
 484     if(lastSpeedSet == setvg)
 485       SetVgroundFpsIC(vg);
 486 }
 487
 488 //******************************************************************************
 489
 490 void FGInitialCondition::SetWindMagKtsIC(double mag) {
 491   wmag=mag*ktstofps;
 492   lastWindSet=setwmd;
 493   calcWindUVW();
 494   if(lastSpeedSet == setvg)
 495       SetVgroundFpsIC(vg);
 496 }
 497
 498 //******************************************************************************
 499
 500 void FGInitialCondition::SetWindDirDegIC(double dir) {
 501   wdir=dir*degtorad;
 502   lastWindSet=setwmd;
 503   calcWindUVW();
 504   if(lastSpeedSet == setvg)
 505       SetVgroundFpsIC(vg);
 506 }
 507
 508
 509 //******************************************************************************
 510
 511 void FGInitialCondition::calcWindUVW(void) {
 512
 513     switch(lastWindSet) {
 514       case setwmd:
 515         wnorth=wmag*cos(wdir);
 516         weast=wmag*sin(wdir);
 517       break;
 518       case setwhc:
 519         wnorth=whead*cos(psi) + wcross*cos(psi+M_PI/2);
 520         weast=whead*sin(psi) + wcross*sin(psi+M_PI/2);
 521       break;
 522       case setwned:
 523       break;
 524     }
 525     uw=wnorth*ctheta*cpsi +
 526        weast*ctheta*spsi -
 527        wdown*stheta;
 528     vw=wnorth*( sphi*stheta*cpsi - cphi*spsi ) +
 529         weast*( sphi*stheta*spsi + cphi*cpsi ) +
 530        wdown*sphi*ctheta;
 531     ww=wnorth*(cphi*stheta*cpsi + sphi*spsi) +
 532        weast*(cphi*stheta*spsi - sphi*cpsi) +
 533        wdown*cphi*ctheta;
 534
 535
 536     /* cout << "FGInitialCondition::calcWindUVW: wnorth, weast, wdown "
 537          << wnorth << ", " << weast << ", " << wdown << endl;
 538     cout << "FGInitialCondition::calcWindUVW: theta, phi, psi "
 539           << theta << ", " << phi << ", " << psi << endl;
 540     cout << "FGInitialCondition::calcWindUVW: uw, vw, ww "
 541           << uw << ", " << vw << ", " << ww << endl; */
 542
 543 }
 544
 545 //******************************************************************************
 546
 547 void FGInitialCondition::SetAltitudeASLFtIC(double tt)
 548 {
 549   altitudeASL=tt;
 550   fdmex->GetPropagate()->SetAltitudeASL(altitudeASL);
 551   fdmex->GetAtmosphere()->Run();
 552   //lets try to make sure the user gets what they intended
 553
 554   switch(lastSpeedSet) {
 555   case setned:
 556   case setuvw:
 557   case setvt:
 558     SetVtrueKtsIC(vt*fpstokts);
 559     break;
 560   case setvc:
 561     SetVcalibratedKtsIC(vc*fpstokts);
 562     break;
 563   case setve:
 564     SetVequivalentKtsIC(ve*fpstokts);
 565     break;
 566   case setmach:
 567     SetMachIC(mach);
 568     break;
 569   case setvg:
 570     SetVgroundFpsIC(vg);
 571     break;
 572   }
 573 }
 574
 575 //******************************************************************************
 576
 577 void FGInitialCondition::SetAltitudeAGLFtIC(double tt)
 578 {
 579   SetAltitudeASLFtIC(terrain_elevation + tt);
 580 }
 581
 582 //******************************************************************************
 583
 584 void FGInitialCondition::SetSeaLevelRadiusFtIC(double tt)
 585 {
 586   sea_level_radius = tt;
 587 }
 588
 589 //******************************************************************************
 590
 591 void FGInitialCondition::SetTerrainElevationFtIC(double tt)
 592 {
 593   terrain_elevation=tt;
 594 }
 595
 596 //******************************************************************************
 597
 598 void FGInitialCondition::calcUVWfromNED(void)
 599 {
 600   u=vnorth*ctheta*cpsi +
 601      veast*ctheta*spsi -
 602      vdown*stheta;
 603   v=vnorth*( sphi*stheta*cpsi - cphi*spsi ) +
 604      veast*( sphi*stheta*spsi + cphi*cpsi ) +
 605      vdown*sphi*ctheta;
 606   w=vnorth*( cphi*stheta*cpsi + sphi*spsi ) +
 607      veast*( cphi*stheta*spsi - sphi*cpsi ) +
 608      vdown*cphi*ctheta;
 609 }
 610
 611 //******************************************************************************
 612
 613 bool FGInitialCondition::getMachFromVcas(double *Mach,double vcas) {
 614
 615   bool result=false;
 616   double guess=1.5;
 617   xlo=xhi=0;
 618   xmin=0;xmax=50;
 619   sfunc=&FGInitialCondition::calcVcas;
 620   if(findInterval(vcas,guess)) {
 621     if(solve(&mach,vcas))
 622       result=true;
 623   }
 624   return result;
 625 }
 626
 627 //******************************************************************************
 628
 629 bool FGInitialCondition::getAlpha(void) {
 630   bool result=false;
 631   double guess=theta-gamma;
 632
 633   if(vt < 0.01) return 0;
 634
 635   xlo=xhi=0;
 636   xmin=fdmex->GetAerodynamics()->GetAlphaCLMin();
 637   xmax=fdmex->GetAerodynamics()->GetAlphaCLMax();
 638   sfunc=&FGInitialCondition::GammaEqOfAlpha;
 639   if(findInterval(0,guess)){
 640     if(solve(&alpha,0)){
 641       result=true;
 642       salpha=sin(alpha);
 643       calpha=cos(alpha);
 644     }
 645   }
 646   calcWindUVW();
 647   return result;
 648 }
 649
 650 //******************************************************************************
 651
 652 bool FGInitialCondition::getTheta(void) {
 653   bool result=false;
 654   double guess=alpha+gamma;
 655
 656   if(vt < 0.01) return 0;
 657
 658   xlo=xhi=0;
 659   xmin=-89;xmax=89;
 660   sfunc=&FGInitialCondition::GammaEqOfTheta;
 661   if(findInterval(0,guess)){
 662     if(solve(&theta,0)){
 663       result=true;
 664       stheta=sin(theta);
 665       ctheta=cos(theta);
 666     }
 667   }
 668   calcWindUVW();
 669   return result;
 670 }
 671
 672 //******************************************************************************
 673
 674 double FGInitialCondition::GammaEqOfTheta(double Theta) {
 675   double a,b,c;
 676   double sTheta,cTheta;
 677
 678   //theta=Theta; stheta=sin(theta); ctheta=cos(theta);
 679   sTheta=sin(Theta); cTheta=cos(Theta);
 680   calcWindUVW();
 681   a=wdown + vt*calpha*cbeta + uw;
 682   b=vt*sphi*sbeta + vw*sphi;
 683   c=vt*cphi*salpha*cbeta + ww*cphi;
 684   return vt*sgamma - ( a*sTheta - (b+c)*cTheta);
 685 }
 686
 687 //******************************************************************************
 688
 689 double FGInitialCondition::GammaEqOfAlpha(double Alpha) {
 690   double a,b,c;
 691   double sAlpha,cAlpha;
 692   sAlpha=sin(Alpha); cAlpha=cos(Alpha);
 693   a=wdown + vt*cAlpha*cbeta + uw;
 694   b=vt*sphi*sbeta + vw*sphi;
 695   c=vt*cphi*sAlpha*cbeta + ww*cphi;
 696
 697   return vt*sgamma - ( a*stheta - (b+c)*ctheta );
 698 }
 699
 700 //******************************************************************************
 701
 702 double FGInitialCondition::calcVcas(double Mach) {
 703
 704   double p=fdmex->GetAtmosphere()->GetPressure();
 705   double psl=fdmex->GetAtmosphere()->GetPressureSL();
 706   double rhosl=fdmex->GetAtmosphere()->GetDensitySL();
 707   double pt,A,B,D,vcas;
 708
 709   if (Mach < 0) Mach=0;
 710   if (Mach < 1)    //calculate total pressure assuming isentropic flow
 711     pt=p*pow((1 + 0.2*Mach*Mach),3.5);
 712   else {
 713     // shock in front of pitot tube, we'll assume its normal and use
 714     // the Rayleigh Pitot Tube Formula, i.e. the ratio of total
 715     // pressure behind the shock to the static pressure in front of
 716     // the normal shock assumption should not be a bad one -- most supersonic
 717     // aircraft place the pitot probe out front so that it is the forward
 718     // most point on the aircraft.  The real shock would, of course, take
 719     // on something like the shape of a rounded-off cone but, here again,
 720     // the assumption should be good since the opening of the pitot probe
 721     // is very small and, therefore, the effects of the shock curvature
 722     // should be small as well. AFAIK, this approach is fairly well accepted
 723     // within the aerospace community
 724
 725     B = 5.76*Mach*Mach/(5.6*Mach*Mach - 0.8);
 726
 727     // The denominator above is zero for Mach ~ 0.38, for which
 728     // we'll never be here, so we're safe
 729
 730     D = (2.8*Mach*Mach-0.4)*0.4167;
 731     pt = p*pow(B,3.5)*D;
 732   }
 733
 734   A = pow(((pt-p)/psl+1),0.28571);
 735   vcas = sqrt(7*psl/rhosl*(A-1));
 736   //cout << "calcVcas: vcas= " << vcas*fpstokts << " mach= " << Mach << " pressure: " << pt << endl;
 737   return vcas;
 738 }
 739
 740 //******************************************************************************
 741
 742 bool FGInitialCondition::findInterval(double x,double guess) {
 743   //void find_interval(inter_params &ip,eqfunc f,double y,double constant, int &flag){
 744
 745   int i=0;
 746   bool found=false;
 747   double flo,fhi,fguess;
 748   double lo,hi,step;
 749   step=0.1;
 750   fguess=(this->*sfunc)(guess)-x;
 751   lo=hi=guess;
 752   do {
 753     step=2*step;
 754     lo-=step;
 755     hi+=step;
 756     if(lo < xmin) lo=xmin;
 757     if(hi > xmax) hi=xmax;
 758     i++;
 759     flo=(this->*sfunc)(lo)-x;
 760     fhi=(this->*sfunc)(hi)-x;
 761     if(flo*fhi <=0) {  //found interval with root
 762       found=true;
 763       if(flo*fguess <= 0) {  //narrow interval down a bit
 764         hi=lo+step;    //to pass solver interval that is as
 765         //small as possible
 766       }
 767       else if(fhi*fguess <= 0) {
 768         lo=hi-step;
 769       }
 770     }
 771     //cout << "findInterval: i=" << i << " Lo= " << lo << " Hi= " << hi << endl;
 772   }
 773   while((found == 0) && (i <= 100));
 774   xlo=lo;
 775   xhi=hi;
 776   return found;
 777 }
 778
 779 //******************************************************************************
 780
 781 bool FGInitialCondition::solve(double *y,double x)
 782 {
 783   double x1,x2,x3,f1,f2,f3,d,d0;
 784   double eps=1E-5;
 785   double const relax =0.9;
 786   int i;
 787   bool success=false;
 788
 789   //initializations
 790   d=1;
 791   x2 = 0;
 792   x1=xlo;x3=xhi;
 793   f1=(this->*sfunc)(x1)-x;
 794   f3=(this->*sfunc)(x3)-x;
 795   d0=fabs(x3-x1);
 796
 797   //iterations
 798   i=0;
 799   while ((fabs(d) > eps) && (i < 100)) {
 800     d=(x3-x1)/d0;
 801     x2 = x1-d*d0*f1/(f3-f1);
 802
 803     f2=(this->*sfunc)(x2)-x;
 804     //cout << "solve x1,x2,x3: " << x1 << "," << x2 << "," << x3 << endl;
 805     //cout << "                " << f1 << "," << f2 << "," << f3 << endl;
 806
 807     if(fabs(f2) <= 0.001) {
 808       x1=x3=x2;
 809     } else if(f1*f2 <= 0.0) {
 810       x3=x2;
 811       f3=f2;
 812       f1=relax*f1;
 813     } else if(f2*f3 <= 0) {
 814       x1=x2;
 815       f1=f2;
 816       f3=relax*f3;
 817     }
 818     //cout << i << endl;
 819     i++;
 820   }//end while
 821   if(i < 100) {
 822     success=true;
 823     *y=x2;
 824   }
 825
 826   //cout << "Success= " << success << " Vcas: " << vcas*fpstokts << " Mach: " << x2 << endl;
 827   return success;
 828 }
 829
 830 //******************************************************************************
 831
 832 double FGInitialCondition::GetWindDirDegIC(void) const {
 833   if(weast != 0.0)
 834     return atan2(weast,wnorth)*radtodeg;
 835   else if(wnorth > 0)
 836     return 0.0;
 837   else
 838     return 180.0;
 839 }
 840
 841 //******************************************************************************
 842
 843 bool FGInitialCondition::Load(string rstfile, bool useStoredPath)
 844 {
 845   string sep = "/";
 846   if( useStoredPath ) {
 847     init_file_name = fdmex->GetFullAircraftPath() + sep + rstfile + ".xml";
 848   } else {
 849     init_file_name = rstfile;
 850   }
 851
 852   document = LoadXMLDocument(init_file_name);
 853
 854   // Make sure that the document is valid
 855   if (!document) {
 856     cerr << "File: " << init_file_name << " could not be read." << endl;
 857     exit(-1);
 858   }
 859
 860   if (document->GetName() != string("initialize")) {
 861     cerr << "File: " << init_file_name << " is not a reset file." << endl;
 862     exit(-1);
 863   }
 864
 865   double version = document->GetAttributeValueAsNumber("version");
 866   if (version == HUGE_VAL) {
 867     return Load_v1(); // Default to the old version
 868   } else if (version >= 3.0) {
 869     cerr << "Only initialization file formats 1 and 2 are currently supported" << endl;
 870     exit (-1);
 871   } else if (version >= 2.0) {
 872     return Load_v2();
 873   } else if (version >= 1.0) {
 874     return Load_v1();
 875   }
 876
 877 }
 878
 879 //******************************************************************************
 880
 881 bool FGInitialCondition::Load_v1(void)
 882 {
 883   int n;
 884
 885   if (document->FindElement("latitude"))
 886     SetLatitudeDegIC(document->FindElementValueAsNumberConvertTo("latitude", "DEG"));
 887   if (document->FindElement("longitude"))
 888     SetLongitudeDegIC(document->FindElementValueAsNumberConvertTo("longitude", "DEG"));
 889   if (document->FindElement("elevation"))
 890     SetTerrainElevationFtIC(document->FindElementValueAsNumberConvertTo("elevation", "FT"));
 891
 892   if (document->FindElement("altitude")) // This is feet above ground level
 893     SetAltitudeAGLFtIC(document->FindElementValueAsNumberConvertTo("altitude", "FT"));
 894   else if (document->FindElement("altitudeAGL")) // This is feet above ground level
 895     SetAltitudeAGLFtIC(document->FindElementValueAsNumberConvertTo("altitudeAGL", "FT"));
 896   else if (document->FindElement("altitudeMSL")) // This is feet above sea level
 897     SetAltitudeASLFtIC(document->FindElementValueAsNumberConvertTo("altitudeMSL", "FT"));
 898
 899   if (document->FindElement("ubody"))
 900     SetUBodyFpsIC(document->FindElementValueAsNumberConvertTo("ubody", "FT/SEC"));
 901   if (document->FindElement("vbody"))
 902     SetVBodyFpsIC(document->FindElementValueAsNumberConvertTo("vbody", "FT/SEC"));
 903   if (document->FindElement("wbody"))
 904     SetWBodyFpsIC(document->FindElementValueAsNumberConvertTo("wbody", "FT/SEC"));
 905   if (document->FindElement("vnorth"))
 906     SetVNorthFpsIC(document->FindElementValueAsNumberConvertTo("vnorth", "FT/SEC"));
 907   if (document->FindElement("veast"))
 908     SetVEastFpsIC(document->FindElementValueAsNumberConvertTo("veast", "FT/SEC"));
 909   if (document->FindElement("vdown"))
 910     SetVDownFpsIC(document->FindElementValueAsNumberConvertTo("vdown", "FT/SEC"));
 911   if (document->FindElement("winddir"))
 912     SetWindDirDegIC(document->FindElementValueAsNumberConvertTo("winddir", "DEG"));
 913   if (document->FindElement("vwind"))
 914     SetWindMagKtsIC(document->FindElementValueAsNumberConvertTo("vwind", "KTS"));
 915   if (document->FindElement("hwind"))
 916     SetHeadWindKtsIC(document->FindElementValueAsNumberConvertTo("hwind", "KTS"));
 917   if (document->FindElement("xwind"))
 918     SetCrossWindKtsIC(document->FindElementValueAsNumberConvertTo("xwind", "KTS"));
 919   if (document->FindElement("vc"))
 920     SetVcalibratedKtsIC(document->FindElementValueAsNumberConvertTo("vc", "KTS"));
 921   if (document->FindElement("vt"))
 922     SetVtrueKtsIC(document->FindElementValueAsNumberConvertTo("vt", "KTS"));
 923   if (document->FindElement("mach"))
 924     SetMachIC(document->FindElementValueAsNumber("mach"));
 925   if (document->FindElement("phi"))
 926     SetPhiDegIC(document->FindElementValueAsNumberConvertTo("phi", "DEG"));
 927   if (document->FindElement("theta"))
 928     SetThetaDegIC(document->FindElementValueAsNumberConvertTo("theta", "DEG"));
 929   if (document->FindElement("psi"))
 930     SetPsiDegIC(document->FindElementValueAsNumberConvertTo("psi", "DEG"));
 931   if (document->FindElement("alpha"))
 932     SetAlphaDegIC(document->FindElementValueAsNumberConvertTo("alpha", "DEG"));
 933   if (document->FindElement("beta"))
 934     SetBetaDegIC(document->FindElementValueAsNumberConvertTo("beta", "DEG"));
 935   if (document->FindElement("gamma"))
 936     SetFlightPathAngleDegIC(document->FindElementValueAsNumberConvertTo("gamma", "DEG"));
 937   if (document->FindElement("roc"))
 938     SetClimbRateFpsIC(document->FindElementValueAsNumberConvertTo("roc", "FT/SEC"));
 939   if (document->FindElement("vground"))
 940     SetVgroundKtsIC(document->FindElementValueAsNumberConvertTo("vground", "KTS"));
 941   if (document->FindElement("targetNlf"))
 942   {
 943     SetTargetNlfIC(document->FindElementValueAsNumber("targetNlf"));
 944   }
 945
 946   // Check to see if any engines are specified to be initialized in a running state
 947   FGPropulsion* propulsion = fdmex->GetPropulsion();
 948   Element* running_elements = document->FindElement("running");
 949   while (running_elements) {
 950     n = int(running_elements->GetDataAsNumber());
 951     propulsion->InitRunning(n);
 952     running_elements = document->FindNextElement("running");
 953   }
 954
 955   fdmex->RunIC();
 956
 957   return true;
 958 }
 959
 960 //******************************************************************************
 961
 962 bool FGInitialCondition::Load_v2(void)
 963 {
 964   int n;
 965   FGColumnVector3 vLoc, vOrient;
 966   bool result = true;
 967   FGInertial* Inertial = fdmex->GetInertial();
 968   FGPropagate* Propagate = fdmex->GetPropagate();
 969
 970   if (document->FindElement("earth_position_angle")) {
 971     double epa = document->FindElementValueAsNumberConvertTo("earth_position_angle", "RAD");
 972     Inertial->SetEarthPositionAngle(epa);
 973   }
 974
 975   // Initialize vehicle position
 976   //
 977   // Allowable frames:
 978   // - ECI (Earth Centered Inertial)
 979   // - ECEF (Earth Centered, Earth Fixed)
 980
 981   Element* position = document->FindElement("position");
 982   if (position) {
 983     vLoc = position->FindElementTripletConvertTo("FT");
 984     string frame = position->GetAttributeValue("frame");
 985     frame = to_lower(frame);
 986     if (frame == "eci") {
 987       // Need to transform vLoc to ECEF for storage and use in FGLocation.
 988       vLoc = Propagate->GetTi2ec()*vLoc;
 989     } else if (frame == "ecef") {
 990       // Move vLoc query until after lat/lon/alt query to eliminate spurious warning msgs.
 991       if (vLoc.Magnitude() == 0.0) {
 992         Propagate->SetLatitudeDeg(position->FindElementValueAsNumberConvertTo("latitude", "DEG"));
 993         Propagate->SetLongitudeDeg(position->FindElementValueAsNumberConvertTo("longitude", "DEG"));
 994         if (position->FindElement("radius")) {
 995           radius_to_vehicle = position->FindElementValueAsNumberConvertTo("radius", "FT");
 996           SetAltitudeASLFtIC(radius_to_vehicle - sea_level_radius);
 997         } else if (position->FindElement("altitudeAGL")) {
 998           SetAltitudeAGLFtIC(position->FindElementValueAsNumberConvertTo("altitudeAGL", "FT"));
 999         } else if (position->FindElement("altitudeMSL")) {
1000           SetAltitudeASLFtIC(position->FindElementValueAsNumberConvertTo("altitudeMSL", "FT"));
1001         } else {
1002           cerr << endl << "  No altitude or radius initial condition is given." << endl;
1003           result = false;
1004         }
1005       }
1006     } else {
1007       cerr << endl << "  Neither ECI nor ECEF frame is specified for initial position." << endl;
1008       result = false;
1009     }
1010   } else {
1011     cerr << endl << "  Initial position not specified in this initialization file." << endl;
1012     result = false;
1013   }
1014
1015   // End of position initialization
1016
1017   // Initialize vehicle orientation
1018   // Allowable frames
1019   // - ECI (Earth Centered Inertial)
1020   // - ECEF (Earth Centered, Earth Fixed)
1021   // - Local
1022   //
1023   // Need to convert the provided orientation to an ECI orientation, using
1024   // the given orientation and knowledge of the Earth position angle.
1025   // This could be done using matrices (where in the subscript "b/a",
1026   // it is meant "b with respect to a", and where b=body frame,
1027   // i=inertial frame, and e=ecef frame) as:
1028   //
1029   // C_b/i =  C_b/e * C_e/i
1030   //
1031   // Using quaternions (note reverse ordering compared to matrix representation):
1032   //
1033   // Q_b/i = Q_e/i * Q_b/e
1034   //
1035   // Use the specific matrices as needed. The above example of course is for the whole
1036   // body to inertial orientation.
1037   // The new orientation angles can be extracted from the matrix or the quaternion.
1038   // ToDo: Do we need to deal with normalization of the quaternions here?
1039
1040   Element* orientation_el = document->FindElement("orientation");
1041   if (orientation_el) {
1042     string frame = orientation_el->GetAttributeValue("frame");
1043     frame = to_lower(frame);
1044     vOrient = orientation_el->FindElementTripletConvertTo("RAD");
1045     FGQuaternion QuatI2Body;
1046     if (frame == "eci") {
1047
1048       QuatI2Body = FGQuaternion(vOrient);
1049
1050     } else if (frame == "ecef") {
1051
1052       // In this case we are given the Euler angles representing the orientation of
1053       // the body with respect to the ECEF system, represented by the C_b/e Matrix.
1054       // We want the body orientation with respect to the inertial system:
1055       //
1056       // C_b/i =  C_b/e * C_e/i
1057       //
1058       // Using quaternions (note reverse ordering compared to matrix representation):
1059       //
1060       // Q_b/i = Q_e/i * Q_b/e
1061
1062       FGQuaternion QuatEC2Body(vOrient); // Store relationship of Body frame wrt ECEF frame, Q_b/e
1063       FGQuaternion QuatI2EC = Propagate->GetTi2ec(); // Get Q_e/i from matrix
1064       QuatI2Body = QuatI2EC * QuatEC2Body; // Q_b/i = Q_e/i * Q_b/e
1065
1066     } else if (frame == "local") {
1067
1068       // In this case, we are supplying the Euler angles for the vehicle with
1069       // respect to the local (NED frame), also called the navigation frame.
1070       // This is the most common way of initializing the orientation of
1071       // aircraft. The matrix representation is:
1072       //
1073       // C_b/i = C_b/n * C_n/e * C_e/i
1074       //
1075       // Or, using quaternions (note reverse ordering compared to matrix representation):
1076       //
1077       // Q_b/i = Q_e/i * Q_n/e * Q_b/n
1078
1079       FGQuaternion QuatLocal2Body = FGQuaternion(vOrient); // Store relationship of Body frame wrt local (NED) frame, Q_b/n
1080       FGQuaternion QuatEC2Local = Propagate->GetTec2l();   // Get Q_n/e from matrix
1081       FGQuaternion QuatI2EC = Propagate->GetTi2ec(); // Get Q_e/i from matrix
1082       QuatI2Body = QuatI2EC * QuatEC2Local * QuatLocal2Body; // Q_b/i = Q_e/i * Q_n/e * Q_b/n
1083
1084     } else {
1085
1086       cerr << endl << fgred << "  Orientation frame type: \"" << frame
1087            << "\" is not supported!" << reset << endl << endl;
1088       result = false;
1089
1090     }
1091
1092     Propagate->SetInertialOrientation(QuatI2Body);
1093   }
1094
1095   // Initialize vehicle velocity
1096   // Allowable frames
1097   // - ECI (Earth Centered Inertial)
1098   // - ECEF (Earth Centered, Earth Fixed)
1099   // - Local
1100   // - Body
1101   // The vehicle will be defaulted to (0,0,0) in the Body frame if nothing is provided.
1102
1103   Element* velocity_el = document->FindElement("velocity");
1104   FGColumnVector3 vInertialVelocity;
1105   FGColumnVector3 vInitVelocity = FGColumnVector3(0.0, 0.0, 0.0);
1106   if (velocity_el) {
1107
1108     string frame = velocity_el->GetAttributeValue("frame");
1109     frame = to_lower(frame);
1110     FGColumnVector3 vInitVelocity = velocity_el->FindElementTripletConvertTo("FT/SEC");
1111     FGColumnVector3 omega_cross_r = Inertial->omega() * Propagate->GetInertialPosition();
1112
1113     if (frame == "eci") {
1114       vInertialVelocity = vInitVelocity;
1115     } else if (frame == "ecef") {
1116       FGMatrix33 mTec2i = Propagate->GetTec2i(); // Get C_i/e
1117       vInertialVelocity = mTec2i * vInitVelocity + omega_cross_r;
1118     } else if (frame == "local") {
1119       FGMatrix33 mTl2i = Propagate->GetTl2i();
1120       vInertialVelocity = mTl2i * vInitVelocity + omega_cross_r;
1121     } else if (frame == "body") {
1122       FGMatrix33 mTb2i = Propagate->GetTb2i();
1123       vInertialVelocity = mTb2i * vInitVelocity + omega_cross_r;
1124     } else {
1125
1126       cerr << endl << fgred << "  Velocity frame type: \"" << frame
1127            << "\" is not supported!" << reset << endl << endl;
1128       result = false;
1129
1130     }
1131
1132   } else {
1133
1134     FGMatrix33 mTb2i = Propagate->GetTb2i();
1135     vInertialVelocity = mTb2i * vInitVelocity + (Inertial->omega() * Propagate->GetInertialPosition());
1136
1137   }
1138
1139   Propagate->SetInertialVelocity(vInertialVelocity);
1140
1141   // Allowable frames
1142   // - ECI (Earth Centered Inertial)
1143   // - ECEF (Earth Centered, Earth Fixed)
1144   // - Body
1145
1146   FGColumnVector3 vInertialRate;
1147   Element* attrate_el = document->FindElement("attitude_rate");
1148   if (attrate_el) {
1149
1150     string frame = attrate_el->GetAttributeValue("frame");
1151     frame = to_lower(frame);
1152     FGColumnVector3 vAttRate = attrate_el->FindElementTripletConvertTo("RAD/SEC");
1153
1154     if (frame == "eci") {
1155       vInertialRate = vAttRate;
1156     } else if (frame == "ecef") {
1157 //      vInertialRate = vAttRate + Inertial->omega();
1158     } else if (frame == "body") {
1159     //Todo: determine local frame rate
1160       FGMatrix33 mTb2l = Propagate->GetTb2l();
1161 //      vInertialRate = mTb2l*vAttRate + Inertial->omega();
1162     } else if (!frame.empty()) { // misspelling of frame
1163
1164       cerr << endl << fgred << "  Attitude rate frame type: \"" << frame
1165            << "\" is not supported!" << reset << endl << endl;
1166       result = false;
1167
1168     } else if (frame.empty()) {
1169
1170     }
1171
1172   } else { // Body frame attitude rate assumed 0 relative to local.
1173 /*
1174     //Todo: determine local frame rate
1175
1176     FGMatrix33 mTi2l = Propagate->GetTi2l();
1177     vVel = mTi2l * vInertialVelocity;
1178
1179     // Compute the local frame ECEF velocity
1180     vVel = Tb2l * VState.vUVW;
1181
1182     FGColumnVector3 vOmegaLocal = FGColumnVector3(
1183        radInv*vVel(eEast),
1184       -radInv*vVel(eNorth),
1185       -radInv*vVel(eEast)*VState.vLocation.GetTanLatitude() );
1186 */
1187   }
1188
1189   // Check to see if any engines are specified to be initialized in a running state
1190   FGPropulsion* propulsion = fdmex->GetPropulsion();
1191   Element* running_elements = document->FindElement("running");
1192   while (running_elements) {
1193     n = int(running_elements->GetDataAsNumber());
1194     propulsion->InitRunning(n);
1195     running_elements = document->FindNextElement("running");
1196   }
1197
1198   // fdmex->RunIC();
1199
1200   return result;
1201 }
1202
1203 //******************************************************************************
1204
1205 void FGInitialCondition::bind(void){
1206   PropertyManager->Tie("ic/vc-kts", this,
1207                        &FGInitialCondition::GetVcalibratedKtsIC,
1208                        &FGInitialCondition::SetVcalibratedKtsIC,
1209                        true);
1210   PropertyManager->Tie("ic/ve-kts", this,
1211                        &FGInitialCondition::GetVequivalentKtsIC,
1212                        &FGInitialCondition::SetVequivalentKtsIC,
1213                        true);
1214   PropertyManager->Tie("ic/vg-kts", this,
1215                        &FGInitialCondition::GetVgroundKtsIC,
1216                        &FGInitialCondition::SetVgroundKtsIC,
1217                        true);
1218   PropertyManager->Tie("ic/vt-kts", this,
1219                        &FGInitialCondition::GetVtrueKtsIC,
1220                        &FGInitialCondition::SetVtrueKtsIC,
1221                        true);
1222   PropertyManager->Tie("ic/mach", this,
1223                        &FGInitialCondition::GetMachIC,
1224                        &FGInitialCondition::SetMachIC,
1225                        true);
1226   PropertyManager->Tie("ic/roc-fpm", this,
1227                        &FGInitialCondition::GetClimbRateFpmIC,
1228                        &FGInitialCondition::SetClimbRateFpmIC,
1229                        true);
1230   PropertyManager->Tie("ic/gamma-deg", this,
1231                        &FGInitialCondition::GetFlightPathAngleDegIC,
1232                        &FGInitialCondition::SetFlightPathAngleDegIC,
1233                        true);
1234   PropertyManager->Tie("ic/alpha-deg", this,
1235                        &FGInitialCondition::GetAlphaDegIC,
1236                        &FGInitialCondition::SetAlphaDegIC,
1237                        true);
1238   PropertyManager->Tie("ic/beta-deg", this,
1239                        &FGInitialCondition::GetBetaDegIC,
1240                        &FGInitialCondition::SetBetaDegIC,
1241                        true);
1242   PropertyManager->Tie("ic/theta-deg", this,
1243                        &FGInitialCondition::GetThetaDegIC,
1244                        &FGInitialCondition::SetThetaDegIC,
1245                        true);
1246   PropertyManager->Tie("ic/phi-deg", this,
1247                        &FGInitialCondition::GetPhiDegIC,
1248                        &FGInitialCondition::SetPhiDegIC,
1249                        true);
1250   PropertyManager->Tie("ic/psi-true-deg", this,
1251                        &FGInitialCondition::GetPsiDegIC );
1252   PropertyManager->Tie("ic/lat-gc-deg", this,
1253                        &FGInitialCondition::GetLatitudeDegIC,
1254                        &FGInitialCondition::SetLatitudeDegIC,
1255                        true);
1256   PropertyManager->Tie("ic/long-gc-deg", this,
1257                        &FGInitialCondition::GetLongitudeDegIC,
1258                        &FGInitialCondition::SetLongitudeDegIC,
1259                        true);
1260   PropertyManager->Tie("ic/h-sl-ft", this,
1261                        &FGInitialCondition::GetAltitudeASLFtIC,
1262                        &FGInitialCondition::SetAltitudeASLFtIC,
1263                        true);
1264   PropertyManager->Tie("ic/h-agl-ft", this,
1265                        &FGInitialCondition::GetAltitudeAGLFtIC,
1266                        &FGInitialCondition::SetAltitudeAGLFtIC,
1267                        true);
1268   PropertyManager->Tie("ic/sea-level-radius-ft", this,
1269                        &FGInitialCondition::GetSeaLevelRadiusFtIC,
1270                        &FGInitialCondition::SetSeaLevelRadiusFtIC,
1271                        true);
1272   PropertyManager->Tie("ic/terrain-elevation-ft", this,
1273                        &FGInitialCondition::GetTerrainElevationFtIC,
1274                        &FGInitialCondition::SetTerrainElevationFtIC,
1275                        true);
1276   PropertyManager->Tie("ic/vg-fps", this,
1277                        &FGInitialCondition::GetVgroundFpsIC,
1278                        &FGInitialCondition::SetVgroundFpsIC,
1279                        true);
1280   PropertyManager->Tie("ic/vt-fps", this,
1281                        &FGInitialCondition::GetVtrueFpsIC,
1282                        &FGInitialCondition::SetVtrueFpsIC,
1283                        true);
1284   PropertyManager->Tie("ic/vw-bx-fps", this,
1285                        &FGInitialCondition::GetWindUFpsIC);
1286   PropertyManager->Tie("ic/vw-by-fps", this,
1287                        &FGInitialCondition::GetWindVFpsIC);
1288   PropertyManager->Tie("ic/vw-bz-fps", this,
1289                        &FGInitialCondition::GetWindWFpsIC);
1290   PropertyManager->Tie("ic/vw-north-fps", this,
1291                        &FGInitialCondition::GetWindNFpsIC);
1292   PropertyManager->Tie("ic/vw-east-fps", this,
1293                        &FGInitialCondition::GetWindEFpsIC);
1294   PropertyManager->Tie("ic/vw-down-fps", this,
1295                        &FGInitialCondition::GetWindDFpsIC);
1296   PropertyManager->Tie("ic/vw-mag-fps", this,
1297                        &FGInitialCondition::GetWindFpsIC);
1298   PropertyManager->Tie("ic/vw-dir-deg", this,
1299                        &FGInitialCondition::GetWindDirDegIC,
1300                        &FGInitialCondition::SetWindDirDegIC,
1301                        true);
1302
1303   PropertyManager->Tie("ic/roc-fps", this,
1304                        &FGInitialCondition::GetClimbRateFpsIC,
1305                        &FGInitialCondition::SetClimbRateFpsIC,
1306                        true);
1307   PropertyManager->Tie("ic/u-fps", this,
1308                        &FGInitialCondition::GetUBodyFpsIC,
1309                        &FGInitialCondition::SetUBodyFpsIC,
1310                        true);
1311   PropertyManager->Tie("ic/v-fps", this,
1312                        &FGInitialCondition::GetVBodyFpsIC,
1313                        &FGInitialCondition::SetVBodyFpsIC,
1314                        true);
1315   PropertyManager->Tie("ic/w-fps", this,
1316                        &FGInitialCondition::GetWBodyFpsIC,
1317                        &FGInitialCondition::SetWBodyFpsIC,
1318                        true);
1319   PropertyManager->Tie("ic/vn-fps", this,
1320                        &FGInitialCondition::GetVNorthFpsIC,
1321                        &FGInitialCondition::SetVNorthFpsIC,
1322                        true);
1323   PropertyManager->Tie("ic/ve-fps", this,
1324                        &FGInitialCondition::GetVEastFpsIC,
1325                        &FGInitialCondition::SetVEastFpsIC,
1326                        true);
1327   PropertyManager->Tie("ic/vd-fps", this,
1328                        &FGInitialCondition::GetVDownFpsIC,
1329                        &FGInitialCondition::SetVDownFpsIC,
1330                        true);
1331   PropertyManager->Tie("ic/gamma-rad", this,
1332                        &FGInitialCondition::GetFlightPathAngleRadIC,
1333                        &FGInitialCondition::SetFlightPathAngleRadIC,
1334                        true);
1335   PropertyManager->Tie("ic/alpha-rad", this,
1336                        &FGInitialCondition::GetAlphaRadIC,
1337                        &FGInitialCondition::SetAlphaRadIC,
1338                        true);
1339   PropertyManager->Tie("ic/theta-rad", this,
1340                        &FGInitialCondition::GetThetaRadIC,
1341                        &FGInitialCondition::SetThetaRadIC,
1342                        true);
1343   PropertyManager->Tie("ic/beta-rad", this,
1344                        &FGInitialCondition::GetBetaRadIC,
1345                        &FGInitialCondition::SetBetaRadIC,
1346                        true);
1347   PropertyManager->Tie("ic/phi-rad", this,
1348                        &FGInitialCondition::GetPhiRadIC,
1349                        &FGInitialCondition::SetPhiRadIC,
1350                        true);
1351   PropertyManager->Tie("ic/psi-true-rad", this,
1352                        &FGInitialCondition::GetPsiRadIC);
1353   PropertyManager->Tie("ic/lat-gc-rad", this,
1354                        &FGInitialCondition::GetLatitudeRadIC,
1355                        &FGInitialCondition::SetLatitudeRadIC,
1356                        true);
1357   PropertyManager->Tie("ic/long-gc-rad", this,
1358                        &FGInitialCondition::GetLongitudeRadIC,
1359                        &FGInitialCondition::SetLongitudeRadIC,
1360                        true);
1361   PropertyManager->Tie("ic/p-rad_sec", this,
1362                        &FGInitialCondition::GetPRadpsIC,
1363                        &FGInitialCondition::SetPRadpsIC,
1364                        true);
1365   PropertyManager->Tie("ic/q-rad_sec", this,
1366                        &FGInitialCondition::GetQRadpsIC,
1367                        &FGInitialCondition::SetQRadpsIC,
1368                        true);
1369   PropertyManager->Tie("ic/r-rad_sec", this,
1370                        &FGInitialCondition::GetRRadpsIC,
1371                        &FGInitialCondition::SetRRadpsIC,
1372                        true);
1373
1374   typedef int (FGInitialCondition::*iPMF)(void) const;
1375   PropertyManager->Tie("simulation/write-state-file",
1376                        this,
1377                        (iPMF)0,
1378                        &FGInitialCondition::WriteStateFile);
1379
1380 }
1381
1382 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1383 //    The bitmasked value choices are as follows:
1384 //    unset: In this case (the default) JSBSim would only print
1385 //       out the normally expected messages, essentially echoing
1386 //       the config files as they are read. If the environment
1387 //       variable is not set, debug_lvl is set to 1 internally
1388 //    0: This requests JSBSim not to output any messages
1389 //       whatsoever.
1390 //    1: This value explicity requests the normal JSBSim
1391 //       startup messages
1392 //    2: This value asks for a message to be printed out when
1393 //       a class is instantiated
1394 //    4: When this value is set, a message is displayed when a
1395 //       FGModel object executes its Run() method
1396 //    8: When this value is set, various runtime state variables
1397 //       are printed out periodically
1398 //    16: When set various parameters are sanity checked and
1399 //       a message is printed out when they go out of bounds
1400
1401 void FGInitialCondition::Debug(int from)
1402 {
1403   if (debug_lvl <= 0) return;
1404
1405   if (debug_lvl & 1) { // Standard console startup message output
1406   }
1407   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
1408     if (from == 0) cout << "Instantiated: FGInitialCondition" << endl;
1409     if (from == 1) cout << "Destroyed:    FGInitialCondition" << endl;
1410   }
1411   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
1412   }
1413   if (debug_lvl & 8 ) { // Runtime state variables
1414   }
1415   if (debug_lvl & 16) { // Sanity checking
1416   }
1417   if (debug_lvl & 64) {
1418     if (from == 0) { // Constructor
1419       cout << IdSrc << endl;
1420       cout << IdHdr << endl;
1421     }
1422   }
1423 }
1424 }