text: /* Distributed under the Apache License, Version 2.0.
     See accompanying NOTICE file for details.*/
  
  #include "EngineTest.h"
  #include "Controller/Controller.h"
  #include "utils/TimingProfile.h"
  #include "utils/testing/SETestReport.h"
  #include "utils/testing/SETestCase.h"
  #include "utils/testing/SETestSuite.h"
  
  #include "patient/SEPatient.h"
  #include "circuit/fluid/SEFluidCircuit.h"
  #include "substance/SESubstanceFraction.h"
  #include "compartment/SECompartmentManager.h"
  #include "compartment/fluid/SEGasCompartmentGraph.h"
  #include "compartment/fluid/SELiquidCompartmentGraph.h"
  #include "compartment/substances/SEGasSubstanceQuantity.h"
  #include "properties/SEScalarPressure.h"
  #include "properties/SEScalarVolume.h"
  #include "properties/SEScalarVolumePerTime.h"
  #include "properties/SEScalarFlowCompliance.h"
  #include "properties/SEScalarTime.h"
  #include "properties/SEScalarInverseVolume.h"
  #include "properties/SEScalarMassPerVolume.h"
  #include "properties/SEHistogramFractionVsLength.h"
  #include "properties/SEScalar0To1.h"
  #include "properties/SEScalarNegative1To1.h"
  #include "properties/SEScalarLength.h"
  
  void PulseEngineTest::AerosolTest(const std::string& sOutputDirectory)
  {
    m_Logger->ResetLogFile(sOutputDirectory + "/AerosolTest.log");
  
    SETestReport testReport = SETestReport(m_Logger);
  
    //// Create a suite
    //SETestSuite& albuteroluite = testReport.CreateTestSuite();
    //albuteroluite.SetName("Albuterol");
    //SESubstanceManager subMgr(m_Logger);
    //subMgr.LoadSubstanceDirectory();
    //SESubstance* albuterol = subMgr.GetSubstance("Albuterol");
    //SizeIndependentDepositionEfficencyCoefficientsTest(albuteroluite, *albuterol, 0.043737, 0.045737, 0.090432, 0.3115);
    //DepositionFractionTest(albuteroluite, *albuterol, 0.0805701, 0.0843696, 0.0701852, 0.111879);
  
    // Create a suite
    SETestSuite& normalDistributionSuite = testReport.CreateTestSuite();
    normalDistributionSuite.SetName("NormalDistribution");
    SESubstance normalDistributedSubstance(m_Logger);
    normalDistributedSubstance.GetAerosolization().GetBronchioleModifier().SetValue(0);//Need something here...
    normalDistributedSubstance.GetAerosolization().GetInflammationCoefficient().SetValue(0.5);//Need something here...
    SEHistogramFractionVsLength& concentrations = normalDistributedSubstance.GetAerosolization().GetParticulateSizeDistribution();
    // Boundary -1 is 0, i.e. no diameter
    concentrations.GetLength().push_back(1.0e-4); // Boundary 0
    concentrations.GetLength().push_back(1.0e-3); // Boundary 1
    concentrations.GetLength().push_back(1.0e-2); // Boundary 2
    concentrations.GetLength().push_back(1.0e-1); // Boundary 3
    concentrations.GetLength().push_back(1.0);    // Boundary 4
    concentrations.GetLength().push_back(1.0e1);  // Boundary 5
    concentrations.GetLength().push_back(1.0e2);  // Boundary 6
                                                  // Boundary n is infinity
    concentrations.SetLengthUnit(LengthUnit::um);
  
                                                   // Bin -1 contains all of the particles that are less than 1.0e-4 um aerodynamic diameter
    concentrations.GetFraction().push_back(0.015); // Bin 0 (Bounded by Boundary 0 and Boundary 1)
    concentrations.GetFraction().push_back(0.035); // Bin 1 (Bounded by Boundary 1 and Boundary 2)
    concentrations.GetFraction().push_back(0.9);   // Bin 2 (Bounded by Boundary 2 and Boundary 3)
    concentrations.GetFraction().push_back(0.035); // Bin 3 (Bounded by Boundary 3 and Boundary 4)
    concentrations.GetFraction().push_back(0.015); // Bin 4 (Bounded by Boundary 4 and Boundary 5)
    concentrations.GetFraction().push_back(0);     // Bin 5 (Bounded by Boundary 5 and Boundary 6)
                                                   // Bin n contains 
  
    SizeIndependentDepositionEfficencyCoefficientsTest(normalDistributionSuite, normalDistributedSubstance, 0.043737,0.045737,0.090432,0.3115);
    DepositionFractionTest(normalDistributionSuite, normalDistributedSubstance, 0.0805701, 0.0843696, 0.0701852, 0.111879);
  
    // Create another suite...
    // Create a suite
    SETestSuite& monodispersedSuite = testReport.CreateTestSuite();
    monodispersedSuite.SetName("Monodispersed");
    SESubstance monodispersedSubstance(m_Logger);
    monodispersedSubstance.GetAerosolization().GetBronchioleModifier().SetValue(0);//Need something here...
    monodispersedSubstance.GetAerosolization().GetInflammationCoefficient().SetValue(0.5);//Need something here...
    SEHistogramFractionVsLength& monoConcentrations = monodispersedSubstance.GetAerosolization().GetParticulateSizeDistribution();
    // Boundary -1 is 0, i.e. no diameter
    monoConcentrations.GetLength().push_back(1.0e-4); // Boundary 0
    monoConcentrations.GetLength().push_back(1.0e-3); // Boundary 1
    monoConcentrations.GetLength().push_back(1.0e-2); // Boundary 2
    monoConcentrations.GetLength().push_back(1.0e-1); // Boundary 3
    monoConcentrations.GetLength().push_back(1.0);    // Boundary 4
    monoConcentrations.GetLength().push_back(1.0e1);  // Boundary 5
    monoConcentrations.GetLength().push_back(1.0e2);  // Boundary 6
                                                  // Boundary n is infinity
    monoConcentrations.SetLengthUnit(LengthUnit::um);
  
    monoConcentrations.GetFraction().push_back(0);     // Bin 0 (Bounded by Boundary 0 and Boundary 1)
    monoConcentrations.GetFraction().push_back(0);     // Bin 1 (Bounded by Boundary 1 and Boundary 2)
    monoConcentrations.GetFraction().push_back(1.0);   // Bin 2 (Bounded by Boundary 2 and Boundary 3)
    monoConcentrations.GetFraction().push_back(0);     // Bin 3 (Bounded by Boundary 3 and Boundary 4)
    monoConcentrations.GetFraction().push_back(0);     // Bin 4 (Bounded by Boundary 4 and Boundary 5)
    monoConcentrations.GetFraction().push_back(0);     // Bin 5 (Bounded by Boundary 5 and Boundary 6)
  
    SizeIndependentDepositionEfficencyCoefficientsTest(monodispersedSuite, monodispersedSubstance, 0.024864, 0.025498, 0.083839, 0.32856);
    DepositionFractionTest(monodispersedSuite, monodispersedSubstance, 0.0471459, 0.0482577, 0.0660969, 0.117988);
  
    // Create a suite
    SETestSuite& mono2Suite = testReport.CreateTestSuite();
    mono2Suite.SetName("Monodispersed_2");
    SESubstance mono2Substance(m_Logger);
    mono2Substance.GetAerosolization().GetBronchioleModifier().SetValue(0);//Need something here...
    mono2Substance.GetAerosolization().GetInflammationCoefficient().SetValue(0.5);//Need something here...
    SEHistogramFractionVsLength& mono2Concentrations = mono2Substance.GetAerosolization().GetParticulateSizeDistribution();
    // Boundary -1 is 0, i.e. no diameter
    mono2Concentrations.GetLength().push_back(1.0e-4); // Boundary 0
    mono2Concentrations.GetLength().push_back(1.0e-3); // Boundary 1
    mono2Concentrations.GetLength().push_back(1.0e-2); // Boundary 2
    mono2Concentrations.GetLength().push_back(1.0e-1); // Boundary 3
    mono2Concentrations.GetLength().push_back(1.0);    // Boundary 4
    mono2Concentrations.GetLength().push_back(1.0e1);  // Boundary 5
    mono2Concentrations.GetLength().push_back(1.0e2);  // Boundary 6
                                                      // Boundary n is infinity
    mono2Concentrations.SetLengthUnit(LengthUnit::um);
  
    mono2Concentrations.GetFraction().push_back(0);     // Bin 0 (Bounded by Boundary 0 and Boundary 1)
    mono2Concentrations.GetFraction().push_back(0);     // Bin 1 (Bounded by Boundary 1 and Boundary 2)
    mono2Concentrations.GetFraction().push_back(0);   // Bin 2 (Bounded by Boundary 2 and Boundary 3)
    mono2Concentrations.GetFraction().push_back(0);     // Bin 3 (Bounded by Boundary 3 and Boundary 4)
    mono2Concentrations.GetFraction().push_back(0);     // Bin 4 (Bounded by Boundary 4 and Boundary 5)
    mono2Concentrations.GetFraction().push_back(1.0);     // Bin 5 (Bounded by Boundary 5 and Boundary 6)
  
    SizeIndependentDepositionEfficencyCoefficientsTest(mono2Suite, mono2Substance, 0.26841, 0.36689, 0.0012658, 0.0014995);
    DepositionFractionTest(mono2Suite, mono2Substance, 0.366641, 0.504666, 0.000931153, 0.000542435);
  
    // Create a suite
    SETestSuite& mono3Suite = testReport.CreateTestSuite();
    mono3Suite.SetName("Monodispersed_3");
    SESubstance mono3Substance(m_Logger);
    mono3Substance.GetAerosolization().GetBronchioleModifier().SetValue(0);//Need something here...
    mono3Substance.GetAerosolization().GetInflammationCoefficient().SetValue(0.5);//Need something here...
    SEHistogramFractionVsLength& mono3Concentrations = mono3Substance.GetAerosolization().GetParticulateSizeDistribution();
    // Boundary -1 is 0, i.e. no diameter
    mono3Concentrations.GetLength().push_back(1.0e-4); // Boundary 0
    mono3Concentrations.GetLength().push_back(1.0e-3); // Boundary 1
    mono3Concentrations.GetLength().push_back(1.0e-2); // Boundary 2
    mono3Concentrations.GetLength().push_back(1.0e-1); // Boundary 3
    mono3Concentrations.GetLength().push_back(1.0);    // Boundary 4
    mono3Concentrations.GetLength().push_back(1.0e1);  // Boundary 5
    mono3Concentrations.GetLength().push_back(1.0e2);  // Boundary 6
                                                       // Boundary n is infinity
    mono3Concentrations.SetLengthUnit(LengthUnit::um);
  
    mono3Concentrations.GetFraction().push_back(1.0);   // Bin 0 (Bounded by Boundary 0 and Boundary 1)
    mono3Concentrations.GetFraction().push_back(0);     // Bin 1 (Bounded by Boundary 1 and Boundary 2)
    mono3Concentrations.GetFraction().push_back(0);     // Bin 2 (Bounded by Boundary 2 and Boundary 3)
    mono3Concentrations.GetFraction().push_back(0);     // Bin 3 (Bounded by Boundary 3 and Boundary 4)
    mono3Concentrations.GetFraction().push_back(0);     // Bin 4 (Bounded by Boundary 4 and Boundary 5)
    mono3Concentrations.GetFraction().push_back(0);     // Bin 5 (Bounded by Boundary 5 and Boundary 6)
  
    SizeIndependentDepositionEfficencyCoefficientsTest(mono3Suite, mono3Substance, 0.43617, 0.77358, 0.21788, 0.005489);
    DepositionFractionTest(mono3Suite, mono3Substance, 0.468578, 0.753432, 0.1188, 0.00197536);
  
    // Create a suite
    SETestSuite& zhangDispersion = testReport.CreateTestSuite();
    zhangDispersion.SetName("ZhangDispersion");
    SESubstance zhangSubstance(m_Logger);
    zhangSubstance.GetAerosolization().GetBronchioleModifier().SetValue(0);//Need something here...
    zhangSubstance.GetAerosolization().GetInflammationCoefficient().SetValue(0.5);//Need something here...
    SEHistogramFractionVsLength& zhangConcentrations = zhangSubstance.GetAerosolization().GetParticulateSizeDistribution();
    // This histogram based on a digitization of Figure 4b in @cite zhang2012chemical
    zhangConcentrations.GetLength().push_back(10.308639); // Boundary 0
    zhangConcentrations.GetLength().push_back(16.263765); // Boundary 1
    zhangConcentrations.GetLength().push_back(25.854805); // Boundary 2
    zhangConcentrations.GetLength().push_back(41.101852); // Boundary 3
    zhangConcentrations.GetLength().push_back(65.83879);    // Boundary 4
    zhangConcentrations.GetLength().push_back(103.87276);  // Boundary 5
    zhangConcentrations.GetLength().push_back(165.1284);  // Boundary 6
    zhangConcentrations.GetLength().push_back(260.5203);  // Boundary 7
    zhangConcentrations.GetLength().push_back(417.31314);  // Boundary 8
    zhangConcentrations.GetLength().push_back(663.4103);  // Boundary 9
    zhangConcentrations.GetLength().push_back(1062.6803);  // Boundary 10
                                                       // Boundary n is infinity
    zhangConcentrations.SetLengthUnit(LengthUnit::um);
  
    zhangConcentrations.GetFraction().push_back(0.0028992);       // Bin 0 (Bounded by Boundary 0 and Boundary 1)
    zhangConcentrations.GetFraction().push_back(0.018143616);     // Bin 1 (Bounded by Boundary 1 and Boundary 2)
    zhangConcentrations.GetFraction().push_back(0.069563083);     // Bin 2 (Bounded by Boundary 2 and Boundary 3)
    zhangConcentrations.GetFraction().push_back(0.170665942);     // Bin 3 (Bounded by Boundary 3 and Boundary 4)
    zhangConcentrations.GetFraction().push_back(0.251295383);     // Bin 4 (Bounded by Boundary 4 and Boundary 5)
    zhangConcentrations.GetFraction().push_back(0.245390533);     // Bin 5 (Bounded by Boundary 5 and Boundary 6)
    zhangConcentrations.GetFraction().push_back(0.145956494);     // Bin 6 (Bounded by Boundary 6 and Boundary 7)
    zhangConcentrations.GetFraction().push_back(0.058885618);     // Bin 7 (Bounded by Boundary 7 and Boundary 8)
    zhangConcentrations.GetFraction().push_back(0.020538896);     // Bin 8 (Bounded by Boundary 8 and Boundary 9)
    zhangConcentrations.GetFraction().push_back(0.016661235);     // Bin 9 (Bounded by Boundary 9 and Boundary 10)
  
    SizeIndependentDepositionEfficencyCoefficientsTest(zhangDispersion, zhangSubstance, 0.25368, 0.3399, 0.00013825, 0.00022882);
    DepositionFractionTest(zhangDispersion, zhangSubstance, 0.353503, 0.478626, 0.000102784, 8.27909e-05);
  
    testReport.WriteFile(sOutputDirectory + "/AerosolTestReport.pba");
  }
  
  void PulseEngineTest::SizeIndependentDepositionEfficencyCoefficientsTest(SETestSuite& suite, SESubstance& substance, 
                                                                              double expectedMouthCoeff, double expectedCarinaCoeff, double expectedDeadSpaceCoeff, double expectedAlveoliCoeff)
  {  
    double PercentTolerance = 0.1;
    // Set up our test report
    TimingProfile pTimer;
    pTimer.Start("Test");
    SETestCase& tc1 = suite.CreateTestCase();
    tc1.SetName(substance.GetName()+"SIDECo");
    
    PulseController pc(m_Logger);
    const SizeIndependentDepositionEfficencyCoefficient& SIDECoeff =  pc.GetSubstances().GetSizeIndependentDepositionEfficencyCoefficient(substance);
    m_ss << "Mouth: " << SIDECoeff.GetMouth();
    Info(m_ss);
    if (GeneralMath::PercentTolerance(SIDECoeff.GetMouth(), expectedMouthCoeff) > PercentTolerance)
    {
      m_ss << " Carina SIDE Coefficient is not correct : " << SIDECoeff.GetMouth() << " expected " << expectedMouthCoeff;
      tc1.AddFailure(m_ss);
    }
    m_ss << "Carina: " << SIDECoeff.GetCarina();
    Info(m_ss);
    if (GeneralMath::PercentTolerance(SIDECoeff.GetCarina(), expectedCarinaCoeff) > PercentTolerance)
    {
      m_ss << " Deadspace SIDE Coefficient is not correct : " << SIDECoeff.GetCarina() << " expected " << expectedCarinaCoeff;
      tc1.AddFailure(m_ss);
    }
    m_ss << "Deadspace: " << SIDECoeff.GetDeadSpace();
    Info(m_ss);
    if (GeneralMath::PercentTolerance(SIDECoeff.GetDeadSpace(), expectedDeadSpaceCoeff) > PercentTolerance)
    {
      m_ss << " Mouth SIDE Coefficient is not correct : " << SIDECoeff.GetDeadSpace() << " expected " << expectedDeadSpaceCoeff;
      tc1.AddFailure(m_ss);
    }
    m_ss << "Alveoli: " << SIDECoeff.GetAlveoli();
    Info(m_ss);
    if (GeneralMath::PercentTolerance(SIDECoeff.GetAlveoli(), expectedAlveoliCoeff) > PercentTolerance)
    {
      m_ss << " Alveoli SIDE Coefficient is not correct : " << SIDECoeff.GetAlveoli() << " expected " << expectedAlveoliCoeff;
      tc1.AddFailure(m_ss);
    }
    tc1.GetDuration().SetValue(pTimer.GetElapsedTime_s("Test"), TimeUnit::s); 
  }
  
  void PulseEngineTest::DepositionFractionTest(SETestSuite& suite, SESubstance& substance,
    double expectedMouthDepFrac, double expectedCarinaDepFrac, double expectedDeadSpaceDepFrac, double expectedAlveoliDepFrac)
  {
    TimingProfile pTimer;
    pTimer.Start("Test");
    double PercentTolerance = 0.1;
    SETestCase& tc = suite.CreateTestCase();
    tc.SetName(substance.GetName()+"DepositionFraction");
  
    PulseController pc(m_Logger);
    pc.GetPatient().LoadFile("./patients/StandardMale.pba");
    pc.SetupPatient();
    pc.m_Config->EnableRenal(cdm::eSwitch::Off);
    pc.m_Config->EnableTissue(cdm::eSwitch::Off);
    pc.CreateCircuitsAndCompartments();
    pc.GetSubstances().InitializeGasCompartments();
    SEEnvironmentalConditions& env = pc.GetEnvironment().GetConditions(); 
    SELiquidCompartment* ambient = pc.GetCompartments().GetLiquidCompartment(pulse::EnvironmentCompartment::Ambient);
  
    SEFluidCircuit* rCircuit = &pc.GetCircuits().GetRespiratoryCircuit();
    SEGasCompartmentGraph* rGraph = &pc.GetCompartments().GetRespiratoryGraph();
    SELiquidCompartmentGraph* aGraph = &pc.GetCompartments().GetAerosolGraph();
  
    const SizeIndependentDepositionEfficencyCoefficient& SIDECoeff = pc.GetSubstances().GetSizeIndependentDepositionEfficencyCoefficient(substance);
    
    
    SEFluidCircuitPath *driverPath = rCircuit->GetPath(pulse::RespiratoryPath::EnvironmentToRespiratoryMuscle);
    SEGasTransporter    gtxpt(VolumePerTimeUnit::L_Per_s, VolumeUnit::L, VolumeUnit::L, NoUnit::unitless, pc.GetLogger());
    SELiquidTransporter ltxpt(VolumePerTimeUnit::mL_Per_s, VolumeUnit::mL, MassUnit::ug, MassPerVolumeUnit::ug_Per_mL, pc.GetLogger());
    SEFluidCircuitCalculator calc(FlowComplianceUnit::L_Per_cmH2O, VolumePerTimeUnit::L_Per_s, FlowInertanceUnit::cmH2O_s2_Per_L, PressureUnit::cmH2O, VolumeUnit::L, FlowResistanceUnit::cmH2O_s_Per_L, pc.GetLogger());
  
    //Set the reference not pressure to the standard environment
    //This is needed because we're not setting the Environment during initialization in this unit test
    rCircuit->GetNode(pulse::EnvironmentNode::Ambient)->GetNextPressure().Set(env.GetAtmosphericPressure());
    rCircuit->GetNode(pulse::EnvironmentNode::Ambient)->GetPressure().Set(env.GetAtmosphericPressure());
    //Precharge the stomach to prevent negative volume
    rCircuit->GetNode(pulse::RespiratoryNode::Stomach)->GetNextPressure().Set(env.GetAtmosphericPressure());
    rCircuit->GetNode(pulse::RespiratoryNode::Stomach)->GetPressure().Set(env.GetAtmosphericPressure());
  
    SEFluidCircuitPath* Env2Mouth = rCircuit->GetPath(pulse::RespiratoryPath::EnvironmentToMouth);
  
    SELiquidSubstanceQuantity* mouthParticulate = nullptr;
    SELiquidSubstanceQuantity* carinaParticulate = nullptr;
    SELiquidSubstanceQuantity* leftDeadSpaceParticulate = nullptr;
    SELiquidSubstanceQuantity* leftAlveoliParticulate = nullptr;
    SELiquidSubstanceQuantity* rightDeadSpaceParticulate = nullptr;
    SELiquidSubstanceQuantity* rightAlveoliParticulate = nullptr;
  
    //Circuit Analysis Test --------------------------------------------------
    //Execution parameters
    double time = 0;
    double deltaT_s = 1.0 / 50.0;
   
    //Drive waveform parameters
    double period = 5.0; 
    double runTime_s = period*2;
    double alpha = (2 * 3.1415) / (period);
    double driverPressure_cmH2O = 0.0;
    double amplitude_cmH2O = 6.0;
    double yOffset = -12.0;
    double sinusoid = 0;
  
    DataTrack trk;
    std::ofstream file;
    double totalInspiredAir_mL = 0;
    double totalExspiredAir_mL = 0;
    double totalInspiredParticulate_ug=0;
  
    double totalParticulateDeposited_ug = 0;
    double depositedMouthParticulate_ug = 0;
    double depositedCarinaParticulate_ug = 0;
    double depositedLeftDeadSpaceParticulate_ug = 0;
    double depositedLeftAlveoliParticulate_ug = 0;
    double depositedRightDeadSpaceParticulate_ug = 0;
    double depositedRightAlveoliParticulate_ug = 0;
  
    for (unsigned int i = 0; i < runTime_s / deltaT_s; i++)
    {
      sinusoid = cos(alpha * time);
      //PreProcess - Push driver pressure and  variable compliance data into the Circuit
      driverPressure_cmH2O = yOffset + amplitude_cmH2O * sinusoid;   //compute new pressure
      driverPath->GetNextPressureSource().SetValue(driverPressure_cmH2O, PressureUnit::cmH2O);
  
      //Process - Execute the circuit
      calc.Process(*rCircuit, deltaT_s);
      //Execute the substance transport function
      gtxpt.Transport(*rGraph, deltaT_s);
      // Do it again for aerosols
      if (aGraph != nullptr)
        ltxpt.Transport(*aGraph, deltaT_s);
      //convert 'Next' values to current
      calc.PostProcess(*rCircuit);
  
      trk.Track("DriverPressure_cmH2O", time, driverPressure_cmH2O);
      trk.Track("Env2Mouth_mL/s", time, Env2Mouth->GetFlow(VolumePerTimeUnit::mL_Per_s));
      
      trk.Track("TotalMouthPreParticulate_ug", time, mouthParticulate == nullptr ? 0 : mouthParticulate->GetMass(MassUnit::ug));
      trk.Track("MouthAerosol_mL", time, aGraph->GetCompartment(pulse::PulmonaryCompartment::Mouth)->GetVolume(VolumeUnit::mL));
      trk.Track("MouthAir_mL", time, rGraph->GetCompartment(pulse::PulmonaryCompartment::Mouth)->GetVolume(VolumeUnit::mL));
  
      if (i > runTime_s*0.5 / deltaT_s)
      {// Completed a cycle now start tracking Deposition and total inspired Air
  
        if (mouthParticulate == nullptr)
        {
          double aerosolConcentration_g_Per_m3 = 2.5;
          pc.GetSubstances().AddActiveSubstance(substance);
          ambient->GetSubstanceQuantity(substance)->GetConcentration().SetValue(aerosolConcentration_g_Per_m3, MassPerVolumeUnit::g_Per_m3);
          ambient->Balance(BalanceLiquidBy::Concentration);
  
          mouthParticulate = aGraph->GetCompartment(pulse::PulmonaryCompartment::Mouth)->GetSubstanceQuantity(substance);
          carinaParticulate = aGraph->GetCompartment(pulse::PulmonaryCompartment::Carina)->GetSubstanceQuantity(substance);
          leftDeadSpaceParticulate = aGraph->GetCompartment(pulse::PulmonaryCompartment::LeftDeadSpace)->GetSubstanceQuantity(substance);
          leftAlveoliParticulate = aGraph->GetCompartment(pulse::PulmonaryCompartment::LeftAlveoli)->GetSubstanceQuantity(substance);
          rightDeadSpaceParticulate = aGraph->GetCompartment(pulse::PulmonaryCompartment::RightDeadSpace)->GetSubstanceQuantity(substance);
          rightAlveoliParticulate = aGraph->GetCompartment(pulse::PulmonaryCompartment::RightAlveoli)->GetSubstanceQuantity(substance);
  
          // Initialize to the env concentrations
          // aerosolConcentration_g_Per_m3 = 0.; // Or initialize to zero if you want.
          mouthParticulate->GetConcentration().SetValue(aerosolConcentration_g_Per_m3, MassPerVolumeUnit::g_Per_m3); mouthParticulate->Balance(BalanceLiquidBy::Concentration);
          carinaParticulate->GetConcentration().SetValue(aerosolConcentration_g_Per_m3, MassPerVolumeUnit::g_Per_m3); carinaParticulate->Balance(BalanceLiquidBy::Concentration);
          leftDeadSpaceParticulate->GetConcentration().SetValue(aerosolConcentration_g_Per_m3, MassPerVolumeUnit::g_Per_m3); leftDeadSpaceParticulate->Balance(BalanceLiquidBy::Concentration);
          leftAlveoliParticulate->GetConcentration().SetValue(aerosolConcentration_g_Per_m3, MassPerVolumeUnit::g_Per_m3); leftAlveoliParticulate->Balance(BalanceLiquidBy::Concentration);
          rightDeadSpaceParticulate->GetConcentration().SetValue(aerosolConcentration_g_Per_m3, MassPerVolumeUnit::g_Per_m3); rightDeadSpaceParticulate->Balance(BalanceLiquidBy::Concentration);
          rightAlveoliParticulate->GetConcentration().SetValue(aerosolConcentration_g_Per_m3, MassPerVolumeUnit::g_Per_m3); rightAlveoliParticulate->Balance(BalanceLiquidBy::Concentration);
        }
  
        if (Env2Mouth->GetFlow().IsPositive())
        {
          // Grab instantaneous        
          double inspiredAir_mL = Env2Mouth->GetFlow(VolumePerTimeUnit::mL_Per_s) * deltaT_s;
          double inspiredParticulate_ug = Env2Mouth->GetFlow(VolumePerTimeUnit::mL_Per_s)*ambient->GetSubstanceQuantity(substance)->GetConcentration(MassPerVolumeUnit::ug_Per_mL) * deltaT_s;
          // Totals from instantaneous
          totalInspiredAir_mL += inspiredAir_mL;
          totalInspiredParticulate_ug += inspiredParticulate_ug;
        }
        else
        {
          double exspiredAir_mL = Env2Mouth->GetFlow(VolumePerTimeUnit::mL_Per_s) * deltaT_s;
          totalExspiredAir_mL += exspiredAir_mL;
        }
        if(true)
        {         
          // Calculate the deposited mass by concentration and flow
          depositedMouthParticulate_ug = mouthParticulate->GetConcentration(MassPerVolumeUnit::ug_Per_mL)*aGraph->GetCompartment(pulse::PulmonaryCompartment::Mouth)->GetInFlow(VolumePerTimeUnit::mL_Per_s)*deltaT_s*SIDECoeff.GetMouth();
          depositedCarinaParticulate_ug = carinaParticulate->GetConcentration(MassPerVolumeUnit::ug_Per_mL)*aGraph->GetCompartment(pulse::PulmonaryCompartment::Carina)->GetInFlow(VolumePerTimeUnit::mL_Per_s)*deltaT_s*SIDECoeff.GetCarina();
          depositedLeftDeadSpaceParticulate_ug = leftDeadSpaceParticulate->GetConcentration(MassPerVolumeUnit::ug_Per_mL)*aGraph->GetCompartment(pulse::PulmonaryCompartment::LeftDeadSpace)->GetInFlow(VolumePerTimeUnit::mL_Per_s)*deltaT_s*SIDECoeff.GetDeadSpace();
          depositedLeftAlveoliParticulate_ug = leftAlveoliParticulate->GetConcentration(MassPerVolumeUnit::ug_Per_mL)*aGraph->GetCompartment(pulse::PulmonaryCompartment::LeftAlveoli)->GetInFlow(VolumePerTimeUnit::mL_Per_s)*deltaT_s*SIDECoeff.GetAlveoli();
          depositedRightDeadSpaceParticulate_ug = rightDeadSpaceParticulate->GetConcentration(MassPerVolumeUnit::ug_Per_mL)*aGraph->GetCompartment(pulse::PulmonaryCompartment::RightDeadSpace)->GetInFlow(VolumePerTimeUnit::mL_Per_s)*deltaT_s*SIDECoeff.GetDeadSpace();
          depositedRightAlveoliParticulate_ug = rightAlveoliParticulate->GetConcentration(MassPerVolumeUnit::ug_Per_mL)*aGraph->GetCompartment(pulse::PulmonaryCompartment::RightAlveoli)->GetInFlow(VolumePerTimeUnit::mL_Per_s)*deltaT_s*SIDECoeff.GetAlveoli();
          
          // Total deposited
          totalParticulateDeposited_ug += depositedMouthParticulate_ug + depositedCarinaParticulate_ug +
            depositedLeftDeadSpaceParticulate_ug + depositedLeftAlveoliParticulate_ug +
            depositedRightDeadSpaceParticulate_ug + depositedRightAlveoliParticulate_ug;
  
          // Remove the delta from the compartment substance quantity mass
          mouthParticulate->GetMass().IncrementValue(-depositedMouthParticulate_ug, MassUnit::ug); mouthParticulate->Balance(BalanceLiquidBy::Mass);
          carinaParticulate->GetMass().IncrementValue(-depositedCarinaParticulate_ug, MassUnit::ug); carinaParticulate->Balance(BalanceLiquidBy::Mass);
          leftDeadSpaceParticulate->GetMass().IncrementValue(-depositedLeftDeadSpaceParticulate_ug, MassUnit::ug); leftDeadSpaceParticulate->Balance(BalanceLiquidBy::Mass);
          leftAlveoliParticulate->GetMass().IncrementValue(-depositedLeftAlveoliParticulate_ug, MassUnit::ug); leftAlveoliParticulate->Balance(BalanceLiquidBy::Mass);
          rightDeadSpaceParticulate->GetMass().IncrementValue(-depositedRightDeadSpaceParticulate_ug, MassUnit::ug); rightDeadSpaceParticulate->Balance(BalanceLiquidBy::Mass);
          rightAlveoliParticulate->GetMass().IncrementValue(-depositedRightAlveoliParticulate_ug, MassUnit::ug); rightAlveoliParticulate->Balance(BalanceLiquidBy::Mass);
          // Add the delta to the compartment substance quantity mass deposited
          mouthParticulate->GetMassDeposited().IncrementValue(depositedMouthParticulate_ug, MassUnit::ug);
          carinaParticulate->GetMassDeposited().IncrementValue(depositedCarinaParticulate_ug, MassUnit::ug);
          leftDeadSpaceParticulate->GetMassDeposited().IncrementValue(depositedLeftDeadSpaceParticulate_ug, MassUnit::ug);
          leftAlveoliParticulate->GetMassDeposited().IncrementValue(depositedLeftAlveoliParticulate_ug, MassUnit::ug);
          rightDeadSpaceParticulate->GetMassDeposited().IncrementValue(depositedRightDeadSpaceParticulate_ug, MassUnit::ug);
          rightAlveoliParticulate->GetMassDeposited().IncrementValue(depositedRightAlveoliParticulate_ug, MassUnit::ug);        
        }
      }
  
      trk.Track("MouthConcentration_ug_Per_mL", time, mouthParticulate == nullptr ? 0 : mouthParticulate->GetConcentration(MassPerVolumeUnit::ug_Per_mL));
      trk.Track("CarinaConcentration_ug_Per_mL", time, carinaParticulate == nullptr ? 0 : carinaParticulate->GetConcentration(MassPerVolumeUnit::ug_Per_mL));
  
      trk.Track("MouthInFlow_mL_Per_s", time, aGraph->GetCompartment(pulse::PulmonaryCompartment::Mouth)->GetInFlow(VolumePerTimeUnit::mL_Per_s));
      trk.Track("CarinaInFlow_mL_Per_s", time, aGraph->GetCompartment(pulse::PulmonaryCompartment::Carina)->GetInFlow(VolumePerTimeUnit::mL_Per_s));
  
      trk.Track("TotalInspiredAir_mL", time, totalInspiredAir_mL);
      trk.Track("TotalExspiredAir_mL", time, totalExspiredAir_mL);
      trk.Track("TotalInspiredParticulate_ug", time, totalInspiredParticulate_ug);
      trk.Track("TotalDepositedParticulate_ug", time, totalParticulateDeposited_ug);
      trk.Track("TotalMouthPostParticulate_ug", time, mouthParticulate==nullptr?0:mouthParticulate->GetMass(MassUnit::ug));
      trk.Track("MouthParticulateDeposited_ug", time, depositedMouthParticulate_ug);
      trk.Track("TotalMouthParticulateDeposited_ug", time, mouthParticulate == nullptr ? 0 : mouthParticulate->GetMassDeposited(MassUnit::ug));
      trk.Track("CarinaParticulate_ug", time, carinaParticulate == nullptr ? 0 : carinaParticulate->GetMass(MassUnit::ug));
      trk.Track("CarinaParticulateDeposited_ug", time, carinaParticulate == nullptr ? 0 : carinaParticulate->GetMassDeposited(MassUnit::ug));
      trk.Track("LeftDeadSpaceParticulate_ug", time, leftDeadSpaceParticulate == nullptr ? 0 : leftDeadSpaceParticulate->GetMass(MassUnit::ug));
      trk.Track("LeftDeadSpaceParticulateDeposited_ug", time, leftDeadSpaceParticulate == nullptr ? 0 : leftDeadSpaceParticulate->GetMassDeposited(MassUnit::ug));
      trk.Track("LeftAlveoliParticulate_ug", time, leftAlveoliParticulate == nullptr ? 0 : leftAlveoliParticulate->GetMass(MassUnit::ug));
      trk.Track("LeftAlveoliParticulateDeposited_ug", time, leftAlveoliParticulate == nullptr ? 0 : leftAlveoliParticulate->GetMassDeposited(MassUnit::ug));
      trk.Track("RightDeadSpaceParticulate_ug", time, rightDeadSpaceParticulate == nullptr ? 0 : rightDeadSpaceParticulate->GetMass(MassUnit::ug));
      trk.Track("RightDeadSpaceParticulateDeposited_ug", time, rightDeadSpaceParticulate == nullptr ? 0 : rightDeadSpaceParticulate->GetMassDeposited(MassUnit::ug));
      trk.Track("RightAlveoliParticulate_ug", time, rightAlveoliParticulate == nullptr ? 0 : rightAlveoliParticulate->GetMass(MassUnit::ug));
      trk.Track("RightAlveoliParticulateDeposited_ug", time, rightAlveoliParticulate == nullptr ? 0 : rightAlveoliParticulate->GetMassDeposited(MassUnit::ug));
  
      if (i == 0)
        trk.CreateFile(std::string("./test_results/unit_tests/Pulse/"+substance.GetName()+"DepositionFraction.txt").c_str(), file);
      trk.StreamTrackToFile(file);
      
      time += deltaT_s;
    }
  
    // Do all the deltas match?
    double expectedInspiredParticulate_ug = totalInspiredAir_mL * ambient->GetSubstanceQuantity(substance)->GetConcentration(MassPerVolumeUnit::ug_Per_mL);  
    if (GeneralMath::PercentTolerance(totalInspiredParticulate_ug, expectedInspiredParticulate_ug) > PercentTolerance)
    {
      m_ss << " Compartment particulate mass does not expected particulate mass : " << totalInspiredParticulate_ug << " expected " << expectedInspiredParticulate_ug;
      tc.AddFailure(m_ss);
    }
  
    m_ss << "Total Particulate Inhaled : " << totalInspiredParticulate_ug;
    Info(m_ss);
    m_ss << "Total Particulate Deposited : " << totalParticulateDeposited_ug;
    Info(m_ss);
  
    m_ss << "Particulate Deposited in Mouth: " << mouthParticulate->GetMassDeposited(MassUnit::ug);
    Info(m_ss);
    if (GeneralMath::PercentTolerance(expectedMouthDepFrac, mouthParticulate->GetMassDeposited(MassUnit::ug) / totalInspiredParticulate_ug) > PercentTolerance)
    {
      m_ss << " Mouth Particle Deposition is not correct : " << mouthParticulate->GetMassDeposited(MassUnit::ug) / totalInspiredParticulate_ug << " expected " << expectedMouthDepFrac;
      tc.AddFailure(m_ss);
    }
  
    m_ss << "Particulate Deposited in Carina: " << carinaParticulate->GetMassDeposited(MassUnit::ug);
    Info(m_ss);
    if (GeneralMath::PercentTolerance(expectedCarinaDepFrac, carinaParticulate->GetMassDeposited(MassUnit::ug) / totalInspiredParticulate_ug) > PercentTolerance)
    {
      m_ss << " Carina Particle Deposition is not correct : " << carinaParticulate->GetMassDeposited(MassUnit::ug) / totalInspiredParticulate_ug << " expected " << expectedCarinaDepFrac;
      tc.AddFailure(m_ss);
    }
  
    m_ss << "Particulate Deposited in Left Dead Space : " << leftDeadSpaceParticulate->GetMassDeposited(MassUnit::ug);
    Info(m_ss);
    if (GeneralMath::PercentTolerance(expectedDeadSpaceDepFrac, leftDeadSpaceParticulate->GetMassDeposited(MassUnit::ug) / totalInspiredParticulate_ug) > PercentTolerance)
    {
      m_ss << " Left Dead Space Particle Deposition is not correct : " << leftDeadSpaceParticulate->GetMassDeposited(MassUnit::ug) / totalInspiredParticulate_ug << " expected " << expectedDeadSpaceDepFrac;
      tc.AddFailure(m_ss);
    }
  
    m_ss << "Particulate Deposited in Left Alveoli : " << leftAlveoliParticulate->GetMassDeposited(MassUnit::ug);
    Info(m_ss);
    if (GeneralMath::PercentTolerance(expectedAlveoliDepFrac, leftAlveoliParticulate->GetMassDeposited(MassUnit::ug) / totalInspiredParticulate_ug) > PercentTolerance)
    {
      m_ss << " Left Alveoli Particle Deposition is not correct : " << leftAlveoliParticulate->GetMassDeposited(MassUnit::ug) / totalInspiredParticulate_ug << " expected " << expectedAlveoliDepFrac;
      tc.AddFailure(m_ss);
    }
  
    m_ss << "Particulate Deposited in Right Dead Space : " << rightDeadSpaceParticulate->GetMassDeposited(MassUnit::ug);
    Info(m_ss);
    if (GeneralMath::PercentTolerance(expectedDeadSpaceDepFrac, rightDeadSpaceParticulate->GetMassDeposited(MassUnit::ug) / totalInspiredParticulate_ug) > PercentTolerance)
    {
      m_ss << " Right Dead Space Particle Deposition is not correct : " << rightDeadSpaceParticulate->GetMassDeposited(MassUnit::ug) / totalInspiredParticulate_ug << " expected " << expectedDeadSpaceDepFrac;
      tc.AddFailure(m_ss);
    }
  
    m_ss << "Particulate Deposited in Right Alveoli : " << rightAlveoliParticulate->GetMassDeposited(MassUnit::ug);
    Info(m_ss);
    if (GeneralMath::PercentTolerance(expectedAlveoliDepFrac, rightAlveoliParticulate->GetMassDeposited(MassUnit::ug) / totalInspiredParticulate_ug) > PercentTolerance)
    {
      m_ss << " Right Alveoli Particle Deposition is not correct : " << rightAlveoliParticulate->GetMassDeposited(MassUnit::ug) / totalInspiredParticulate_ug << " expected " << expectedAlveoliDepFrac;
      tc.AddFailure(m_ss);
    }
  
    tc.GetDuration().SetValue(pTimer.GetElapsedTime_s("Test"), TimeUnit::s);
  }