{These functions perform blood gas calculations.  They are based on the
 papers of Gabel '81, Kelman '67, Lutz '75 and Loeppky ' 83.

 By C. S. Tritt
 01/09/86 vs. 3.02}

procedure phpco2si(var slope, intercep: real; be: real; var error: boolean);

{Returns the slope and intercept of the ph vs. ln(Pco2 / 42.5 mmHg) plot}

begin
  slope := -0.2884;
  if (be < -10.0) or (be > 20.0) then error := true else error := false;
  intercep := 7.3674 + 0.01320 * be;
end; {phpco2si}

function fpkcalc(temp, ph: real; var error: boolean): real;

{Returns the pk of H2CO3 corrected for temperature (C) and pH.  Uses Kelman's
 method.}

const

  normph = 7.4;
  normt = 38.0;
  k1 = 6.086;
  k2 = 0.042;
  k3 = 0.00472;
  k4 = 0.00139;

var

  deltaph, deltat: real;

  begin
    deltaph := normph - ph;
    deltat := normt - temp;
    if (abs(deltaph) > 2.5) or (abs(deltat) > 3.0) then
    begin
      error := true;
      fpkcalc := 6.1;
    end
    else
    begin
      error := false;
      fpkcalc := k1 + k2 * deltaph + deltat * (k3 + k4 * deltaph);
    end;
  end; {fpkcalc}

function fsolco2(temp: real; var error: boolean): real;

{Returns the solubility (mM/l/mmHg) of CO2 in plasma corrected for temperature
   (C).  Uses the method of Kelman}

const

  normt = 37.0;
  k1 = 0.0307;
  k2 = 0.00057;
  k3 = 0.00002;

var

  deltat: real;

begin
  deltat := normt - temp;
  if abs(deltat) > 3.0 then
  begin
    error := true;
    fsolco2 := 0.03;
  end
  else
  begin
    error := false;
    fsolco2 := k1 + k2 * deltat + k3 * sqr(deltat);
  end;
end;

function fbicarb(pco2, ph, pkco2, solco2: real; var error: boolean): real;

{Returns bicarbinate concentration in mVal/l for PCO2 in mmHg and solco2 in
   mM/l/mmHg.  Uses the method of Lutz.}

const

  zer = 0.0;
  ln10 = 2.30259;

var

  phlesspk: real;

begin
  phlesspk := ph - pkco2;
  if (pco2 < 0.0) or (abs(phlesspk) > 5.0) then
  begin
    error := true;
    fbicarb := 0.0;
  end
  else
  begin
    error := false;
    fbicarb := solco2 * pco2 * exp(ln10 * phlesspk);
  end;
end; { fbicarb}

function fbaseex(bicarb, ph: real; var error: boolean): real;

{Returns calculated base excess and bicarb in mVal/l.  Uses the method of
   Lutz.}

const

  hco3fact = 0.775;
  phfact = 26.88;
  constant = 217.9;
  zer = 0.0;
  nin = 99.9;

begin
  if (bicarb < 0.0) or (bicarb > 75.0) or (abs(ph - 7.4) > 2.5) then
  begin
    error := true;
    fbaseex := 0.0;
  end
  else
  begin
    error := false;
    fbaseex := hco3fact * bicarb + phfact * ph - constant;
  end;
end; {fbaseex}

function fpo2cor(po2, ph, be, temp, bohr: real; var error: boolean): real;

{Returns PO2, in mmHg, corrected for Bohr effect, pH, BE (mVal/l) and
 temperature (C).  Uses the method of Lutz.}

const

  normph = 7.4;
  befac = 0.0013;
  tempfac = 0.031;
  normtemp = 37.0;
  zer = 0.0;
  nin = 99.9;
  ln10 = 2.30259;

var

  deltat, deltaph, term: real;

begin
  deltat := temp - normtemp;
  deltaph := ph - normph;
  term :=  bohr * deltaph + befac * be - tempfac * deltat;
  if po2 < 0.0 then po2 := 0.0;
  if (abs(deltaph) > 2.5) or (abs(be) > 25.0) or (abs(deltat) > 3.0) or
    (term > 10.0) or (po2 = 0.0) then
  begin
    error := true;
    fpo2cor := po2;
  end
  else
  begin
    error := false;
    fpo2cor := po2 * exp(ln10 * term);
  end;
end; {fpo2cor}

type

parmlist = array[1..10] of real;

function fo2sat(po2: real; var parm: parmlist; var error: boolean): real;

