flightgear/src/FDM/JSBSim/models/FGGasCell.cpp

/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 Header:       FGGasCell.h
 Author:       Anders Gidenstam
 Date started: 01/21/2006

 ----- Copyright (C) 2006 - 2008  Anders Gidenstam (anders(at)gidenstam.org) --

 This program is free software; you can redistribute it and/or modify it under
 the terms of the GNU Lesser General Public License as published by the Free Software
 Foundation; either version 2 of the License, or (at your option) any later
 version.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
 details.

 You should have received a copy of the GNU Lesser General Public License along with
 this program; if not, write to the Free Software Foundation, Inc., 59 Temple
 Place - Suite 330, Boston, MA  02111-1307, USA.

 Further information about the GNU Lesser General Public License can also be found on
 the world wide web at http://www.gnu.org.

FUNCTIONAL DESCRIPTION
--------------------------------------------------------------------------------
See header file.

HISTORY
--------------------------------------------------------------------------------
01/21/2006  AG   Created

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
INCLUDES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

#include <FGFDMExec.h>
#include <models/FGAuxiliary.h>
#include <models/FGAtmosphere.h>
#include <models/FGInertial.h>
#include <models/FGMassBalance.h>
#include "FGGasCell.h"

using std::cerr;
using std::endl;
using std::cout;

