/* Distributed under the Apache License, Version 2.0.
   See accompanying NOTICE file for details.*/
#include "EngineTest.h"
#include "CommonDataModel.h"
#include "controller/Controller.h"
#include "controller/Substances.h"
#include "controller/Circuits.h"
#include "controller/Compartments.h"
#include "physiology/Tissue.h"
#include "physiology/Saturation.h"
#include "PulseConfiguration.h"

#include "patient/SEPatient.h"
#include "substance/SESubstance.h"
#include "substance/SESubstanceManager.h"
#include "circuit/fluid/SEFluidCircuit.h"
#include "circuit/fluid/SEFluidCircuitCalculator.h"
#include "compartment/SECompartmentManager.h"
#include "compartment/fluid/SEFluidCompartmentLink.h"
#include "compartment/fluid/SELiquidCompartmentGraph.h"
#include "compartment/fluid/SELiquidCompartment.h"
#include "utils/GeneralMath.h"
#include "utils/DataTrack.h"
#include "properties/SEScalar0To1.h"
#include "properties/SEScalarMass.h"
#include "properties/SEScalarMassPerVolume.h"
#include "properties/SEScalarMassPerAmount.h"
#include "properties/SEScalarVolume.h"
#include "properties/SEScalarAmountPerVolume.h"

double TotalHbMols(SELiquidCompartmentGraph& Graph, SESubstance& Hb, SESubstance& HbO2, SESubstance& HbO2CO2, SESubstance& HbCO2)
{
  double Hb_g_Per_mol = Hb.GetMolarMass(MassPerAmountUnit::g_Per_mol);
  double HbO2_g_Per_mol = HbO2.GetMolarMass(MassPerAmountUnit::g_Per_mol);
  double HbCO2_g_Per_mol = HbCO2.GetMolarMass(MassPerAmountUnit::g_Per_mol);
  double HbO2CO2_g_Per_mol = HbO2CO2.GetMolarMass(MassPerAmountUnit::g_Per_mol);

  double totalHb_g = 0.0;
  double totalHbO2_g = 0.0;
  double totalHbCO2_g = 0.0;
  double totalHBO2CO2_g = 0.0;

  for (SELiquidCompartment* cmpt : Graph.GetCompartments())
  {
    totalHb_g += cmpt->GetSubstanceQuantity(Hb)->GetMass(MassUnit::g);
    totalHbO2_g += cmpt->GetSubstanceQuantity(HbO2)->GetMass(MassUnit::g);
    totalHbCO2_g += cmpt->GetSubstanceQuantity(HbCO2)->GetMass(MassUnit::g);
    totalHBO2CO2_g += cmpt->GetSubstanceQuantity(HbO2CO2)->GetMass(MassUnit::g);
  }

  return totalHb_g / Hb_g_Per_mol + totalHbO2_g / HbO2_g_Per_mol + totalHbCO2_g / HbCO2_g_Per_mol + totalHBO2CO2_g / HbO2CO2_g_Per_mol;
}

