src/FDM/JSBSim/models/FGMassBalance.cpp

   1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   2
   3  Module:       FGMassBalance.cpp
   4  Author:       Jon S. Berndt
   5  Date started: 09/12/2000
   6  Purpose:      This module models weight and balance
   7
   8  ------------- Copyright (C) 2000  Jon S. Berndt (jon@jsbsim.org) --------------
   9
  10  This program is free software; you can redistribute it and/or modify it under
  11  the terms of the GNU Lesser General Public License as published by the Free Software
  12  Foundation; either version 2 of the License, or (at your option) any later
  13  version.
  14
  15  This program is distributed in the hope that it will be useful, but WITHOUT
  16  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  17  FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
  18  details.
  19
  20  You should have received a copy of the GNU Lesser General Public License along with
  21  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  22  Place - Suite 330, Boston, MA  02111-1307, USA.
  23
  24  Further information about the GNU Lesser General Public License can also be found on
  25  the world wide web at http://www.gnu.org.
  26
  27 FUNCTIONAL DESCRIPTION
  28 --------------------------------------------------------------------------------
  29
  30 This class models the change in weight and balance of the aircraft due to fuel
  31 burnoff, etc.
  32
  33 HISTORY
  34 --------------------------------------------------------------------------------
  35 09/12/2000  JSB  Created
  36
  37 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  38 INCLUDES
  39 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  40
  41 #include <iostream>
  42 #include <iomanip>
  43 #include <cstdlib>
  44 #include "FGMassBalance.h"
  45 #include "FGFDMExec.h"
  46 #include "input_output/FGPropertyManager.h"
  47
  48 using namespace std;
  49
  50 namespace JSBSim {
  51
  52 static const char *IdSrc = "$Id: FGMassBalance.cpp,v 1.39 2011/11/09 21:58:26 bcoconni Exp $";
  53 static const char *IdHdr = ID_MASSBALANCE;
  54
  55 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  56 CLASS IMPLEMENTATION
  57 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  58
  59
  60 FGMassBalance::FGMassBalance(FGFDMExec* fdmex) : FGModel(fdmex)
  61 {
  62   Name = "FGMassBalance";
  63   Weight = EmptyWeight = Mass = 0.0;
  64
  65   vbaseXYZcg.InitMatrix(0.0);
  66   vXYZcg.InitMatrix(0.0);
  67   vLastXYZcg.InitMatrix(0.0);
  68   vDeltaXYZcg.InitMatrix(0.0);
  69   baseJ.InitMatrix();
  70   mJ.InitMatrix();
  71   mJinv.InitMatrix();
  72   pmJ.InitMatrix();
  73
  74   bind();
  75
  76   Debug(0);
  77 }
  78
  79 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  80
  81 FGMassBalance::~FGMassBalance()
  82 {
  83   for (unsigned int i=0; i<PointMasses.size(); i++) delete PointMasses[i];
  84   PointMasses.clear();
  85
  86   Debug(1);
  87 }
  88
  89 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  90
  91 bool FGMassBalance::InitModel(void)
  92 {
  93   vLastXYZcg.InitMatrix(0.0);
  94   vDeltaXYZcg.InitMatrix(0.0);
  95
  96   return true;
  97 }
  98
  99 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 100
 101 bool FGMassBalance::Load(Element* el)
 102 {
 103   Element *element;
 104   string element_name = "";
 105   double bixx, biyy, bizz, bixy, bixz, biyz;
 106
 107   FGModel::Load(el); // Perform base class Load.
 108
 109   bixx = biyy = bizz = bixy = bixz = biyz = 0.0;
 110   if (el->FindElement("ixx"))
 111     bixx = el->FindElementValueAsNumberConvertTo("ixx", "SLUG*FT2");
 112   if (el->FindElement("iyy"))
 113     biyy = el->FindElementValueAsNumberConvertTo("iyy", "SLUG*FT2");
 114   if (el->FindElement("izz"))
 115     bizz = el->FindElementValueAsNumberConvertTo("izz", "SLUG*FT2");
 116   if (el->FindElement("ixy"))
 117     bixy = el->FindElementValueAsNumberConvertTo("ixy", "SLUG*FT2");
 118   if (el->FindElement("ixz"))
 119     bixz = el->FindElementValueAsNumberConvertTo("ixz", "SLUG*FT2");
 120   if (el->FindElement("iyz"))
 121     biyz = el->FindElementValueAsNumberConvertTo("iyz", "SLUG*FT2");
 122   SetAircraftBaseInertias(FGMatrix33(  bixx,  -bixy,  bixz,
 123                                       -bixy,  biyy,  -biyz,
 124                                        bixz,  -biyz,  bizz ));
 125   if (el->FindElement("emptywt")) {
 126     EmptyWeight = el->FindElementValueAsNumberConvertTo("emptywt", "LBS");
 127   }
 128
 129   element = el->FindElement("location");
 130   while (element) {
 131     element_name = element->GetAttributeValue("name");
 132     if (element_name == "CG") vbaseXYZcg = element->FindElementTripletConvertTo("IN");
 133     element = el->FindNextElement("location");
 134   }
 135
 136 // Find all POINTMASS elements that descend from this METRICS branch of the
 137 // config file.
 138
 139   element = el->FindElement("pointmass");
 140   while (element) {
 141     AddPointMass(element);
 142     element = el->FindNextElement("pointmass");
 143   }
 144
 145   double ChildFDMWeight = 0.0;
 146   for (int fdm=0; fdm<FDMExec->GetFDMCount(); fdm++) {
 147     if (FDMExec->GetChildFDM(fdm)->mated) ChildFDMWeight += FDMExec->GetChildFDM(fdm)->exec->GetMassBalance()->GetWeight();
 148   }
 149
 150   Weight = EmptyWeight + in.TanksWeight + GetTotalPointMassWeight()
 151     + in.GasMass*slugtolb + ChildFDMWeight;
 152
 153   Mass = lbtoslug*Weight;
 154
 155   PostLoad(el, PropertyManager);
 156
 157   Debug(2);
 158   return true;
 159 }
 160
 161 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 162
 163 bool FGMassBalance::Run(bool Holding)
 164 {
 165   double denom, k1, k2, k3, k4, k5, k6;
 166   double Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
 167
 168   if (FGModel::Run(Holding)) return true;
 169   if (Holding) return false;
 170
 171   RunPreFunctions();
 172
 173   double ChildFDMWeight = 0.0;
 174   for (int fdm=0; fdm<FDMExec->GetFDMCount(); fdm++) {
 175     if (FDMExec->GetChildFDM(fdm)->mated) ChildFDMWeight += FDMExec->GetChildFDM(fdm)->exec->GetMassBalance()->GetWeight();
 176   }
 177
 178   Weight = EmptyWeight + in.TanksWeight + GetTotalPointMassWeight()
 179     + in.GasMass*slugtolb + ChildFDMWeight;
 180
 181   Mass = lbtoslug*Weight;
 182
 183 // Calculate new CG
 184
 185   vXYZcg = (EmptyWeight*vbaseXYZcg
 186             + GetPointMassMoment()
 187             + in.TanksMoment
 188             + in.GasMoment) / Weight;
 189
 190   // Track frame-by-frame delta CG, and move the EOM-tracked location
 191   // by this amount.
 192   if (vLastXYZcg.Magnitude() == 0.0) vLastXYZcg = vXYZcg;
 193   vDeltaXYZcg = vXYZcg - vLastXYZcg;
 194   vDeltaXYZcgBody = StructuralToBody(vLastXYZcg) - StructuralToBody(vXYZcg);
 195   vLastXYZcg = vXYZcg;
 196   FDMExec->GetPropagate()->NudgeBodyLocation(vDeltaXYZcgBody);
 197
 198 // Calculate new total moments of inertia
 199
 200   // At first it is the base configuration inertia matrix ...
 201   mJ = baseJ;
 202   // ... with the additional term originating from the parallel axis theorem.
 203   mJ += GetPointmassInertia( lbtoslug * EmptyWeight, vbaseXYZcg );
 204   // Then add the contributions from the additional pointmasses.
 205   mJ += CalculatePMInertias();
 206   mJ += in.TankInertia;
 207   mJ += in.GasInertia;
 208
 209   Ixx = mJ(1,1);
 210   Iyy = mJ(2,2);
 211   Izz = mJ(3,3);
 212   Ixy = -mJ(1,2);
 213   Ixz = -mJ(1,3);
 214   Iyz = -mJ(2,3);
 215
 216 // Calculate inertia matrix inverse (ref. Stevens and Lewis, "Flight Control & Simulation")
 217
 218   k1 = (Iyy*Izz - Iyz*Iyz);
 219   k2 = (Iyz*Ixz + Ixy*Izz);
 220   k3 = (Ixy*Iyz + Iyy*Ixz);
 221
 222   denom = 1.0/(Ixx*k1 - Ixy*k2 - Ixz*k3 );
 223   k1 = k1*denom;
 224   k2 = k2*denom;
 225   k3 = k3*denom;
 226   k4 = (Izz*Ixx - Ixz*Ixz)*denom;
 227   k5 = (Ixy*Ixz + Iyz*Ixx)*denom;
 228   k6 = (Ixx*Iyy - Ixy*Ixy)*denom;
 229
 230   mJinv.InitMatrix( k1, k2, k3,
 231                     k2, k4, k5,
 232                     k3, k5, k6 );
 233
 234   RunPostFunctions();
 235
 236   Debug(0);
 237
 238   return false;
 239 }
 240
 241 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 242
 243 void FGMassBalance::AddPointMass(Element* el)
 244 {
 245   double radius=0, length=0;
 246   Element* loc_element = el->FindElement("location");
 247   string pointmass_name = el->GetAttributeValue("name");
 248   if (!loc_element) {
 249     cerr << "Pointmass " << pointmass_name << " has no location." << endl;
 250     exit(-1);
 251   }
 252
 253   double w = el->FindElementValueAsNumberConvertTo("weight", "LBS");
 254   FGColumnVector3 vXYZ = loc_element->FindElementTripletConvertTo("IN");
 255
 256   PointMass *pm = new PointMass(w, vXYZ);
 257   pm->SetName(pointmass_name);
 258
 259   Element* form_element = el->FindElement("form");
 260   if (form_element) {
 261     string shape = form_element->GetAttributeValue("shape");
 262     Element* radius_element = form_element->FindElement("radius");
 263     Element* length_element = form_element->FindElement("length");
 264     if (radius_element) radius = form_element->FindElementValueAsNumberConvertTo("radius", "FT");
 265     if (length_element) length = form_element->FindElementValueAsNumberConvertTo("length", "FT");
 266     if (shape == "tube") {
 267       pm->SetPointMassShapeType(PointMass::esTube);
 268       pm->SetRadius(radius);
 269       pm->SetLength(length);
 270       pm->CalculateShapeInertia();
 271     } else if (shape == "cylinder") {
 272       pm->SetPointMassShapeType(PointMass::esCylinder);
 273       pm->SetRadius(radius);
 274       pm->SetLength(length);
 275       pm->CalculateShapeInertia();
 276     } else if (shape == "sphere") {
 277       pm->SetPointMassShapeType(PointMass::esSphere);
 278       pm->SetRadius(radius);
 279       pm->CalculateShapeInertia();
 280     } else if (shape == "ball") {
 281       pm->SetPointMassShapeType(PointMass::esBall);
 282       pm->SetRadius(radius);
 283       pm->CalculateShapeInertia();
 284     } else {
 285     }
 286   }
 287
 288   pm->bind(PropertyManager, PointMasses.size());
 289   PointMasses.push_back(pm);
 290 }
 291
 292 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 293
 294 double FGMassBalance::GetTotalPointMassWeight(void) const
 295 {
 296   double PM_total_weight = 0.0;
 297
 298   for (unsigned int i=0; i<PointMasses.size(); i++) {
 299     PM_total_weight += PointMasses[i]->Weight;
 300   }
 301   return PM_total_weight;
 302 }
 303
 304 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 305
 306 const FGColumnVector3& FGMassBalance::GetPointMassMoment(void)
 307 {
 308   PointMassCG.InitMatrix();
 309
 310   for (unsigned int i=0; i<PointMasses.size(); i++) {
 311     PointMassCG += PointMasses[i]->Weight*PointMasses[i]->Location;
 312   }
 313   return PointMassCG;
 314 }
 315
 316 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 317
 318 const FGMatrix33& FGMassBalance::CalculatePMInertias(void)
 319 {
 320   unsigned int size;
 321
 322   size = PointMasses.size();
 323   if (size == 0) return pmJ;
 324
 325   pmJ = FGMatrix33();
 326
 327   for (unsigned int i=0; i<size; i++) {
 328     pmJ += GetPointmassInertia( lbtoslug * PointMasses[i]->Weight, PointMasses[i]->Location );
 329     pmJ += PointMasses[i]->GetPointMassInertia();
 330   }
 331
 332   return pmJ;
 333 }
 334
 335 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 336
 337 FGColumnVector3 FGMassBalance::StructuralToBody(const FGColumnVector3& r) const
 338 {
 339   // Under the assumption that in the structural frame the:
 340   //
 341   // - X-axis is directed afterwards,
 342   // - Y-axis is directed towards the right,
 343   // - Z-axis is directed upwards,
 344   //
 345   // (as documented in http://jsbsim.sourceforge.net/JSBSimCoordinates.pdf)
 346   // we have to subtract first the center of gravity of the plane which
 347   // is also defined in the structural frame:
 348   //
 349   //   FGColumnVector3 cgOff = r - vXYZcg;
 350   //
 351   // Next, we do a change of units:
 352   //
 353   //   cgOff *= inchtoft;
 354   //
 355   // And then a 180 degree rotation is done about the Y axis so that the:
 356   //
 357   // - X-axis is directed forward,
 358   // - Y-axis is directed towards the right,
 359   // - Z-axis is directed downward.
 360   //
 361   // This is needed because the structural and body frames are 180 degrees apart.
 362
 363   return FGColumnVector3(inchtoft*(vXYZcg(1)-r(1)),
 364                          inchtoft*(r(2)-vXYZcg(2)),
 365                          inchtoft*(vXYZcg(3)-r(3)));
 366 }
 367
 368 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 369
 370 void FGMassBalance::bind(void)
 371 {
 372   typedef double (FGMassBalance::*PMF)(int) const;
 373   PropertyManager->Tie("inertia/mass-slugs", this,
 374                        &FGMassBalance::GetMass);
 375   PropertyManager->Tie("inertia/weight-lbs", this,
 376                        &FGMassBalance::GetWeight);
 377   PropertyManager->Tie("inertia/empty-weight-lbs", this,
 378                        &FGMassBalance::GetEmptyWeight);
 379   PropertyManager->Tie("inertia/cg-x-in", this,1,
 380                        (PMF)&FGMassBalance::GetXYZcg);
 381   PropertyManager->Tie("inertia/cg-y-in", this,2,
 382                        (PMF)&FGMassBalance::GetXYZcg);
 383   PropertyManager->Tie("inertia/cg-z-in", this,3,
 384                        (PMF)&FGMassBalance::GetXYZcg);
 385 }
 386
 387 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 388 //
 389 // This function binds properties for each pointmass object created.
 390 //
 391 void FGMassBalance::PointMass::bind(FGPropertyManager* PropertyManager, int num) {
 392   string tmp = CreateIndexedPropertyName("inertia/pointmass-weight-lbs", num);
 393   PropertyManager->Tie( tmp.c_str(), this, &PointMass::GetPointMassWeight,
 394                                        &PointMass::SetPointMassWeight);
 395
 396   tmp = CreateIndexedPropertyName("inertia/pointmass-location-X-inches", num);
 397   PropertyManager->Tie( tmp.c_str(), this, eX, &PointMass::GetPointMassLocation,
 398                                            &PointMass::SetPointMassLocation);
 399   tmp = CreateIndexedPropertyName("inertia/pointmass-location-Y-inches", num);
 400   PropertyManager->Tie( tmp.c_str(), this, eY, &PointMass::GetPointMassLocation,
 401                                            &PointMass::SetPointMassLocation);
 402   tmp = CreateIndexedPropertyName("inertia/pointmass-location-Z-inches", num);
 403   PropertyManager->Tie( tmp.c_str(), this, eZ, &PointMass::GetPointMassLocation,
 404                                            &PointMass::SetPointMassLocation);
 405 }
 406
 407 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 408
 409 void FGMassBalance::GetMassPropertiesReport(void) const
 410 {
 411   cout << endl << fgblue << highint
 412        << "  Mass Properties Report (English units: lbf, in, slug-ft^2)"
 413        << reset << endl;
 414   cout << "                                  " << underon << "    Weight    CG-X    CG-Y"
 415        << "    CG-Z         Ixx         Iyy         Izz" << underoff << endl;
 416   cout.precision(1);
 417   cout << highint << setw(34) << left << "    Base Vehicle " << normint
 418        << right << setw(10) << EmptyWeight << setw(8) << vbaseXYZcg(eX) << setw(8)
 419        << vbaseXYZcg(eY) << setw(8) << vbaseXYZcg(eZ) << setw(12) << baseJ(1,1)
 420        << setw(12) << baseJ(2,2) << setw(12) << baseJ(3,3) << endl;
 421
 422   for (unsigned int i=0;i<PointMasses.size();i++) {
 423     PointMass* pm = PointMasses[i];
 424     double pmweight = pm->GetPointMassWeight();
 425     cout << highint << left << setw(4) << i << setw(30) << pm->GetName() << normint
 426          << right << setw(10) << pmweight << setw(8) << pm->GetLocation()(eX)
 427          << setw(8) << pm->GetLocation()(eY) << setw(8) << pm->GetLocation()(eZ)
 428          << setw(12) << pm->GetPointMassMoI(1,1) << setw(12) << pm->GetPointMassMoI(2,2)
 429          << setw(12) << pm->GetPointMassMoI(3,3) << endl;
 430   }
 431
 432   cout << FDMExec->GetPropulsionTankReport();
 433
 434   cout << underon << setw(104) << " " << underoff << endl;
 435   cout << highint << left << setw(30) << "    Total: " << right << setw(14) << Weight
 436        << setw(8) << vXYZcg(eX)
 437        << setw(8) << vXYZcg(eY)
 438        << setw(8) << vXYZcg(eZ)
 439        << setw(12) << mJ(1,1)
 440        << setw(12) << mJ(2,2)
 441        << setw(12) << mJ(3,3)
 442        << normint << endl;
 443
 444   cout.setf(ios_base::fixed);
 445 }
 446
 447 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 448 //    The bitmasked value choices are as follows:
 449 //    unset: In this case (the default) JSBSim would only print
 450 //       out the normally expected messages, essentially echoing
 451 //       the config files as they are read. If the environment
 452 //       variable is not set, debug_lvl is set to 1 internally
 453 //    0: This requests JSBSim not to output any messages
 454 //       whatsoever.
 455 //    1: This value explicity requests the normal JSBSim
 456 //       startup messages
 457 //    2: This value asks for a message to be printed out when
 458 //       a class is instantiated
 459 //    4: When this value is set, a message is displayed when a
 460 //       FGModel object executes its Run() method
 461 //    8: When this value is set, various runtime state variables
 462 //       are printed out periodically
 463 //    16: When set various parameters are sanity checked and
 464 //       a message is printed out when they go out of bounds
 465
 466 void FGMassBalance::Debug(int from)
 467 {
 468   if (debug_lvl <= 0) return;
 469
 470   if (debug_lvl & 1) { // Standard console startup message output
 471     if (from == 2) { // Loading
 472       cout << endl << "  Mass and Balance:" << endl;
 473       cout << "    baseIxx: " << baseJ(1,1) << " slug-ft2" << endl;
 474       cout << "    baseIyy: " << baseJ(2,2) << " slug-ft2" << endl;
 475       cout << "    baseIzz: " << baseJ(3,3) << " slug-ft2" << endl;
 476       cout << "    baseIxy: " << baseJ(1,2) << " slug-ft2" << endl;
 477       cout << "    baseIxz: " << baseJ(1,3) << " slug-ft2" << endl;
 478       cout << "    baseIyz: " << baseJ(2,3) << " slug-ft2" << endl;
 479       cout << "    Empty Weight: " << EmptyWeight << " lbm" << endl;
 480       cout << "    CG (x, y, z): " << vbaseXYZcg << endl;
 481       // ToDo: Need to add point mass outputs here
 482       for (unsigned int i=0; i<PointMasses.size(); i++) {
 483         cout << "    Point Mass Object: " << PointMasses[i]->Weight << " lbs. at "
 484                    << "X, Y, Z (in.): " << PointMasses[i]->Location(eX) << "  "
 485                    << PointMasses[i]->Location(eY) << "  "
 486                    << PointMasses[i]->Location(eZ) << endl;
 487       }
 488     }
 489   }
 490   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 491     if (from == 0) cout << "Instantiated: FGMassBalance" << endl;
 492     if (from == 1) cout << "Destroyed:    FGMassBalance" << endl;
 493   }
 494   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 495   }
 496   if (debug_lvl & 8 ) { // Runtime state variables
 497   }
 498   if (debug_lvl & 16) { // Sanity checking
 499     if (from == 2) {
 500       if (EmptyWeight <= 0.0 || EmptyWeight > 1e9)
 501         cout << "MassBalance::EmptyWeight out of bounds: " << EmptyWeight << endl;
 502       if (Weight <= 0.0 || Weight > 1e9)
 503         cout << "MassBalance::Weight out of bounds: " << Weight << endl;
 504       if (Mass <= 0.0 || Mass > 1e9)
 505         cout << "MassBalance::Mass out of bounds: " << Mass << endl;
 506     }
 507   }
 508   if (debug_lvl & 64) {
 509     if (from == 0) { // Constructor
 510       cout << IdSrc << endl;
 511       cout << IdHdr << endl;
 512     }
 513   }
 514 }
 515 }