namespace JSBSim {

static const char *IdSrc = "$Id$";
static const char *IdHdr = ID_GASCELL;

/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CLASS IMPLEMENTATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/* Constants. */
const double FGGasCell::R = 3.4071;              // [lbs ft/(mol Rankine)]
const double FGGasCell::M_air = 0.0019186;       // [slug/mol]
const double FGGasCell::M_hydrogen = 0.00013841; // [slug/mol]
const double FGGasCell::M_helium = 0.00027409;   // [slug/mol]

FGGasCell::FGGasCell(FGFDMExec* exec, Element* el, int num) : FGForce(exec)
{
  string token;
  Element* element;

  Auxiliary = exec->GetAuxiliary();
  Atmosphere = exec->GetAtmosphere();
  PropertyManager = exec->GetPropertyManager();
  Inertial = exec->GetInertial();
  MassBalance = exec->GetMassBalance();

  gasCellJ = FGMatrix33();
  gasCellM = FGColumnVector3();

  Buoyancy = MaxVolume = MaxOverpressure = Temperature = Pressure =
    Contents = Volume = dVolumeIdeal = 0.0;
  Xradius = Yradius = Zradius = Xwidth = Ywidth = Zwidth = 0.0;
  ValveCoefficient = ValveOpen = 0.0;
  CellNum = num;

  // NOTE: In the local system X points north, Y points east and Z points down.
  SetTransformType(FGForce::tLocalBody);

  type = el->GetAttributeValue("type");
  if      (type == "HYDROGEN") Type = ttHYDROGEN;
  else if (type == "HELIUM")   Type = ttHELIUM;
  else if (type == "AIR")      Type = ttAIR;
  else                         Type = ttUNKNOWN;

  element = el->FindElement("location");
  if (element) {
    vXYZ = element->FindElementTripletConvertTo("IN");
  } else {
    cerr << "Fatal Error: No location found for this gas cell." << endl;
    exit(-1);
  }
  if ((el->FindElement("x_radius") || el->FindElement("x_width")) &&
      (el->FindElement("y_radius") || el->FindElement("y_width")) &&
      (el->FindElement("z_radius") || el->FindElement("z_width"))) {

    if (el->FindElement("x_radius")) {
      Xradius = el->FindElementValueAsNumberConvertTo("x_radius", "FT");
    }
    if (el->FindElement("y_radius")) {
      Yradius = el->FindElementValueAsNumberConvertTo("y_radius", "FT");
    }
    if (el->FindElement("z_radius")) {
      Zradius = el->FindElementValueAsNumberConvertTo("z_radius", "FT");
    }

    if (el->FindElement("x_width")) {
      Xwidth = el->FindElementValueAsNumberConvertTo("x_width", "FT");
    }
    if (el->FindElement("y_width")) {
      Ywidth = el->FindElementValueAsNumberConvertTo("y_width", "FT");
    }
    if (el->FindElement("z_width")) {
      Zwidth = el->FindElementValueAsNumberConvertTo("z_width", "FT");
    }

    // The volume is a (potentially) extruded ellipsoid.
    // However, currently only a few combinations of radius and width are
    // fully supported.
    if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
        (Xwidth  == 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
      // Ellipsoid volume.
      MaxVolume = 4.0  * M_PI * Xradius * Yradius * Zradius / 3.0;
    } else if  ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
                (Xwidth  != 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
      // Cylindrical volume.
      MaxVolume = M_PI * Yradius * Zradius * Xwidth;
    } else {
      cerr << "Warning: Unsupported gas cell shape." << endl;
      MaxVolume = 
        (4.0  * M_PI * Xradius * Yradius * Zradius / 3.0 +
         M_PI * Yradius * Zradius * Xwidth +
         M_PI * Xradius * Zradius * Ywidth +
         M_PI * Xradius * Yradius * Zwidth +
         2.0  * Xradius * Ywidth * Zwidth +
         2.0  * Yradius * Xwidth * Zwidth +
         2.0  * Zradius * Xwidth * Ywidth +
         Xwidth * Ywidth * Zwidth);
    }
  } else {
    cerr << "Fatal Error: Gas cell shape must be given." << endl;
    exit(-1);
  }
  if (el->FindElement("max_overpressure")) {
    MaxOverpressure = el->FindElementValueAsNumberConvertTo("max_overpressure",
                                                            "LBS/FT2");
  }
  if (el->FindElement("fullness")) {
    const double Fullness = el->FindElementValueAsNumber("fullness");
    if (0 <= Fullness) { 
      Volume = Fullness * MaxVolume; 
    } else {
      cerr << "Warning: Invalid initial gas cell fullness value." << endl;
    }
  }  
  if (el->FindElement("valve_coefficient")) {
    ValveCoefficient =
      el->FindElementValueAsNumberConvertTo("valve_coefficient",
                                            "FT4*SEC/SLUG");
    ValveCoefficient = max(ValveCoefficient, 0.0);
  }

  // Initialize state
  SetLocation(vXYZ);

  if (Temperature == 0.0) {
    Temperature = Atmosphere->GetTemperature();
  }
  if (Pressure == 0.0) {
    Pressure = Atmosphere->GetPressure();
  }
  if (Volume != 0.0) {
    // Calculate initial gas content.
    Contents = Pressure * Volume / (R * Temperature);
    
    // Clip to max allowed value.
    const double IdealPressure = Contents * R * Temperature / MaxVolume;
    if (IdealPressure > Pressure + MaxOverpressure) {
      Contents = (Pressure + MaxOverpressure) * MaxVolume / (R * Temperature);
      Pressure = Pressure + MaxOverpressure;
    } else {
      Pressure = max(IdealPressure, Pressure);
    }
  } else {
    // Calculate initial gas content.
    Contents = Pressure * MaxVolume / (R * Temperature);
  }

  Volume = Contents * R * Temperature / Pressure;
  Mass = Contents * M_gas();

  // Bind relevant properties
  char property_name[80];
  snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/max_volume-ft3",
           CellNum);
  PropertyManager->Tie( property_name, &MaxVolume );
  PropertyManager->SetWritable( property_name, false );
  snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/temp-R",
           CellNum);
  PropertyManager->Tie( property_name, &Temperature );
  snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/pressure-psf",
           CellNum);
  PropertyManager->Tie( property_name, &Pressure );
  snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/volume-ft3",
           CellNum);
  PropertyManager->Tie( property_name, &Volume );
  snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/buoyancy-lbs",
           CellNum);
  PropertyManager->Tie( property_name, &Buoyancy );
  snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/contents-mol",
           CellNum);
  PropertyManager->Tie( property_name, &Contents );
  snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/valve_open",
           CellNum);
  PropertyManager->Tie( property_name, &ValveOpen );

  Debug(0);

  // Read heat transfer coefficients
  if (Element* heat = el->FindElement("heat")) {
    Element* function_element = heat->FindElement("function");
    while (function_element) {
      HeatTransferCoeff.push_back(new FGFunction(PropertyManager,
                                                 function_element));
      function_element = heat->FindNextElement("function");
    }
  }

  // Load ballonets if there are any
  if (Element* ballonet_element = el->FindElement("ballonet")) {
    while (ballonet_element) {
      Ballonet.push_back(new FGBallonet(exec,
                                        ballonet_element,
                                        Ballonet.size(),
                                        this));
      ballonet_element = el->FindNextElement("ballonet");
    }
  }

}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

