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