Updates to respiratory driver unit test.

Hardening respiratory compliance calculation to handle ~0 pressure.

engine/cpp/physiology/Respiratory.cpp

@@ -2136,6 +2136,7 @@ void Respiratory::UpdateChestWallCompliances()
 
     double minCompliance_L_Per_cmH2O = 0.01; //Minimum possible compliance
     double sideCompliance_L_Per_cmH2O = minCompliance_L_Per_cmH2O;
+    double approxZero = 1e-10;
     if (lungVolume_L > residualVolume_L && lungVolume_L < vitalCapacity_L + residualVolume_L)
     {
       double pressureCornerUpper_cmH2O = (vitalCapacity_L - functionalResidualCapacity_L) / healthySideCompliance_L_Per_cmH2O;
@@ -2144,7 +2145,15 @@ void Respiratory::UpdateChestWallCompliances()
       double d = (pressureCornerUpper_cmH2O - c) / 2.0;
       double expectedPressure_cmH2O = d * log((lungVolume_L - residualVolume_L) / (residualVolume_L + vitalCapacity_L - lungVolume_L)) + c;
       double volumeAtZeroPressure = residualVolume_L + (vitalCapacity_L / (1.0 + exp(c / d)));
-      sideCompliance_L_Per_cmH2O = (lungVolume_L - volumeAtZeroPressure) / expectedPressure_cmH2O;
+
+      sideCompliance_L_Per_cmH2O = healthySideCompliance_L_Per_cmH2O;
+
+      if (!(expectedPressure_cmH2O < approxZero && expectedPressure_cmH2O > -approxZero))
+      {
+        //Not dividing by ~0
+        sideCompliance_L_Per_cmH2O = (lungVolume_L - volumeAtZeroPressure) / expectedPressure_cmH2O;
+      }
+
       if (sideCompliance_L_Per_cmH2O == 0.0)
       {
         sideCompliance_L_Per_cmH2O = healthySideCompliance_L_Per_cmH2O;
@@ -3023,8 +3032,8 @@ void Respiratory::Debugging(SEFluidCircuit& RespirationCircuit)
   double leftPleuralPressure = m_LeftPleural->GetNextPressure(PressureUnit::cmH2O);
   double leftPleuralVolume = m_LeftPleural->GetNextVolume(VolumeUnit::L);
   double leftFlow = m_LeftAlveolarDeadSpaceToLeftAlveoli->GetNextFlow(VolumePerTimeUnit::L_Per_s);
-  double leftChestWallCompliance_L_Per_cmH2O = m_RightPleuralToRespiratoryMuscle->GetNextCompliance(VolumePerPressureUnit::L_Per_cmH2O);
-  double leftLungCompliance_L_Per_cmH2O = m_RightAlveoliToRightPleuralConnection->GetNextCompliance(VolumePerPressureUnit::L_Per_cmH2O);
+  double leftChestWallCompliance_L_Per_cmH2O = m_LeftPleuralToRespiratoryMuscle->GetNextCompliance(VolumePerPressureUnit::L_Per_cmH2O);
+  double leftLungCompliance_L_Per_cmH2O = m_LeftAlveoliToLeftPleuralConnection->GetNextCompliance(VolumePerPressureUnit::L_Per_cmH2O);
 
   double rightAlveoliPressure = m_RightAlveoli->GetNextPressure(PressureUnit::cmH2O);
   double rightAlveoliVolume = m_RightAlveoli->GetNextVolume(VolumeUnit::L);
@@ -3057,6 +3066,6 @@ void Respiratory::Debugging(SEFluidCircuit& RespirationCircuit)
   double leftSideCompliance_L_Per_cmH2O = 1.0 / (1.0 / leftChestWallCompliance_L_Per_cmH2O + 1.0 / leftLungCompliance_L_Per_cmH2O);
   double rightSideCompliance_L_Per_cmH2O = 1.0 / (1.0 / rightChestWallCompliance_L_Per_cmH2O + 1.0 / rightLungCompliance_L_Per_cmH2O);
 
-  double total_Compliance_L_Per_cmH2O = leftSideCompliance_L_Per_cmH2O + rightSideCompliance_L_Per_cmH2O;
-  m_data.GetDataTrack().Probe("total_Compliance_L_Per_cmH2O", total_Compliance_L_Per_cmH2O);
+  double totalCompliance_L_Per_cmH2O = leftSideCompliance_L_Per_cmH2O + rightSideCompliance_L_Per_cmH2O;
+  m_data.GetDataTrack().Probe("totalCompliance_L_Per_cmH2O", totalCompliance_L_Per_cmH2O);
 }
\ No newline at end of file

test/engine/cpp/RespiratoryCircuit.cpp

@@ -252,6 +252,7 @@ void PulseEngineTest::RespiratoryDriverTest(const std::string& sTestDirectory)
   bool bRVReached = false;
   int iTime = 0;
   double AmbientPresure_cmH2O = 1033.23; // = 1 atm
+  double calculatedTotalCompliance_L_Per_cmH2O = 0.0;
 
   RespCircuit.GetNode(pulse::EnvironmentNode::Ambient)->GetNextPressure().SetValue(AmbientPresure_cmH2O, PressureUnit::cmH2O);
 
@@ -290,28 +291,53 @@ void PulseEngineTest::RespiratoryDriverTest(const std::string& sTestDirectory)
         if (bRVReached)
         {
           //Output values
+          double driverPressure = driverPressurePath->GetPressureSource(PressureUnit::cmH2O);
+
           trk1.Track("LungVolume_L", iTime, TotalVolume_L);
-          trk1.Track("DriverPressure_cmH2O", iTime, driverPressurePath->GetPressureSource(PressureUnit::cmH2O));
+          trk1.Track("DriverPressure_cmH2O", iTime, driverPressure);
+
+          double leftAlveoliPressure = leftAlveoliNode->GetPressure(PressureUnit::cmH2O);
+          double leftAlveoliVolume = leftAlveoliNode->GetVolume(VolumeUnit::L);
+          double leftPleuralPressure = leftPleuralNode->GetPressure(PressureUnit::cmH2O);
+          double leftPleuralVolume = leftPleuralNode->GetVolume(VolumeUnit::L);
+          double leftChestWallCompliance_L_Per_cmH2O = leftPleuralToRespiratoryMuscle->GetCompliance(VolumePerPressureUnit::L_Per_cmH2O);
+          double leftLungCompliance_L_Per_cmH2O = leftPleuralToRespiratoryMuscle->GetCompliance(VolumePerPressureUnit::L_Per_cmH2O);
+
+          double rightAlveoliPressure = rightAlveoliNode->GetPressure(PressureUnit::cmH2O);
+          double rightAlveoliVolume = rightAlveoliNode->GetVolume(VolumeUnit::L);
+          double rightPleuralPressure = rightPleuralNode->GetPressure(PressureUnit::cmH2O);
+          double rightPleuralVolume = rightPleuralNode->GetVolume(VolumeUnit::L);
+          double rightChestWallCompliance_L_Per_cmH2O = rightPleuralToRespiratoryMuscle->GetCompliance(VolumePerPressureUnit::L_Per_cmH2O);
+          double rightLungCompliance_L_Per_cmH2O = rightPleuralToRespiratoryMuscle->GetCompliance(VolumePerPressureUnit::L_Per_cmH2O);
+
+          double leftSideCompliance_L_Per_cmH2O = 1.0 / (1.0 / leftChestWallCompliance_L_Per_cmH2O + 1.0 / leftLungCompliance_L_Per_cmH2O);
+          double rightSideCompliance_L_Per_cmH2O = 1.0 / (1.0 / rightChestWallCompliance_L_Per_cmH2O + 1.0 / rightLungCompliance_L_Per_cmH2O);
 
-          //Plug the mouth and relax the muscles (just like the real life test)
-          mouthToCarina->GetNextResistance().SetValue(10e100, PressureTimePerVolumeUnit::cmH2O_s_Per_L);
-          driverPressurePath->GetNextPressureSource().SetValue(0.0, PressureUnit::cmH2O);
-          //Process the circuit
-          calc.Process(RespCircuit, deltaT_s);
-          //Advance time
-          calc.PostProcess(RespCircuit);
+          double totalCompliance_L_Per_cmH2O = leftSideCompliance_L_Per_cmH2O + rightSideCompliance_L_Per_cmH2O;
+
+          double approxZero = 1e-10;
+          if (!(driverPressure < approxZero && driverPressure > -approxZero))
+          {
+            calculatedTotalCompliance_L_Per_cmH2O = -(TotalVolume_L - pc.GetCurrentPatient().GetFunctionalResidualCapacity(VolumeUnit::L)) / driverPressure;
+          }           
 
           trk1.Track("ChestWallPressure_cmH2O", iTime, muscleNode->GetPressure(PressureUnit::cmH2O) - AmbientPresure_cmH2O);
-          trk1.Track("RightAlveoliVolume_L", iTime, rightAlveoliNode->GetVolume(VolumeUnit::L));
-          trk1.Track("LeftAlveoliVolume_L", iTime, leftAlveoliNode->GetVolume(VolumeUnit::L));
-          trk1.Track("RightAlveoliPressure_cmH2O", iTime, rightAlveoliNode->GetPressure(PressureUnit::cmH2O) - AmbientPresure_cmH2O);
-          trk1.Track("LeftAlveoliPressure_cmH2O", iTime, leftAlveoliNode->GetPressure(PressureUnit::cmH2O) - AmbientPresure_cmH2O);
-          trk1.Track("RightPleuralVolume_L", iTime, rightPleuralNode->GetVolume(VolumeUnit::L));
-          trk1.Track("LeftPleuralVolume_L", iTime, leftPleuralNode->GetVolume(VolumeUnit::L));
-          trk1.Track("RightPleuralPressure_cmH2O", iTime, rightPleuralNode->GetPressure(PressureUnit::cmH2O) - AmbientPresure_cmH2O);
-          trk1.Track("LeftPleuralPressure_cmH2O", iTime, leftPleuralNode->GetPressure(PressureUnit::cmH2O) - AmbientPresure_cmH2O);
-          trk1.Track("rightTranspulmonaryPressure_cmH2O", iTime, rightAlveoliNode->GetPressure(PressureUnit::cmH2O) - rightPleuralNode->GetPressure(PressureUnit::cmH2O));
-          trk1.Track("LeftTranspulmonaryPressure_cmH2O", iTime, leftAlveoliNode->GetPressure(PressureUnit::cmH2O) - leftPleuralNode->GetPressure(PressureUnit::cmH2O));
+          trk1.Track("RightAlveoliVolume_L", iTime, rightAlveoliVolume);
+          trk1.Track("LeftAlveoliVolume_L", iTime, leftAlveoliVolume);
+          trk1.Track("RightAlveoliPressure_cmH2O", iTime, rightAlveoliPressure - AmbientPresure_cmH2O);
+          trk1.Track("LeftAlveoliPressure_cmH2O", iTime, leftAlveoliPressure - AmbientPresure_cmH2O);
+          trk1.Track("RightPleuralVolume_L", iTime, rightPleuralVolume);
+          trk1.Track("LeftPleuralVolume_L", iTime, leftPleuralVolume);
+          trk1.Track("RightPleuralPressure_cmH2O", iTime, rightPleuralPressure - AmbientPresure_cmH2O);
+          trk1.Track("LeftPleuralPressure_cmH2O", iTime, leftPleuralPressure - AmbientPresure_cmH2O);
+          trk1.Track("RightTranspulmonaryPressure_cmH2O", iTime, rightAlveoliPressure - rightPleuralPressure);
+          trk1.Track("LeftTranspulmonaryPressure_cmH2O", iTime, leftAlveoliPressure - leftPleuralPressure);
+
+          trk1.Track("leftSideCompliance_L_Per_cmH2O", iTime, leftSideCompliance_L_Per_cmH2O);
+          trk1.Track("rightSideCompliance_L_Per_cmH2O", iTime, rightSideCompliance_L_Per_cmH2O);
+          trk1.Track("totalCompliance_L_Per_cmH2O", iTime, totalCompliance_L_Per_cmH2O);
+
+          trk1.Track("calculatedTotalCompliance_L_Per_cmH2O", iTime, calculatedTotalCompliance_L_Per_cmH2O);
 
           iTime++;
         }