FGGasCell::~FGGasCell()
{
  unsigned int i;

  for (i = 0; i < HeatTransferCoeff.size(); i++) delete HeatTransferCoeff[i];
  HeatTransferCoeff.clear();

  for (i = 0; i < Ballonet.size(); i++) delete Ballonet[i];
  Ballonet.clear();

  Debug(1);
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

void FGGasCell::Calculate(double dt)
{
  const double AirTemperature = Atmosphere->GetTemperature();  // [Rankine]
  const double AirPressure    = Atmosphere->GetPressure();     // [lbs/ft<66>]
  const double AirDensity     = Atmosphere->GetDensity();      // [slug/ft<66>]
  const double g = Inertial->gravity();                        // [lbs/slug]

  const double OldTemperature = Temperature;
  const double OldPressure    = Pressure;
  unsigned int i;
  const unsigned int no_ballonets = Ballonet.size();

  //-- Read ballonet state --
  // NOTE: This model might need a more proper integration technique. 
  double BallonetsVolume = 0.0;
  double BallonetsHeatFlow = 0.0;
  for (i = 0; i < no_ballonets; i++) {
    BallonetsVolume   += Ballonet[i]->GetVolume();
    BallonetsHeatFlow += Ballonet[i]->GetHeatFlow();
  }

  //-- Gas temperature --

  if (HeatTransferCoeff.size() > 0) {
    // The model is based on the ideal gas law.
    // However, it does look a bit fishy. Please verify.
    //   dT/dt = dU / (Cv n R)
    double dU = 0.0;
    for (i = 0; i < HeatTransferCoeff.size(); i++) {
      dU += HeatTransferCoeff[i]->GetValue();
    }
    // Don't include dt when accounting for adiabatic expansion/contraction.
    // The rate of adiabatic cooling looks about right: ~5.4 Rankine/1000ft. 
    if (Contents > 0) {
      Temperature +=
        (dU * dt - Pressure * dVolumeIdeal - BallonetsHeatFlow) /
        (Cv_gas() * Contents * R);
    } else {
      Temperature = AirTemperature;
    }
  } else {
    // No simulation of complex temperature changes.
    // Note: Making the gas cell behave adiabatically might be a better
    // option.
    Temperature = AirTemperature;
  }

  //-- Pressure --
  const double IdealPressure =
    Contents * R * Temperature / (MaxVolume - BallonetsVolume);
  if (IdealPressure > AirPressure + MaxOverpressure) {
    Pressure = AirPressure + MaxOverpressure;
  } else {
    Pressure = max(IdealPressure, AirPressure);
  }

  //-- Manual valving --

  // FIXME: Presently the effect of manual valving is computed using
  //        an ad hoc formula which might not be a good representation
  //        of reality.
  if ((ValveCoefficient > 0.0) && (ValveOpen > 0.0)) {
    // First compute the difference in pressure between the gas in the
    // cell and the air above it.
    // FixMe: CellHeight should depend on current volume.
    const double CellHeight = 2 * Zradius + Zwidth;                   // [ft]
    const double GasMass    = Contents * M_gas();                     // [slug]
    const double GasVolume  = Contents * R * Temperature / Pressure;  // [ft<66>]
    const double GasDensity = GasMass / GasVolume;
    const double DeltaPressure =
      Pressure + CellHeight * g * (AirDensity - GasDensity) - AirPressure;
    const double VolumeValved =
      ValveOpen * ValveCoefficient * DeltaPressure * dt;
    Contents =
      max(0.0, Contents - Pressure * VolumeValved / (R * Temperature));
  }

  //-- Update ballonets. --
  // Doing that here should give them the opportunity to react to the
  // new pressure.
  BallonetsVolume = 0.0;
  for (i = 0; i < no_ballonets; i++) {
    Ballonet[i]->Calculate(dt);
    BallonetsVolume += Ballonet[i]->GetVolume();
  }

  //-- Automatic safety valving. --
  if (Contents * R * Temperature / (MaxVolume - BallonetsVolume) >
      AirPressure + MaxOverpressure) {
    // Gas is automatically valved. Valving capacity is assumed to be infinite.
    // FIXME: This could/should be replaced by damage to the gas cell envelope.
    Contents =
      (AirPressure + MaxOverpressure) *
      (MaxVolume - BallonetsVolume) / (R * Temperature);
  }

  //-- Volume --
  Volume = Contents * R * Temperature / Pressure + BallonetsVolume;
  dVolumeIdeal =
    Contents * R * (Temperature / Pressure - OldTemperature / OldPressure);

  //-- Current buoyancy --
  // The buoyancy is computed using the atmospheres local density.
  Buoyancy = Volume * AirDensity * g;
  
  // Note: This is gross buoyancy. The weight of the gas itself and
  // any ballonets is not deducted here as the effects of the gas mass
  // is handled by FGMassBalance.
  vFn.InitMatrix(0.0, 0.0, - Buoyancy);

  // Compute the inertia of the gas cell.
  // Consider the gas cell as a shape of uniform density.
  // FIXME: If the cell isn't ellipsoid or cylindrical the inertia will
  //        be wrong.
  gasCellJ = FGMatrix33();
  const double mass = Contents * M_gas();
  double Ixx, Iyy, Izz;
  if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
      (Xwidth  == 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
    // Ellipsoid volume.
    Ixx = (1.0 / 5.0) * mass * (Yradius*Yradius + Zradius*Zradius);
    Iyy = (1.0 / 5.0) * mass * (Xradius*Xradius + Zradius*Zradius);
    Izz = (1.0 / 5.0) * mass * (Xradius*Xradius + Yradius*Yradius);     
  } else if  ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
              (Xwidth  != 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
    // Cylindrical volume (might not be valid with an elliptical cross-section).
    Ixx = (1.0 / 2.0) * mass * Yradius * Zradius;
    Iyy =
      (1.0 / 4.0) * mass * Yradius * Zradius +
      (1.0 / 12.0) * mass * Xwidth * Xwidth;
    Izz =
      (1.0 / 4.0) * mass * Yradius * Zradius +
      (1.0 / 12.0) * mass * Xwidth * Xwidth;
  } else {
    // Not supported. Revert to pointmass model.
    Ixx = Iyy = Izz = 0.0;
  }
  // The volume is symmetric, so Ixy = Ixz = Iyz = 0.
  gasCellJ(1,1) = Ixx;
  gasCellJ(2,2) = Iyy;
  gasCellJ(3,3) = Izz;
  Mass = mass;
  gasCellM.InitMatrix();
  gasCellM(eX) +=
    GetXYZ(eX) * Mass*slugtolb;
  gasCellM(eY) +=
    GetXYZ(eY) * Mass*slugtolb;
  gasCellM(eZ) +=
    GetXYZ(eZ) * Mass*slugtolb;

  if (no_ballonets > 0) {
    // Add the mass, moment and inertia of any ballonets.
    const FGColumnVector3 p = MassBalance->StructuralToBody( GetXYZ() );

    for (i = 0; i < no_ballonets; i++) {
      Mass += Ballonet[i]->GetMass();
       
      // Add ballonet moments.
      gasCellM(eX) +=
        Ballonet[i]->GetXYZ(eX) * Ballonet[i]->GetMass()*slugtolb;
      gasCellM(eY) +=
        Ballonet[i]->GetXYZ(eY) * Ballonet[i]->GetMass()*slugtolb;
      gasCellM(eZ) +=
        Ballonet[i]->GetXYZ(eZ) * Ballonet[i]->GetMass()*slugtolb;
      
      // Moments of inertia must be converted to the gas cell frame here.
      FGColumnVector3 v =
        MassBalance->StructuralToBody( Ballonet[i]->GetXYZ() ) - p;
      // Body basis is in FT. 
      const double mass = Ballonet[i]->GetMass();
      gasCellJ += Ballonet[i]->GetInertia() +
        FGMatrix33( 0,                - mass*v(1)*v(2), - mass*v(1)*v(3),
                    - mass*v(2)*v(1), 0,                - mass*v(2)*v(3),
                    - mass*v(3)*v(1), - mass*v(3)*v(2), 0 );
    }
  }
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//    The bitmasked value choices are as follows:
//    unset: In this case (the default) JSBSim would only print
//       out the normally expected messages, essentially echoing
//       the config files as they are read. If the environment
//       variable is not set, debug_lvl is set to 1 internally
//    0: This requests JSBSim not to output any messages
//       whatsoever.
//    1: This value explicity requests the normal JSBSim
//       startup messages
//    2: This value asks for a message to be printed out when
//       a class is instantiated
//    4: When this value is set, a message is displayed when a
//       FGModel object executes its Run() method
//    8: When this value is set, various runtime state variables
//       are printed out periodically
//    16: When set various parameters are sanity checked and
//       a message is printed out when they go out of bounds

void FGGasCell::Debug(int from)
{
  if (debug_lvl <= 0) return;

  if (debug_lvl & 1) { // Standard console startup message output
    if (from == 0) { // Constructor
      cout << "    Gas cell holds " << Contents << " mol " <<
        type << endl;
      cout << "      Cell location (X, Y, Z) (in.): " << vXYZ(eX) << ", " <<
        vXYZ(eY) << ", " << vXYZ(eZ) << endl;
      cout << "      Maximum volume: " << MaxVolume << " ft3" << endl;
      cout << "      Relief valve release pressure: " << MaxOverpressure << 
        " lbs/ft2" << endl;
      cout << "      Manual valve coefficient: " << ValveCoefficient << 
        " ft4*sec/slug" << endl;
      cout << "      Initial temperature: " << Temperature << " Rankine" <<
        endl;
      cout << "      Initial pressure: " << Pressure << " lbs/ft2" << endl;
      cout << "      Initial volume: " << Volume << " ft3" << endl;
      cout << "      Initial mass: " << GetMass() << " slug mass" << endl;
      cout << "      Initial weight: " << GetMass()*lbtoslug << " lbs force" <<
        endl;
      cout << "      Heat transfer: " << endl;
    }
  }
  if (debug_lvl & 2 ) { // Instantiation/Destruction notification
    if (from == 0) cout << "Instantiated: FGGasCell" << endl;
    if (from == 1) cout << "Destroyed:    FGGasCell" << endl;
  }
  if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
  }
  if (debug_lvl & 8 ) { // Runtime state variables    
      cout << "      " << type << " cell holds " << Contents << " mol " << endl;
      cout << "      Temperature: " << Temperature << " Rankine" << endl;
      cout << "      Pressure: " << Pressure << " lbs/ft2" << endl;
      cout << "      Volume: " << Volume << " ft3" << endl;
      cout << "      Mass: " << GetMass() << " slug mass" << endl;
      cout << "      Weight: " << GetMass()*lbtoslug << " lbs force" << endl;
  }
  if (debug_lvl & 16) { // Sanity checking
  }
  if (debug_lvl & 64) {
    if (from == 0) { // Constructor
      cout << IdSrc << endl;
      cout << IdHdr << endl;
    }
  }
}

/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CLASS IMPLEMENTATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
const double FGBallonet::R = 3.4071;              // [lbs ft/(mol Rankine)]
const double FGBallonet::M_air = 0.0019186;       // [slug/mol]
const double FGBallonet::Cv_air = 5.0/2.0;        // [??]

FGBallonet::FGBallonet(FGFDMExec* exec, Element* el, int num, FGGasCell* parent)
{
  string token;
  Element* element;

  Auxiliary = exec->GetAuxiliary();
  Atmosphere = exec->GetAtmosphere();
  PropertyManager = exec->GetPropertyManager();
  Inertial = exec->GetInertial();

  ballonetJ = FGMatrix33();

  MaxVolume = MaxOverpressure = Temperature = Pressure =
    Contents = Volume = dVolumeIdeal = dU = 0.0;
  Xradius = Yradius = Zradius = Xwidth = Ywidth = Zwidth = 0.0;
  ValveCoefficient = ValveOpen = 0.0;
  BlowerInput = NULL;
  CellNum = num;
  Parent = parent;

  // NOTE: In the local system X points north, Y points east and Z points down.
  element = el->FindElement("location");
  if (element) {
    vXYZ = element->FindElementTripletConvertTo("IN");
  } else {
    cerr << "Fatal Error: No location found for this ballonet." << endl;
    exit(-1);
  }
  if ((el->FindElement("x_radius") || el->FindElement("x_width")) &&
      (el->FindElement("y_radius") || el->FindElement("y_width")) &&
      (el->FindElement("z_radius") || el->FindElement("z_width"))) {

    if (el->FindElement("x_radius")) {
      Xradius = el->FindElementValueAsNumberConvertTo("x_radius", "FT");
    }
    if (el->FindElement("y_radius")) {
      Yradius = el->FindElementValueAsNumberConvertTo("y_radius", "FT");
    }
    if (el->FindElement("z_radius")) {
      Zradius = el->FindElementValueAsNumberConvertTo("z_radius", "FT");
    }

    if (el->FindElement("x_width")) {
      Xwidth = el->FindElementValueAsNumberConvertTo("x_width", "FT");
    }
    if (el->FindElement("y_width")) {
      Ywidth = el->FindElementValueAsNumberConvertTo("y_width", "FT");
    }
    if (el->FindElement("z_width")) {
      Zwidth = el->FindElementValueAsNumberConvertTo("z_width", "FT");
    }

    // The volume is a (potentially) extruded ellipsoid.
    // FIXME: However, currently only a few combinations of radius and
    //        width are fully supported.
    if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
        (Xwidth  == 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
      // Ellipsoid volume.
      MaxVolume = 4.0  * M_PI * Xradius * Yradius * Zradius / 3.0;
    } else if  ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
                (Xwidth  != 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
      // Cylindrical volume.
      MaxVolume = M_PI * Yradius * Zradius * Xwidth;
    } else {
      cerr << "Warning: Unsupported ballonet shape." << endl;
      MaxVolume = 
        (4.0  * M_PI * Xradius * Yradius * Zradius / 3.0 +
         M_PI * Yradius * Zradius * Xwidth +
         M_PI * Xradius * Zradius * Ywidth +
         M_PI * Xradius * Yradius * Zwidth +
         2.0  * Xradius * Ywidth * Zwidth +
         2.0  * Yradius * Xwidth * Zwidth +
         2.0  * Zradius * Xwidth * Ywidth +
         Xwidth * Ywidth * Zwidth);
    }
  } else {
    cerr << "Fatal Error: Ballonet shape must be given." << endl;
    exit(-1);
  }
  if (el->FindElement("max_overpressure")) {
    MaxOverpressure = el->FindElementValueAsNumberConvertTo("max_overpressure",
                                                            "LBS/FT2");
  }
  if (el->FindElement("fullness")) {
    const double Fullness = el->FindElementValueAsNumber("fullness");
    if (0 <= Fullness) { 
      Volume = Fullness * MaxVolume; 
    } else {
      cerr << "Warning: Invalid initial ballonet fullness value." << endl;
    }
  }  
  if (el->FindElement("valve_coefficient")) {
    ValveCoefficient =
      el->FindElementValueAsNumberConvertTo("valve_coefficient",
                                            "FT4*SEC/SLUG");
    ValveCoefficient = max(ValveCoefficient, 0.0);
  }

  // Initialize state
  if (Temperature == 0.0) {
    Temperature = Parent->GetTemperature();
  }
  if (Pressure == 0.0) {
    Pressure = Parent->GetPressure();
  }
  if (Volume != 0.0) {
    // Calculate initial air content.
    Contents = Pressure * Volume / (R * Temperature);
    
    // Clip to max allowed value.
    const double IdealPressure = Contents * R * Temperature / MaxVolume;
    if (IdealPressure > Pressure + MaxOverpressure) {
      Contents = (Pressure + MaxOverpressure) * MaxVolume / (R * Temperature);
      Pressure = Pressure + MaxOverpressure;
    } else {
      Pressure = max(IdealPressure, Pressure);
    }
  } else {
    // Calculate initial air content.
    Contents = Pressure * MaxVolume / (R * Temperature);
  }

  Volume = Contents * R * Temperature / Pressure;

  // Bind relevant properties
  char property_name[80];
  snprintf(property_name, 80,
           "buoyant_forces/gas-cell[%d]/ballonet[%d]/max_volume-ft3",
           Parent->GetIndex(),
           CellNum);
  PropertyManager->Tie( property_name, &MaxVolume );
  PropertyManager->SetWritable( property_name, false );
  snprintf(property_name, 80,
           "buoyant_forces/gas-cell[%d]/ballonet[%d]/temp-R",
           Parent->GetIndex(),
           CellNum);
  PropertyManager->Tie( property_name, &Temperature );
  snprintf(property_name, 80,
           "buoyant_forces/gas-cell[%d]/ballonet[%d]/pressure-psf",
           Parent->GetIndex(),
           CellNum);
  PropertyManager->Tie( property_name, &Pressure );
  snprintf(property_name, 80,
           "buoyant_forces/gas-cell[%d]/ballonet[%d]/volume-ft3",
           Parent->GetIndex(),
           CellNum);
  PropertyManager->Tie( property_name, &Volume );
  snprintf(property_name, 80,
           "buoyant_forces/gas-cell[%d]/ballonet[%d]/contents-mol",
           Parent->GetIndex(),
           CellNum);
  PropertyManager->Tie( property_name, &Contents );
  snprintf(property_name, 80,
           "buoyant_forces/gas-cell[%d]/ballonet[%d]/valve_open",
           Parent->GetIndex(),
           CellNum);
  PropertyManager->Tie( property_name, &ValveOpen );

  Debug(0);

  // Read heat transfer coefficients
  if (Element* heat = el->FindElement("heat")) {
    Element* function_element = heat->FindElement("function");
    while (function_element) {
      HeatTransferCoeff.push_back(new FGFunction(PropertyManager,
                                                 function_element));
      function_element = heat->FindNextElement("function");
    }
  }
  // Read blower input function
  if (Element* blower = el->FindElement("blower_input")) {
    Element* function_element = blower->FindElement("function");
    BlowerInput = new FGFunction(PropertyManager,
                                 function_element);
  }
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

FGBallonet::~FGBallonet()
{
  unsigned int i;

  for (i = 0; i < HeatTransferCoeff.size(); i++) delete HeatTransferCoeff[i];
  HeatTransferCoeff.clear();

  delete BlowerInput;
  BlowerInput = NULL;

  Debug(1);
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

void FGBallonet::Calculate(double dt)
{
  const double ParentPressure = Parent->GetPressure();         // [lbs/ft<66>]
  const double AirPressure    = Atmosphere->GetPressure();     // [lbs/ft<66>]

  const double OldTemperature = Temperature;
  const double OldPressure    = Pressure;
  unsigned int i;

  //-- Gas temperature --

  // The model is based on the ideal gas law.
  // However, it does look a bit fishy. Please verify.
  //   dT/dt = dU / (Cv n R)
  dU = 0.0;
  for (i = 0; i < HeatTransferCoeff.size(); i++) {
    dU += HeatTransferCoeff[i]->GetValue();
  }
  // dt is already accounted for in dVolumeIdeal.
  if (Contents > 0) {
    Temperature +=
      (dU * dt - Pressure * dVolumeIdeal) / (Cv_air * Contents * R);
  } else {
    Temperature = Parent->GetTemperature();
  }

  //-- Pressure --
  const double IdealPressure = Contents * R * Temperature / MaxVolume;
  // The pressure is at least that of the parent gas cell.
  Pressure = max(IdealPressure, ParentPressure);

  //-- Blower input --
  if (BlowerInput) {
    const double AddedVolume = BlowerInput->GetValue() * dt;
    if (AddedVolume > 0.0) {
      Contents += Pressure * AddedVolume / (R * Temperature);
    }
  }

  //-- Pressure relief and manual valving --
  // FIXME: Presently the effect of valving is computed using
  //        an ad hoc formula which might not be a good representation
  //        of reality.
  if ((ValveCoefficient > 0.0) &&
      ((ValveOpen > 0.0) || (Pressure > AirPressure + MaxOverpressure))) {
    const double DeltaPressure = Pressure - AirPressure;
    const double VolumeValved =
      ((Pressure > AirPressure + MaxOverpressure) ? 1.0 : ValveOpen) *
      ValveCoefficient * DeltaPressure * dt;
    // FIXME: Too small values of Contents sometimes leads to NaN.
    //        Currently the minimum is restricted to a safe value.
    Contents =
      max(1.0, Contents - Pressure * VolumeValved / (R * Temperature));
  }

  //-- Volume --
  Volume = Contents * R * Temperature / Pressure;
  dVolumeIdeal =
    Contents * R * (Temperature / Pressure - OldTemperature / OldPressure);

  // Compute the inertia of the ballonet.
  // Consider the ballonet as a shape of uniform density.
  // FIXME: If the ballonet isn't ellipsoid or cylindrical the inertia will
  //        be wrong.
  ballonetJ = FGMatrix33();
  const double mass = Contents * M_air;
  double Ixx, Iyy, Izz;
  if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
      (Xwidth  == 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
    // Ellipsoid volume.
    Ixx = (1.0 / 5.0) * mass * (Yradius*Yradius + Zradius*Zradius);
    Iyy = (1.0 / 5.0) * mass * (Xradius*Xradius + Zradius*Zradius);
    Izz = (1.0 / 5.0) * mass * (Xradius*Xradius + Yradius*Yradius);     
  } else if  ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
              (Xwidth  != 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
    // Cylindrical volume (might not be valid with an elliptical cross-section).
    Ixx = (1.0 / 2.0) * mass * Yradius * Zradius;
    Iyy =
      (1.0 / 4.0) * mass * Yradius * Zradius +
      (1.0 / 12.0) * mass * Xwidth * Xwidth;
    Izz =
      (1.0 / 4.0) * mass * Yradius * Zradius +
      (1.0 / 12.0) * mass * Xwidth * Xwidth;
  } else {
    // Not supported. Revert to pointmass model.
    Ixx = Iyy = Izz = 0.0;
  }
  // The volume is symmetric, so Ixy = Ixz = Iyz = 0.
  ballonetJ(1,1) = Ixx;
  ballonetJ(2,2) = Iyy;
  ballonetJ(3,3) = Izz;
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//    The bitmasked value choices are as follows:
//    unset: In this case (the default) JSBSim would only print
//       out the normally expected messages, essentially echoing
//       the config files as they are read. If the environment
//       variable is not set, debug_lvl is set to 1 internally
//    0: This requests JSBSim not to output any messages
//       whatsoever.
//    1: This value explicity requests the normal JSBSim
//       startup messages
//    2: This value asks for a message to be printed out when
//       a class is instantiated
//    4: When this value is set, a message is displayed when a
//       FGModel object executes its Run() method
//    8: When this value is set, various runtime state variables
//       are printed out periodically
//    16: When set various parameters are sanity checked and
//       a message is printed out when they go out of bounds

void FGBallonet::Debug(int from)
{
  if (debug_lvl <= 0) return;

  if (debug_lvl & 1) { // Standard console startup message output
    if (from == 0) { // Constructor
      cout << "      Ballonet holds " << Contents << " mol air" << endl;
      cout << "        Location (X, Y, Z) (in.): " << vXYZ(eX) << ", " <<
        vXYZ(eY) << ", " << vXYZ(eZ) << endl;
      cout << "        Maximum volume: " << MaxVolume << " ft3" << endl;
      cout << "        Relief valve release pressure: " << MaxOverpressure << 
        " lbs/ft2" << endl;
      cout << "        Relief valve coefficient: " << ValveCoefficient << 
        " ft4*sec/slug" << endl;
      cout << "        Initial temperature: " << Temperature << " Rankine" <<
        endl;
      cout << "        Initial pressure: " << Pressure << " lbs/ft2" << endl;
      cout << "        Initial volume: " << Volume << " ft3" << endl;
      cout << "        Initial mass: " << GetMass() << " slug mass" << endl;
      cout << "        Initial weight: " << GetMass()*lbtoslug <<
        " lbs force" << endl;
      cout << "        Heat transfer: " << endl;
    }
  }
  if (debug_lvl & 2 ) { // Instantiation/Destruction notification
    if (from == 0) cout << "Instantiated: FGBallonet" << endl;
    if (from == 1) cout << "Destroyed:    FGBallonet" << endl;
  }
  if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
  }
  if (debug_lvl & 8 ) { // Runtime state variables    
      cout << "        Ballonet holds " << Contents <<
        " mol air" << endl;
      cout << "        Temperature: " << Temperature << " Rankine" << endl;
      cout << "        Pressure: " << Pressure << " lbs/ft2" << endl;
      cout << "        Volume: " << Volume << " ft3" << endl;
      cout << "        Mass: " << GetMass() << " slug mass" << endl;
      cout << "        Weight: " << GetMass()*lbtoslug << " lbs force" << endl;
  }
  if (debug_lvl & 16) { // Sanity checking
  }
  if (debug_lvl & 64) {
    if (from == 0) { // Constructor
      cout << IdSrc << endl;
      cout << IdHdr << endl;
    }
  }
}
}