{Returns percentage hemeglobin saturation. PO2 in mmHg and corrected for Bohr
  effect, pH, BE (mVal/l) and temperature (C).  Uses the method of Lutz.  The
  parameters for the species of interest are passed in the array parmlist.}

const

  cent = 100.0;
  one = 1.0;

var

  o2sat: real;

begin
  if po2 <= 0.0 then
  begin
    error := true;
    fo2sat := 0.0;
  end
  else
  begin
    error := false;
    o2sat := parm[4] * sqr(po2) + parm[5] * po2;
    if o2sat <= parm[2] then fo2sat := o2sat
    else o2sat := parm[6] * po2 - parm[7];
    if o2sat < parm[3] then fo2sat := o2sat
    else fo2sat :=  cent / (exp(ln(parm[8] / po2) * parm[9]) + one);
  end;
end; {fo2sat}

function fo2con(po2, so2, hb: real; var error: boolean): real;

{Returns O2 concentration in blood - ml 02 /100 ml  blood.  PO2 in mmHg,
  SO2 in percent and hb in gm/dl.  Uses Gabel's method.}

const

  solo2 = 0.003;
  caphb = 0.0136;
  zer = 0.0;

begin
  if po2 < 0.0 then
  begin
    error := true;
    fo2con := 0.0;
  end
  else
  begin
    error := false;
    fo2con := solo2 * po2 + caphb * so2 * hb;
  end;
end; {fo2con}

function fco2conk(pco2, ph, o2sat, hct, pk, sol: real; var error: boolean):
  real;

{Returns CO2 concentration in blood - ml CO2 / 100 ml blood. PCO2 in mmHg,
  O2sat in percent, hct in percent and sol in mM/l.  Uses Kelmen's method.}

const

  normph = 7.4;
  doxa = 0.590;
  doxb = 0.2913;
  doxc = 0.0844;
  dra = 0.664;
  drb = 0.2275;
  drc = 0.0938;
  ln10 = 2.30259;
  unitcorr = 2.22;

var

 devph, phlesspk, dox, dr, o2satfrc, d, cp, cc, hctfract: real;

begin
  devph := normph - ph;
  phlesspk := ph - pk;
  if (pco2 < 0.0) or (abs(devph) > 2.5) or (abs(hct - 30.0) > 20.0)  or
    (abs(phlesspk) > 10.0) then
  begin
    error := true;
    if pco2 < 0.0 then pco2 := 0.0;
    if ph < 0.0 then ph := 0.0
      else if ph > 14.0 then ph := 14.0;
    fco2conk := unitcorr * 0.03 * pco2 * ( 1.0 + exp(ln10 * (ph - 6.1)));
  end
  else
  begin
    error := false;
    dox := doxa + doxb * devph - doxc * sqr(devph);
    dr := dra + drb * devph - drc * sqr(devph);
    o2satfrc := o2sat / 100.0;
    d := dox * o2satfrc + dr * (1.0 - o2satfrc);
    cp := sol * pco2 * (1.0 + exp(ln10 * phlesspk));
    cc := d * cp;
    hctfract := hct / 100.0;
    fco2conk := (cc * hctfract + cp * (1.0 - hctfract)) * unitcorr;
  end;
end; {fco2conk}

function fco2conl(pco2, sol, bicarb: real; var error: boolean): real;

{Returns the total CO2 content of blood based on a modification Lutz's
 method.  PCO2 is in mmHg, sol is in mmol/l/mmHg and bicarb is in mEq/l.
 The result is in volume %.}

const

  unitcorr = 2.22;

begin
  if (bicarb < 0.0) or (bicarb > 75.0) or (pco2 < 0.0) then
  begin
    error := true;
    if pco2 < 0.0 then pco2 := 0.0;
    if bicarb < 0.0 then bicarb := 0.0;
    fco2conl := unitcorr * (0.03 * pco2 + bicarb);
  end
  else
  begin
    error := false;
    fco2conl := (bicarb + sol * pco2) * unitcorr;
  end;
end; {fco2conl}

function fco2conf(pco2, ph, o2sat, hb, pk, sol: real; var error: boolean):
  real;

{Returns the total CO2 content of blood based on a modification of
 Loeppky's method.  PCO2 is in mmHg, O2sat is in %, Hb is in g%, and
 Sol is in mmol/l/mmHg.  The result is in volume %.}

const

  unitcorr = 2.226;
  ln10 = 2.30259;
  hbfact = 0.02924;
  satfact1 = 2.244;
  satfact2 = 0.00422;
  phfact = 8.740;

var

  plasco2, phlesspk: real;

begin
  error := false;
  phlesspk := ph - pk;
  if (pco2 < 0.0) or (abs(ph - 7.4) > 5.0) or (abs(hb - 15.0) > 15.0) or
    (abs(phlesspk) > 10.0) then
  begin
    error := true;
    if pco2 < 0.0 then pco2 := 0.0;
    if ph < 0.0 then ph := 0.0
      else if ph > 14.0 then ph := 14.0;
    fco2conf := unitcorr * 0.03 * pco2 * ( 1.0 + exp(ln10 * (ph - 6.1)));
  end
  else
  begin
    error := false;
    plasco2 := sol * pco2 * (1.0 + exp(ln10 * phlesspk));
    fco2conf := (1.0 - hbfact * hb / (satfact1 - satfact2 * o2sat) /
      (phfact - ph)) * plasco2 * unitcorr;
  end;
end; {fco2conf}

function fgasxtop(percent, pressure: real): real;

  {This function converts gas mole percent to partial pressure}

  begin
    fgasxtop := percent / 100.0 * pressure;
  end; {fgasxtop}

procedure bldcon(var o2con, co2con, baseex: real; po2, pco2, ph, temp, hct,
 hbgm100, hbcor: real; co2meth: integer; var parm: parmlist; var mastflag:
 boolean);

  {This function calculates the gas contents of the blood.  For units see the
     functions and procedures called.}

var

  pkco2, solco2, bicarb, corrpo2, o2sat: real;

  errflag: boolean;

begin
  mastflag := false;
  pkco2 := fpkcalc(temp, ph, errflag);
  mastflag := mastflag or errflag;
  solco2 := fsolco2(temp, errflag);
  mastflag := mastflag or errflag;
  bicarb := fbicarb(pco2, ph, pkco2, solco2, errflag);
  mastflag := mastflag or errflag;
  baseex := fbaseex(bicarb, ph, errflag);
  mastflag := mastflag or errflag;
  corrpo2 := fpo2cor(po2, ph, baseex, temp, parm[1], errflag);
  mastflag := mastflag or errflag;
  o2sat := fo2sat(corrpo2, parm, errflag);
  mastflag := mastflag or errflag;
  o2con := fo2con(po2, o2sat, hbcor, errflag);
  mastflag := mastflag or errflag;
  case co2meth of
    1 : begin
          co2con := fco2conl(pco2, solco2, bicarb, errflag);
          mastflag := mastflag or errflag;
        end;
    2 : begin
          co2con := fco2conk(pco2, ph, o2sat, hct, pkco2, solco2, errflag);
          mastflag := mastflag or errflag;
        end;
   3 : begin
          co2con := fco2conf(pco2, ph, o2sat, hbgm100, pkco2, solco2,
            errflag);
          mastflag := mastflag or errflag;
        end;
  end;
end; {bldcon}

function fphcalc(pco2, slope, intercep, be: real): real;

  {This function returns the estimated blood pH at a given PCO2.  A linear
   relationship between pH and ln(PCO2/42.5) is used.  At very low pCO2, a
   limiting pH calculated from the equation of Thews given in Lutz's paper
   is returned.  Found to estimate the pH of pig blood for PCO2 between 10 and
   45 mmHg too high.}

var

  phmax, mxlessin: real;

begin
  phmax := (be + 217.9) / 26.88;
  mxlessin := phmax - intercep;
  if abs(mxlessin) > 10.0 then mxlessin := 10 * mxlessin / abs(mxlessin);
  if pco2 < 42.5 * exp(mxlessin / slope) then fphcalc := phmax
  else fphcalc := slope * ln(pco2 / 42.5) + intercep;

end; {fphcalc}

function fphcalc2(pco2, slope, intercep, be, temp: real; becnflag: boolean;
  var error: boolean): real;

  {This function returns the estimated blood pH at a given PCO2 (mmHg).
     Newton's method is used, assuming a constant base excess and using Thews'
     equation as given in Lutz '75.}

  var

    phold, phnew, solco2, pk, const1, const2, ln10: real;

    done: boolean;

    loops, maxloops: integer;

  begin
    solco2 := fsolco2(temp, error);
    phnew := fphcalc(pco2, slope, intercep, be);
    if becnflag then
    begin
      ln10 := ln(10.0);
      const1 := 0.775 * solco2 * pco2;
      maxloops := 15;
      loops := 0;
      repeat
        loops := loops + 1;
        phold := phnew;
        if phold > 14.0 then
        begin
          phold := 14.0;
          loops := maxloops;
        end;
        if phold < 0.0 then
        begin
          phold := 0.0;
          loops := maxloops;
        end;
        pk := fpkcalc(temp, phold, error);
        if error then loops := maxloops;
        const2 := const1 / exp(ln10 * pk);
        phnew := phold - (const2 * exp(ln10 * phold) + 26.88 * phold - 217.9 -
          be) / (ln10 * const2 * exp(ln10 * phold) + 26.88);
        done := loops = maxloops;
      until done or (abs(phold - phnew) < 0.0005);
      if done then error := true;
    end;
    fphcalc2 := phnew;
  end; {fphcalc2}

function fshunt(cin, cout, couteq: real): real;

{This function calculate and returns the content based shunt for the general
   case.  All contents must be in the same (arbitary) units.}
  begin
    fshunt := (cout - couteq) / (cin - couteq);
  end; {fshunt}

  function fph2o(t: real): real;

  {This function returns the vapor pressure of water (in mmHg) at the given
   temperature (in C.  For use between 30.0 and 50.0 C.}

  begin
    if t > 50.0  then t := 50.0;
    if t < 30.0  then t := 30.0;
    fph2o := 5.9562 * exp(0.05558 * t);
  end; {fph2o}

function fconout(conin, flow, fluxout: real): real;

  {This function returns the outlet concentration of a component given the
   inlet concentration, the total flow rate and the flux of the component out
   of the flow using a mass balance.  All concentrations in volume percent.}

  var

    term: real;

  begin
    term := conin / 100.0 * flow - fluxout;
    if term < 0.0 then fconout := 0.0
      else fconout := term * 100.0 / flow;
  end;

procedure delo2(var def, esto2sat, toler: real; po2, o2con, ph, be, hb, temp:
    real; var parm: parmlist; var masterr: boolean);

  {This procedure returns the difference between the calculated O2 content at
   the given PO2 (mmHg) and the given O2 content.  If the calculated hb - O2
   saturation is near a discountinuity in the Lutz curves the convergance
   tolerance is increased and the error flag is set.  In all cases the O2
   saturation is also returned.  All contents in volume percent.}

  var

    cpo2: real;

    error: boolean;

  begin
    masterr := true;
    if o2con < 0.0 then o2con := 0.0
    else masterr := false;
    cpo2 := fpo2cor(po2, ph, be, temp, parm[1], error);
    masterr := masterr or error;
    esto2sat := fo2sat(cpo2, parm, error);
    masterr := masterr or error;
    if (abs(esto2sat - parm[2]) < parm[10]) or (abs(esto2sat - parm[3]) <
      parm[10]) then
      begin
        toler := 2.0 * parm[10] / parm[6];
        masterr := true;
      end;
    def := fo2con(po2, esto2sat, hb, error) - o2con;
    masterr := masterr or error;
  end; {fdelo2}

procedure cntopo2(var po2, o2sat, delnew: real; delpo2, o2con, ph, be, hb,
  temp: real; var parm: parmlist; var masterr: boolean);

  {This procedure calculates the PO2 (mmHg) and hb - O2 saturation (%) for a
   given O2 content (Vol %).}

  var

    delold, toler: real;

    loops, maxloops: integer;

    error: boolean;

  begin
    if o2con <= 0.0 then
    begin
      po2 := 0.0;
      o2sat := po2;
      delnew := -o2con;
      masterr := true;
    end
    else
    begin
      toler := 0.1;
      masterr := false;
      delo2(delnew, o2sat, toler, po2, o2con, ph, be, hb, temp, parm, error);
      masterr := masterr or error;
      maxloops := 20;
      loops := 0;
      repeat
        loops := loops + 1;
        delold := delnew;
        if delnew < 0.0 then po2 := po2 + delpo2
          else po2 := po2 - delpo2;
        delo2(delnew, o2sat, toler, po2, o2con, ph, be, hb, temp, parm, error);
        masterr := masterr or error;
        if (delnew * delold) < 0.0 then delpo2 := delpo2 / 2.0;
      until (delpo2 < toler) or (loops = maxloops);
      if loops = maxloops then masterr := true
        else if delnew < 0.0  then po2 := po2 + delpo2
          else po2 := po2 - delpo2;
    end;
  end; {cntopo2}

  procedure cntopco2(var pco2, phnew, delnew: real; delpco2, co2con, o2sat,
   hct, be, temp: real; phmeth2: boolean; var masterr: boolean);

  {This procedure calculates the PCO2 (mmHg) and ph for a given CO2 content
    (Vol %) and base excess (mEq/l).}

  var

    delold, phold, sol, pk, slope, intercep: real;

    loops: integer;

    error: boolean;

  begin
    if co2con <= 0.0 then
    begin
      pco2 := 0.0;
      masterr := true;
    end
    else
    begin
      masterr := false;
      phpco2si(slope, intercep, be, error);
      masterr := masterr or error;
      pk := fpkcalc(temp, phnew, error);
      masterr := masterr or error;
      sol:= fsolco2(temp, error);
      masterr := masterr or error;
      delnew := fco2conk(pco2, phnew, o2sat, hct, pk, sol, error) - co2con;
      masterr := masterr or error;
      loops := 0;
      repeat
        loops := loops + 1;
        delold := delnew;
        phold := phnew;
        if delnew < 0.0 then pco2 := pco2 + delpco2
          else pco2 := pco2 - delpco2;
        phnew := fphcalc2(pco2, slope, intercep, be, temp, phmeth2, error);
        masterr := masterr or error;
        pk := fpkcalc(temp, phnew, error);
        masterr := masterr or error;
        delnew := fco2conk(pco2, phnew, o2sat, hct, pk, sol, error) - co2con;
        masterr := masterr or error;
        if (delnew * delold) < 0.0 then delpco2 := delpco2 / 2.0;

      until ((abs(delnew) < 0.1) and (delpco2 < 0.1) and (abs(phnew-phold) <
        0.005)) or (loops = 20);
      if (loops = 20) then masterr := true
        else if delnew < 0.0 then pco2 := pco2 + delpco2
          else pco2 := pco2 - delpco2;
    end;
  end; {cntopco2}

  procedure contops(var po2new, pco2new, phnew, erro2, errco2: real; delpo2,
    delpco2, be, temp, hb, hct, o2con, co2con: real; var parm: parmlist;
    phmeth2: boolean; var masterr: boolean);

  {This procedure calculates PO2, PCO2 and pH for given base excess, oxygen
   content and carbon dioxide content conditions.}

  var

    po2old, pco2old, phold, o2sat: real;

    loopscn, maxloops: integer;

    error: boolean;

  begin
    maxloops := 20;
    loopscn := 0;
    repeat
      masterr := false;
      loopscn := loopscn + 1;
      po2old := po2new;
      pco2old := pco2new;
      phold := phnew;
      cntopo2(po2new, o2sat, erro2, delpo2, o2con, phnew, be, hb, temp, parm,
        error);
      masterr := masterr or error;
      cntopco2(pco2new, phnew, errco2, delpco2, co2con, o2sat, hct, be, temp,
        phmeth2, error);
      masterr := masterr or error;

      delpo2 := 10.0; {Reset step size after first calls}
      delpco2 := 5.0; {Reset step size after first calls}

    until ((abs(po2old - po2new) < 1.0) and (abs(pco2old - pco2new) < 0.1) and
      (abs(phold - phnew) < 0.005)) or (loopscn = maxloops);

    if loopscn = maxloops then masterr := true;
  end; {cntops}

procedure outgas(var gasflwout, gasopco2, gasopo2: real; gasflwin, flxco2,
  flxo2, pcorfact, h2ofact, baropres, qnco2, qno2: real; var error: boolean);

{This procedure calculate the outlet gas flow corrected for temperature,
 pressure and water vapor and the outlet gas pco2 and po2 based on the inlet
 quantities and the transfer rates.}

begin
  gasflwout := (gasflwin + flxco2 - flxo2) * pcorfact * h2ofact;
  gasopco2 := baropres * (flxco2 + qnco2) * pcorfact * h2ofact /
    gasflwout;
  gasopo2 := baropres * (qno2 - flxo2) * pcorfact * h2ofact /
    gasflwout;
  if (gasopco2 < 0.0) or (gasopo2 < 0.0) then
  begin
    gasflwout := gasflwin;
    gasopco2 := baropres * qnco2 * pcorfact * h2ofact / gasflwout;
    gasopo2 := baropres * qno2 * pcorfact * h2ofact / gasflwout;
    error := true;
  end
  else
    error := false;
end; {outgas}