void PulseEngineTest::FourCompartmentTest(bool usingAcidBase, bool usingProductionConsumption, bool usingDiffusion, const std::string& rptDirectory)
{
  DataTrack trk;
  std::string outputName;
  if (!usingAcidBase && !usingProductionConsumption && !usingDiffusion)
  {
    outputName = "/FourCompartmentTestSimple";
  }
  else if (usingAcidBase && !usingProductionConsumption && !usingDiffusion)
  {
    outputName = "/AcidBaseFourCompartmentTest";
  }
  else if (!usingAcidBase && !usingProductionConsumption && usingDiffusion)
  {
    outputName = "/FiveCompartmentTestWithDiffusion";
  }
  else if (usingAcidBase && usingProductionConsumption && !usingDiffusion)
  {
    outputName = "/AcidBaseFourCompartmentTestWithProductionConsumption";
  }
  else if (usingAcidBase && !usingProductionConsumption && usingDiffusion)
  {
    outputName = "/AcidBaseFiveCompartmentTestWithDiffusion";
  }
  else if (usingAcidBase && usingProductionConsumption && usingDiffusion)
  {
    outputName = "/AcidBaseFiveCompartmentTestWithProductionConsumptionAndDiffusion";
  }
  m_Logger->ResetLogFile(rptDirectory + outputName + ".log");
  std::ofstream file;
  SELiquidTransporter txpt(VolumePerTimeUnit::mL_Per_s, VolumeUnit::mL, MassUnit::ug, MassPerVolumeUnit::ug_Per_mL, m_Logger);
  SEFluidCircuitCalculator calc(VolumePerPressureUnit::mL_Per_mmHg, VolumePerTimeUnit::mL_Per_s, PressureTimeSquaredPerVolumeUnit::mmHg_s2_Per_mL, PressureUnit::mmHg, VolumeUnit::mL, PressureTimePerVolumeUnit::mmHg_s_Per_mL, m_Logger);
  PulseController pc(m_Logger);
  Tissue& tsu = (Tissue&)pc.GetTissue();
  SEPatient patient(pc.GetLogger());
  patient.SerializeFromFile("./patients/StandardMale.json", JSON);
  pc.SetupPatient(patient);
  pc.m_Config->EnableRenal(eSwitch::Off);
  pc.m_Config->EnableTissue(eSwitch::Off);
  pc.CreateCircuitsAndCompartments();
  pc.GetSubstances().InitializeGasCompartments();


  //Create a circuit
  double veinsPressure_mmHg = 4.0; //Central Venous Pressure
  double pulmonaryPressure_mmHg = 10.0; //Average Pressure
  double arteriesPressure_mmHg = 90.0; //Mean Arterial Pressure
  double capillariesPressure_mmHg = 30.0; //http://www.cvphysiology.com/Blood%20Pressure/BP019.htm
  double flow_mL_Per_min = 5600.0;

  //http://droualb.faculty.mjc.edu/Course%20Materials/Physiology%20101/Chapter%20Notes/Fall%202011/chapter_14%20Fall%202011.htm
  double totalVolume_L = 5.0;
  double arteriesVolume = 0.23 * totalVolume_L;
  double veinsVolume_L = 0.60 * totalVolume_L;
  double pulmonaryVolume = 0.12 * totalVolume_L;
  double capillariesVolume = 0.05 * totalVolume_L;

  SECircuitManager circuits(m_Logger);
  SEFluidCircuit* Circuit = &circuits.CreateFluidCircuit("Circuit");

  SEFluidCircuitNode& nPulmonary = Circuit->CreateNode("Pulmonary");
  nPulmonary.GetVolumeBaseline().SetValue(pulmonaryVolume, VolumeUnit::L);
  SEFluidCircuitNode& nVeins = Circuit->CreateNode("Veins");
  nVeins.GetVolumeBaseline().SetValue(veinsVolume_L, VolumeUnit::L);
  SEFluidCircuitNode& nCapillaries = Circuit->CreateNode("Capillaries");
  nCapillaries.GetVolumeBaseline().SetValue(capillariesVolume, VolumeUnit::L);
  SEFluidCircuitNode& nArteries = Circuit->CreateNode("Arteries");
  nArteries.GetVolumeBaseline().SetValue(arteriesVolume, VolumeUnit::L);
  SEFluidCircuitNode& nGround = Circuit->CreateNode("Ground");  
  nGround.SetAsReferenceNode();
  nGround.GetPressure().SetValue(0.0, PressureUnit::mmHg);

  SEFluidCircuitPath& pPulmonaryToArteries = Circuit->CreatePath(nPulmonary, nArteries, "PulmonaryToArteries");
  pPulmonaryToArteries.GetFlowSourceBaseline().SetValue(flow_mL_Per_min, VolumePerTimeUnit::mL_Per_min);
  SEFluidCircuitPath& pArteriesToCapillaries = Circuit->CreatePath(nArteries, nCapillaries, "ArteriesToCapillaries");
  pArteriesToCapillaries.GetFlowSourceBaseline().SetValue(flow_mL_Per_min, VolumePerTimeUnit::mL_Per_min);
  SEFluidCircuitPath& pCapillariesToVeins = Circuit->CreatePath(nCapillaries, nVeins, "CapillariesToVeins");
  pCapillariesToVeins.GetFlowSourceBaseline().SetValue(flow_mL_Per_min, VolumePerTimeUnit::mL_Per_min);
  SEFluidCircuitPath& pVeinsToPulmonary = Circuit->CreatePath(nVeins, nPulmonary, "VeinsToPulmonary");
  pVeinsToPulmonary.GetFlowSourceBaseline().SetValue(flow_mL_Per_min, VolumePerTimeUnit::mL_Per_min);

  SEFluidCircuitPath& pGroundToArteries = Circuit->CreatePath(nGround, nArteries, "GroundToArteries");
  pGroundToArteries.GetPressureSourceBaseline().SetValue(arteriesPressure_mmHg, PressureUnit::mmHg);
  SEFluidCircuitPath& pGroundToPulmonary = Circuit->CreatePath(nGround, nPulmonary, "GroundToPulmonary");
  pGroundToPulmonary.GetPressureSourceBaseline().SetValue(pulmonaryPressure_mmHg, PressureUnit::mmHg);
  SEFluidCircuitPath& pGroundToCapillaries = Circuit->CreatePath(nGround, nCapillaries, "GroundToCapillaries");
  pGroundToCapillaries.GetPressureSourceBaseline().SetValue(capillariesPressure_mmHg, PressureUnit::mmHg);
  SEFluidCircuitPath& pGroundToVeins = Circuit->CreatePath(nGround, nVeins, "GroundToVeins");
  pGroundToVeins.GetPressureSourceBaseline().SetValue(veinsPressure_mmHg, PressureUnit::mmHg);

  Circuit->SetNextAndCurrentFromBaselines();
  Circuit->StateChange();

  SELiquidCompartment& cPulmonary = pc.GetCompartments().CreateLiquidCompartment("Pulmonary");
  cPulmonary.MapNode(nPulmonary);
  SELiquidCompartment& cVeins = pc.GetCompartments().CreateLiquidCompartment("Veins");
  cVeins.MapNode(nVeins);
  SELiquidCompartment& cArteries = pc.GetCompartments().CreateLiquidCompartment("Arteries");
  cArteries.MapNode(nArteries);
  SELiquidCompartment& cCapillaries = pc.GetCompartments().CreateLiquidCompartment("Capillaries");
  cCapillaries.MapNode(nCapillaries);

  SELiquidCompartment& cTissue = pc.GetCompartments().CreateLiquidCompartment("Tissue");

  //Create the links
  SELiquidCompartmentLink& lPulmonaryToArteries = pc.GetCompartments().CreateLiquidLink(cPulmonary, cArteries, "PulmonaryToArteries");
  lPulmonaryToArteries.MapPath(pPulmonaryToArteries);
  SELiquidCompartmentLink& lArteriesToCapillaries = pc.GetCompartments().CreateLiquidLink(cArteries, cCapillaries, "ArteriesToCapillaries");
  lArteriesToCapillaries.MapPath(pArteriesToCapillaries);
  SELiquidCompartmentLink& lCapillariesToVeins = pc.GetCompartments().CreateLiquidLink(cCapillaries, cVeins, "CapillariesToVeins");
  lCapillariesToVeins.MapPath(pCapillariesToVeins);
  SELiquidCompartmentLink& lVeinsToPulmonary = pc.GetCompartments().CreateLiquidLink(cVeins, cPulmonary, "VeinsToPulmonary");
  lVeinsToPulmonary.MapPath(pVeinsToPulmonary);

  //Create the graph
  SELiquidCompartmentGraph& Graph = pc.GetCompartments().CreateLiquidGraph("Graph");
  Graph.AddCompartment(cPulmonary);
  Graph.AddCompartment(cArteries);
  Graph.AddCompartment(cVeins);
  Graph.AddCompartment(cCapillaries);
  Graph.AddLink(lPulmonaryToArteries);
  Graph.AddLink(lArteriesToCapillaries);
  Graph.AddLink(lCapillariesToVeins);
  Graph.AddLink(lVeinsToPulmonary);
  Graph.StateChange();


  //Initialize substances
  SESubstance& O2 = pc.GetSubstances().GetO2();
  SESubstance& Hb = pc.GetSubstances().GetHb();
  SESubstance& HbO2 = pc.GetSubstances().GetHbO2();
  SESubstance& HbO2CO2 = pc.GetSubstances().GetHbO2CO2();
  SESubstance& CO2 = pc.GetSubstances().GetCO2();
  SESubstance& HCO3 = pc.GetSubstances().GetHCO3();
  SESubstance& HbCO2 = pc.GetSubstances().GetHbCO2();
  SESubstance& N2 = pc.GetSubstances().GetN2();

  pc.GetSubstances().AddActiveSubstance(O2);
  pc.GetSubstances().AddActiveSubstance(Hb);
  pc.GetSubstances().AddActiveSubstance(HbO2);
  pc.GetSubstances().AddActiveSubstance(HbO2CO2);
  pc.GetSubstances().AddActiveSubstance(CO2);
  pc.GetSubstances().AddActiveSubstance(HCO3);
  pc.GetSubstances().AddActiveSubstance(HbCO2);
  pc.GetSubstances().AddActiveSubstance(N2);

  SEScalarMassPerVolume   albuminConcentration;
  SEScalar0To1        hematocrit;
  SEScalarTemperature     bodyTemp;
  SEScalarAmountPerVolume strongIonDifference;
  SEScalarAmountPerVolume phosphate;

  albuminConcentration.SetValue(45.0, MassPerVolumeUnit::g_Per_L);
  hematocrit.SetValue(0.45);
  bodyTemp.SetValue(37.0, TemperatureUnit::C);
  strongIonDifference.SetValue(40.5, AmountPerVolumeUnit::mmol_Per_L);
  phosphate.SetValue(1.1, AmountPerVolumeUnit::mmol_Per_L);

  double Hb_total_g_Per_dL = hematocrit.GetValue()*34.0;
  double Hb_total_mM = Hb_total_g_Per_dL / Hb.GetMolarMass(MassPerAmountUnit::g_Per_mmol) * 10.0;

  
  pc.GetSaturationCalculator().SetBodyState(albuminConcentration, hematocrit, bodyTemp, strongIonDifference, phosphate);

  // Highly Oxygenated (Aorta/PulmCaps) (90.5/40 O2/CO2 PP)
  double High_CO2_sat = 0.0282123;
  double High_O2_sat = 0.974759;
  double High_CO2_mmol_Per_L = 1.30049;
  double High_HCO3_mmol_Per_L = 25.9377;
  double High_O2_mmol_Per_L = 0.129065;
  double High_pH = 7.39982;

  // Low Oxygenated (Everywhere else) (36.5/45.5 O2/CO2 PP)
  double Low_CO2_sat = 0.162382;
  double Low_O2_sat = 0.789701;
  double Low_CO2_mmol_Per_L = 1.40557;
  double Low_HCO3_mmol_Per_L = 26.1182;
  double Low_O2_mmol_Per_L = 0.0560892;
  double Low_pH = 7.36909;

  pc.GetSubstances().InitializeBloodGases(cPulmonary, Hb_total_mM, High_O2_sat, High_O2_mmol_Per_L, High_CO2_sat, High_CO2_mmol_Per_L, High_HCO3_mmol_Per_L, High_pH, false);
  pc.GetSubstances().InitializeBloodGases(cArteries, Hb_total_mM, High_O2_sat, High_O2_mmol_Per_L, High_CO2_sat, High_CO2_mmol_Per_L, High_HCO3_mmol_Per_L, High_pH, false);
  pc.GetSubstances().InitializeBloodGases(cCapillaries, Hb_total_mM, Low_O2_sat, Low_O2_mmol_Per_L, Low_CO2_sat, Low_CO2_mmol_Per_L, Low_HCO3_mmol_Per_L, Low_pH, false);
  pc.GetSubstances().InitializeBloodGases(cVeins, Hb_total_mM, Low_O2_sat, Low_O2_mmol_Per_L, Low_CO2_sat, Low_CO2_mmol_Per_L, Low_HCO3_mmol_Per_L, Low_pH, false);

  if (usingAcidBase)
  {
    pc.GetSaturationCalculator().CalculateBloodGasDistribution(cPulmonary);
    pc.GetSaturationCalculator().CalculateBloodGasDistribution(cArteries);
    pc.GetSaturationCalculator().CalculateBloodGasDistribution(cCapillaries);
    pc.GetSaturationCalculator().CalculateBloodGasDistribution(cVeins);
  }

  //Tissue diffusion values
  double tissueVolume_L = 61.7;  //sum of all tissue volumes
  double tissueO2Mass_g = tissueVolume_L * cCapillaries.GetSubstanceQuantity(O2)->GetConcentration(MassPerVolumeUnit::g_Per_L);
  double tissueCO2Mass_g = tissueVolume_L * cCapillaries.GetSubstanceQuantity(CO2)->GetConcentration(MassPerVolumeUnit::g_Per_L);
  if (usingDiffusion)
  {
    //cTissue.GetTotalVolume().SetValue(tissueVolume_L, VolumeUnit::L);
    cTissue.GetVolume().SetValue(tissueVolume_L, VolumeUnit::L);
    cTissue.GetSubstanceQuantity(O2)->GetMass().SetValue(tissueO2Mass_g, MassUnit::g);
    cTissue.GetSubstanceQuantity(Hb)->GetMass().SetValue(0, MassUnit::g);
    cTissue.GetSubstanceQuantity(HbO2)->GetMass().SetValue(0, MassUnit::g);
    cTissue.GetSubstanceQuantity(HbO2CO2)->GetMass().SetValue(0, MassUnit::g);
    cTissue.GetSubstanceQuantity(CO2)->GetMass().SetValue(tissueCO2Mass_g, MassUnit::g);
    cTissue.GetSubstanceQuantity(HCO3)->GetMass().SetValue(0, MassUnit::g);
    cTissue.GetSubstanceQuantity(HbCO2)->GetMass().SetValue(0, MassUnit::g);
    cTissue.Balance(BalanceLiquidBy::Mass);
  }
  

  //Production, Consumption / Alveolar Transfer values
  double O2Consumption_ugPerS = 5954.2;
  double CO2Production_ugPerS = 6140.4;

  //Execution parameters
  double time = 0;
  double deltaT_s = 1.0 / 90.0;
  double runTime_min = 3.0;
  std::stringstream ss;

  for (unsigned int i = 0; i < runTime_min * 60.0 / deltaT_s; i++)
  {
    //PreProcess - Mimic Produce, Consume, and Alveolar Transfer
    if (usingProductionConsumption)
    {
      if (usingDiffusion)
      {
        //Do it in the tissues and let it diffuse
        cPulmonary.GetSubstanceQuantity(O2)->GetMass().IncrementValue(O2Consumption_ugPerS * deltaT_s, MassUnit::ug);
        cPulmonary.GetSubstanceQuantity(CO2)->GetMass().IncrementValue(-CO2Production_ugPerS * deltaT_s, MassUnit::ug);
        cTissue.GetSubstanceQuantity(O2)->GetMass().IncrementValue(-O2Consumption_ugPerS * deltaT_s, MassUnit::ug);
        cTissue.GetSubstanceQuantity(CO2)->GetMass().IncrementValue(CO2Production_ugPerS * deltaT_s, MassUnit::ug);

        //See if we took too much out
        double pulmonaryCO2Mass_ug = cPulmonary.GetSubstanceQuantity(CO2)->GetMass(MassUnit::ug);
        double tissueO2Mass_ug = cCapillaries.GetSubstanceQuantity(O2)->GetMass(MassUnit::ug);
        if (pulmonaryCO2Mass_ug < 0.0)
        {
          ss << "Not Enough Pulmonary CO2 at time " << i*deltaT_s;
          cPulmonary.GetSubstanceQuantity(CO2)->GetMass().SetValue(0.0, MassUnit::ug);
          Info(ss.str());
          ss.str("");
          ss.clear();
        }
        if (tissueO2Mass_ug < 0.0)
        {
          ss << "Not Enough Tissue O2 at time " << i*deltaT_s;
          cTissue.GetSubstanceQuantity(O2)->GetMass().SetValue(0.0, MassUnit::ug);
          Info(ss.str());
          ss.str("");
          ss.clear();
        }

        cPulmonary.Balance(BalanceLiquidBy::Mass);
        cTissue.Balance(BalanceLiquidBy::Mass);
      }
      else
      {
        //No tissues, so put straight into the blood
        cPulmonary.GetSubstanceQuantity(O2)->GetMass().IncrementValue(O2Consumption_ugPerS * deltaT_s, MassUnit::ug);
        cPulmonary.GetSubstanceQuantity(CO2)->GetMass().IncrementValue(-CO2Production_ugPerS * deltaT_s, MassUnit::ug);
        cCapillaries.GetSubstanceQuantity(O2)->GetMass().IncrementValue(-O2Consumption_ugPerS * deltaT_s, MassUnit::ug);
        cCapillaries.GetSubstanceQuantity(CO2)->GetMass().IncrementValue(CO2Production_ugPerS * deltaT_s, MassUnit::ug);

        //See if we took too much out
        double pulmonaryCO2Mass_ug = cPulmonary.GetSubstanceQuantity(CO2)->GetMass(MassUnit::ug);
        double capillaryO2Mass_ug = cCapillaries.GetSubstanceQuantity(O2)->GetMass(MassUnit::ug);
        if (pulmonaryCO2Mass_ug < 0.0)
        {
          ss << "Not Enough Pulmonary CO2 at time " << i*deltaT_s;
          cPulmonary.GetSubstanceQuantity(CO2)->GetMass().SetValue(0.0, MassUnit::ug);
          Info(ss.str());
          ss.str("");
          ss.clear();
        }
        if (capillaryO2Mass_ug < 0.0)
        {
          ss << "Not Enough Capillary O2 at time " << i*deltaT_s;
          cCapillaries.GetSubstanceQuantity(O2)->GetMass().SetValue(0.0, MassUnit::ug);
          Info(ss.str());
          ss.str("");
          ss.clear();
        }

        cPulmonary.Balance(BalanceLiquidBy::Mass);
        cCapillaries.Balance(BalanceLiquidBy::Mass);
      }
    }

    //Process - Execute the circuit
    calc.Process(*Circuit, deltaT_s);
    //Execute the substance transport function
    txpt.Transport(Graph, deltaT_s);

  if (usingAcidBase)
  {
    for (SELiquidCompartment* cmpt : Graph.GetCompartments())
    {
      pc.GetSaturationCalculator().CalculateBloodGasDistribution(*cmpt);
    }
  }

    //Diffuse if necessary
    if (usingDiffusion)
    {
      tsu.MoveMassByInstantDiffusion(cCapillaries, cTissue, O2, deltaT_s);
      tsu.MoveMassByInstantDiffusion(cCapillaries, cTissue, CO2, deltaT_s);
    }

    //convert 'Next' values to current
    calc.PostProcess(*Circuit);

    trk.Track(time, *Circuit);
    trk.Track(time, Graph);

    //Total masses
    for (SESubstance* sub : pc.GetCompartments().GetLiquidCompartmentSubstances())
    {
      double totalMass = 0.0;
      std::string name;
      for (SELiquidCompartment* cmpt : Graph.GetCompartments())
      {
        totalMass += cmpt->GetSubstanceQuantity(*sub)->GetMass(MassUnit::g);
        name = cmpt->GetName();
      }
      if(sub->GetName()=="Hemoglobin") trk.Track("Unbound" + sub->GetName() + "_TotalMass_" + "g", time, totalMass);
      else trk.Track(sub->GetName() + "_TotalMass_" + "g", time, totalMass);
    }

    //Total Hb
    double Hb_Mols = TotalHbMols(Graph, Hb, HbO2, HbO2CO2, HbCO2);
    trk.Track("HbTotal_mols", time, Hb_Mols);

    //Total molarity
    for (SELiquidCompartment* cmpt : Graph.GetCompartments())
    {
      SELiquidSubstanceQuantity* sqO2 = cmpt->GetSubstanceQuantity(O2);
      SELiquidSubstanceQuantity* sqHbO2 = cmpt->GetSubstanceQuantity(HbO2);
      SELiquidSubstanceQuantity* sqHbO2CO2 = cmpt->GetSubstanceQuantity(HbO2CO2);

      SELiquidSubstanceQuantity* sqCO2 = cmpt->GetSubstanceQuantity(CO2);
      SELiquidSubstanceQuantity* sqHbCO2 = cmpt->GetSubstanceQuantity(HbCO2);
      SELiquidSubstanceQuantity* sqHCO3 = cmpt->GetSubstanceQuantity(HCO3);
      trk.Track(cmpt->GetName() + "_TotalOxygenMolarConcentration_mmol_per_L", time, sqO2->GetMolarity(AmountPerVolumeUnit::mmol_Per_L) + sqHbO2->GetMolarity(AmountPerVolumeUnit::mmol_Per_L) + sqHbO2CO2->GetMolarity(AmountPerVolumeUnit::mmol_Per_L));
      trk.Track(cmpt->GetName() + "_TotalCarbonDioxideMolarConcentration_mmol_per_L", time, sqCO2->GetMolarity(AmountPerVolumeUnit::mmol_Per_L) + sqHCO3->GetMolarity(AmountPerVolumeUnit::mmol_Per_L) + sqHbCO2->GetMolarity(AmountPerVolumeUnit::mmol_Per_L) + sqHbO2CO2->GetMolarity(AmountPerVolumeUnit::mmol_Per_L));
    }

    //Tissue and total oxygen w/ diffusion
    if (usingDiffusion)
    {
      cTissue.Balance(BalanceLiquidBy::Mass);
      cCapillaries.Balance(BalanceLiquidBy::Mass);
      SELiquidSubstanceQuantity* sqO2Tissue = cTissue.GetSubstanceQuantity(O2);
      trk.Track("Tissue_OxygenMass_ug", time, sqO2Tissue->GetMass(MassUnit::ug));
      trk.Track("Tissue_OxygenConcentration_ug_per_mL", time, sqO2Tissue->GetMass(MassUnit::ug)/cTissue.GetVolume(VolumeUnit::mL));

      double totalMassO2 = 0;

      for (SELiquidCompartment* cmpt : Graph.GetCompartments())
      {
        totalMassO2 += cmpt->GetSubstanceQuantity(O2)->GetMass(MassUnit::ug);
      }

      totalMassO2 += sqO2Tissue->GetMass(MassUnit::ug);
      trk.Track("xTotalOxygenMass_ug", time, totalMassO2);

    }

    if (i == 0)
    {
      trk.CreateFile(std::string(rptDirectory + outputName + ".csv").c_str(), file);
    }
    trk.StreamTrackToFile(file);

    time += deltaT_s;
  }
  file.close();
}

void PulseEngineTest::FourCompartmentTestSimple(const std::string& sOutputDirectory)
{
  FourCompartmentTest(false, false, false, sOutputDirectory);
}
void PulseEngineTest::AcidBaseFourCompartmentTest(const std::string& sOutputDirectory)
{
  FourCompartmentTest(true, false, false, sOutputDirectory);
}
void PulseEngineTest::FiveCompartmentTestWithDiffusion(const std::string& sOutputDirectory)
{
  FourCompartmentTest(false, false, true, sOutputDirectory);
}
void PulseEngineTest::AcidBaseFourCompartmentTestWithProductionConsumption(const std::string& sOutputDirectory)
{
  FourCompartmentTest(true, true, false, sOutputDirectory);
}
void PulseEngineTest::AcidBaseFiveCompartmentTestWithDiffusion(const std::string& sOutputDirectory)
{
  FourCompartmentTest(true, false, true, sOutputDirectory);
}
void PulseEngineTest::AcidBaseFiveCompartmentTestWithProductionConsumptionAndDiffusion(const std::string& sOutputDirectory)
{
  FourCompartmentTest(true, true, true, sOutputDirectory);
}