{$N+}
{$R+}

program lungdp66;
{
  Calculates the deposition in the five regions of the respiratory
  tract using the ICRP 66 respiratory tract model.  Written in
  Turbo Pascal V6.0.

  T. R. La Bone
  22 January 1998
}

uses wincrt;

const
   chi = 1.5;   {particle shape factor}
   rho = 3.0;  {density of particle}

type
  depo_record = record
    eta : double; {filtration efficiency for a region}
    phi : double; {fraction of air reaching a region}
    xi : double;  {ratio of phi for two filters}
    DE : double;  {deposition in a region}
  end;
  depo_type = (ae,th);
  depo_array = array[1..5] of double;
var
  ouf : text;
  dep : depo_array;
  sum : double;
  i,j : integer;
  BB1s,                {slow-clearance fraction deposited in BB}
  bb2s,                {slow-clearance fraction deposited in bb}
  amad,                {activity median aerodynamic diameter}
  s,                   {geometric standard deviation of the amad}
  amtd,                {activity median thermodynamic diameter}
  U,                   {wind speed in m/s}
  V,                   {total air flow rate through nose and mouth}
  Vn,                  {air flow rate through the nose}
  FRC,                 {functional residual capacity of the lung}
  VT,                  {tidal volume of the lung}
  VDet,                {anatomical dead space of the ET region}
  VDbb2,               {anatomical dead space of the bb region}
  VDBB1,               {anatomical dead space of the BB region}
  SFA,                 {scaling factor for the bronchiole}
  SFt,                 {scaling factor for the trachea}
  SFb : double;        {scaling factor for the bronchiolar airways}
  poly_deposition : depo_array; {deposition in five regions of RT}

(**********************************************************************)

function sexp( x:double ) : double;
{ Prevents runtime errors caused by exp underflowing. }
begin
   if x < -2000 then
      sexp := 0.0
   else
      sexp := exp(x);
end;

function rse( a,b:double ) : double;
{ Raises a to the b power. }
var
  sign : integer;
  c : double;
begin
  if (a > 0.0) then
    rse := sexp(b*ln(a))
  else if a = 0.0 then
    rse := 0
  else
    begin
      c := trunc(b);
      if b <> c then begin
        write('error in function rse');
        halt;
      end;
      if odd(trunc(b)) then sign := -1 else sign := 1;
        rse := sign*sexp(b*ln(abs(a)));
    end;
end; {rse}

function log( x:double ) : double;
{ Calculates the common (base 10) logarithm of x. }
begin
  log := ln(x) / ln(10.0);
end; {log}

function alog( x:double ) : double;
{ Calculates 10 raised to the power of x. }
begin
  alog := sexp(x*ln(10));
end; {alog}

function c( d:double ) : double;
{ Slip correction factor, Equation D.6 in ICRP 66. }
begin
  c := 1.0 + 0.0683/d*(2.514+0.8*sexp(-0.55*d/0.0683));
end; {c}

procedure calc_diameter( var dt:double; da:double );
{ Calculates the thermodynamic diameter from the aerodynamic diameter,
  see Equation D.13 in ICRP 66.  Uses Newton's method to find root. }
var
  new_dt,diff : double;
function k_0( da,dt:double ) : double;
begin
  k_0 := sqr(dt)*c(dt) - sqr(da)*c(da)*chi/rho;
end;
function k_1( dt:double ) : double;
begin
  k_1 := 0.8*0.0683*sexp(-0.55*dt/0.0683) + 2.0*dt + 2.514*0.0683
       - 0.55*0.8*dt*sexp(-0.55*dt/0.0683);
end;
begin
  dt := da;
  repeat
    new_dt := dt - k_0(da,dt) / k_1(dt);
    diff := abs(new_dt - dt)/new_dt;
    dt := new_dt;
  until diff < 1.0e-15;
end; {calc_diameter}

procedure nose_deposition(
  var deposition : depo_array;
  da,
  dt,
  U,
  V,
  Vn,
  FRC,
  VT,
  VDet,
  VDbb2,
  VDBB1,
  SFA,
  SFt,
  SFb : double );
{
  The procedure nose_deposition calculates the deposition in the five
  regions of the respiratory track following an inhalation of a
  monodisperse particulate material through the nose.
}
var
  depo : array[1..9] of depo_record;
  eta,a,R,p : array[depo_type,1..9] of double;
  i : integer;
  V_Dbb2,
  V_DBB1,
  D,
  Psith,
  tB1,
  tb2,
  tA : double;
begin
  for i := 1 to 5 do deposition[i] := 0.0;
  V_Dbb2 := VDbb2 * (1 + VT/FRC);
  V_DBB1 := VDBB1 * (1 + VT/FRC);
  D := 1.38e-16*310*c(dt) / (3*Pi*1.88e-8*dt);
  Psith := 1 + 100*sexp(-Sqr(log(100+10/rse(dt,0.9))));
  tB1 := VDBB1/V * (1 + 0.5*VT/FRC);
  tA := ( VT - VDet - (VDBB1+VDbb2)*(1 + VT/FRC) ) / V;
  tb2 := VDbb2/V * (1 + 0.5*VT/FRC);

{ assign parameters from Table 12 of ICRP 66 }
  a[ae,1] := 3.0e-4;
  a[ae,2] := 5.5e-5;
  a[ae,3] := 4.08e-6;
  a[ae,4] := 0.1147;
  a[ae,5] := 0.146*rse(SFA,0.98);
  a[ae,6] := 0.1147;
  a[ae,7] := 2.04e-6;
  a[ae,8] := 5.5e-5;
  a[ae,9] := 3.0e-4;
  a[th,1] := 18.0;
  a[th,2] := 15.1;
  a[th,3] := 22.02*rse(SFt,1.24)*Psith;
  a[th,4] := -76.8 + 167.0*rse(SFb,0.65);
  a[th,5] := 170 + 103*rse(SFA,2.13);
  a[th,6] := -76.8 + 167.0*rse(SFb,0.65);
  a[th,7] := 22.02*rse(SFt,1.24)*Psith;
  a[th,8] := 15.1;
  a[th,9] := 18.0;

  R[ae,1] := Sqr(da)*Vn*rse(SFt,3);
  R[ae,2] := Sqr(da)*Vn*rse(SFt,3);
  R[ae,3] := Sqr(da)*V*rse(SFt,3);
  R[ae,4] := (0.056+rse(tb2,1.5))*rse(da,rse(tb2,-0.25));
  R[ae,5] := Sqr(da)*tA;
  R[ae,6] := (0.056+rse(tb2,1.5))*rse(da,rse(tb2,-0.25));
  R[ae,7] := Sqr(da)*V*rse(SFt,3);
  R[ae,8] := Sqr(da)*Vn*rse(SFt,3);
  R[ae,9] := Sqr(da)*Vn*rse(SFt,3);
  R[th,1] := D*rse(Vn*SFt,-0.25);
  R[th,2] := D*rse(Vn*SFt,-0.25);
  R[th,3] := D*tB1;
  R[th,4] := D*tb2;
  R[th,5] := D*tA;
  R[th,6] := D*tb2;
  R[th,7] := D*tB1;
  R[th,8] := D*rse(Vn*SFt,-0.25);
  R[th,9] := D*rse(Vn*SFt,-0.25);

  p[ae,1] := 1.0;
  p[ae,2] := 1.17;
  p[ae,3] := 1.152;
  p[ae,4] := 1.173;
  p[ae,5] := 0.6495;
  p[ae,6] := 1.173;
  p[ae,7] := 1.152;
  p[ae,8] := 1.17;
  p[ae,9] := 1.0;
  p[th,1] := 0.5;
  p[th,2] := 0.538;
  p[th,3] := 0.6391;
  p[th,4] := 0.5676;
  p[th,5] := 0.6101;
  p[th,6] := 0.5676;
  p[th,7] := 0.6391;
  p[th,8] := 0.538;
  p[th,9] := 0.5;

{ calculate phi, xi, and DE}
  depo[1].phi := 1.0;
  depo[2].phi := 1.0;
  depo[3].phi := 1 - VDet/VT;
  depo[4].phi := 1 - (VDet+V_DBB1)/VT;
  depo[5].phi := 1 - (VDet+V_DBB1+V_Dbb2)/VT;
  depo[6].phi := depo[4].phi;
  depo[7].phi := depo[3].phi;
  depo[8].phi := depo[2].phi;
  depo[9].phi := depo[1].phi;

  eta[ae,1] := 0.5*(1-1/(a[ae,1]*rse(R[ae,1],p[ae,1])+1));
  eta[ae,2] := 1-1/(a[ae,2]*rse(R[ae,2],p[ae,2])+1);
  eta[ae,3] := 1 - sexp(-a[ae,3]*rse(R[ae,3],p[ae,3]));
  eta[ae,4] := 1 - sexp(-a[ae,4]*rse(R[ae,4],p[ae,4]));
  eta[ae,5] := 1 - sexp(-a[ae,5]*rse(R[ae,5],p[ae,5]));
  eta[ae,6] := eta[ae,4];
  eta[ae,7] := 1 - sexp(-a[ae,7]*rse(R[ae,7],p[ae,7]));
  eta[ae,8] := eta[ae,2];
  eta[ae,9] := eta[ae,1];
  eta[th,1] := 0.5*(1-sexp(-a[th,1]*rse(R[th,1],p[th,1])));
  eta[th,2] := 1 - sexp(-a[th,2]*rse(R[th,2],p[th,2]));
  eta[th,3] := 1 - sexp(-a[th,3]*rse(R[th,3],p[th,3]));
  eta[th,4] := 1 - sexp(-a[th,4]*rse(R[th,4],p[th,4]));
  eta[th,5] := 1 - sexp(-a[th,5]*rse(R[th,5],p[th,5]));
  eta[th,6] := eta[th,4];
  eta[th,7] := eta[th,3];
  eta[th,8] := eta[th,2];
  eta[th,9] := eta[th,1];

  depo[1].xi := 1.0;
  for i := 2 to 9 do
    depo[i].xi := depo[i].phi / depo[i-1].phi;

  for i := 1 to 9 do { Equation 10}
    depo[i].eta := Sqrt( Sqr(eta[ae,i]) + Sqr(eta[th,i]));

  depo[1].DE := depo[1].eta { Equation 8}
    * (1 - 0.5*(1-1/(7.6e-4*rse(da,2.8)+1)) + 1e-5*rse(U,2.75)*sexp(0.055*da));

  for i := 2 to 9 do
    depo[i].DE := depo[i-1].DE * depo[i].eta * depo[i].xi * (1/depo[i-1].eta - 1);

  depo[1].DE := depo[1].DE + depo[9].DE;
  depo[2].DE := depo[2].DE + depo[8].DE;
  depo[3].DE := depo[3].DE + depo[7].DE;
  depo[4].DE := depo[4].DE + depo[6].DE;
  for i := 1 to 5 do
    deposition[i] := depo[i].DE;
end; {nose_deposition}

function integral( x,m,s:double ) : double;
{ Calculates integral of the normal PDF between -infinity and x,
  given x, the mean m, and the standard eviation s. }
var
  p,z,r : double;
begin
  z := (x-m)/s;
  r := 1/(1+0.2316419*abs(z));
  p := sexp(-z*z/2) * r*(0.319381530 + r*(-0.356563782 + r*(1.781477937
     + r*(-1.821255978 + r*1.330274429)))) / sqrt(2*Pi);
  if z > 0.0 then p := 1 - p;
  integral := p;
end; {integral}

procedure nose_polydeposition(
  var poly_deposition : depo_array;
  var BB1s : double;
  var bb2s : double;
  amad,
  amtd,
  s,
  U,
  V,
  Vn,
  FRC,
  VT,
  VDet,
  VDbb2,
  VDBB1,
  SFA,
  SFt,
  SFb : double );
{
  The procedure nose_polydeposition calculates the deposition in
  the five regions of the respiratory track following an inhalation
  of a polydisperse particulate material through the nose.
}
const
  number_increments = 500;
  z = 4.0;
var
  lower_limit,upper_limit,increment : double;
  da,dt,fs : double;
  prob : double;
  i,j : integer;
  deposition : depo_array;
begin
  BB1s := 0.0;
  bb2s := 0.0;
  lower_limit := ln(amad) - z*ln(s);
  upper_limit := ln(amad) + z*ln(s);
  increment := (upper_limit - lower_limit) / number_increments;
  for i := 1 to 5 do poly_deposition[i] := 0.0;
  for i := 1 to number_increments do begin
    prob := integral(lower_limit+increment,ln(amad),ln(s))
          - integral(lower_limit,ln(amad),ln(s));
    da := exp(lower_limit + increment/2);
    calc_diameter(dt,da);
    nose_deposition(deposition,da,dt,U,V,Vn,FRC,VT,VDet,VDbb2,VDBB1,SFA,SFt,SFb);
    for j := 1 to 5 do
      poly_deposition[j] := poly_deposition[j] + prob*deposition[j];
    if da <= 2.5*sqrt(rho/chi) then
       fs := 0.5
    else
       fs := 0.5*sexp( -ln(2)/1.1*(dt-2.5));
    BB1s := BB1s + fs*prob*deposition[3];
    bb2s := bb2s + fs*prob*deposition[4];
    lower_limit := lower_limit + increment;
  end;
end; {nose_polydeposition }

procedure sigma_g( var s:double; amad,amtd:double );
{ Calculates the geometric standard deviation of a polydisperse aerosol,
  see Equation 16 of ICRP 66. }
begin
   if amad > 1.0 then
      s := 2.5
   else if amad < 0.0005 then
      s := 1.0
   else
      s := 1 + 1.5*(1- 1/(100*rse(amtd,1.5)+1));
end; {sigma_g}

(**********************************************************************)

begin
  clrscr;

{for light exercise}
  V := 833.0;
  Vn := V;
  VT := 1250.0;

  amad := 5.0;
  U := 0.0;
  FRC := 3301.0;
  VDet := 50.0;
  VDbb2 := 47.0;
  VDBB1 := 49.0;
  SFA := 1.0;
  SFt := 1.0;
  SFb := 1.0;

  calc_diameter(amtd,amad);
  sigma_g(s,amad,amtd);

  for i := 1 to 5 do poly_deposition[i] := 0.0;
  nose_polydeposition(dep,BB1s,bb2s,amad,amtd,s,U,V,Vn,FRC,VT,VDet,VDbb2
  ,VDBB1,SFA,SFt,SFb);
  for i := 1 to 5 do poly_deposition[i] := dep[i];

  writeln('ET1 : ',(100*poly_deposition[1]):12:2);
  writeln('ET2 : ',(100*poly_deposition[2]):12:2);
  writeln('BB  : ',(100*poly_deposition[3]):12:2);
  writeln('bb  : ',(100*poly_deposition[4]):12:2);
  writeln('AI  : ',(100*poly_deposition[5]):12:2);

  sum := 0.0;
  for i := 1 to 5 do
    sum := sum + poly_deposition[i];
  writeln('Sum : ',100*sum:12:2);
  writeln;
  writeln('Fs BB : ',(BB1s/poly_deposition[3]):12:5);
  writeln('Fs bb : ',(bb2s/poly_deposition[4]):12:5);

  assign(ouf,'lungdp66.out');
  rewrite(ouf);
  writeln(ouf,'ET1 : ',(100*poly_deposition[1]):12:2);
  writeln(ouf,'ET2 : ',(100*poly_deposition[2]):12:2);
  writeln(ouf,'BB  : ',(100*poly_deposition[3]):12:2);
  writeln(ouf,'bb  : ',(100*poly_deposition[4]):12:2);
  writeln(ouf,'AI  : ',(100*poly_deposition[5]):12:2);
  writeln(ouf,'Sum : ',100*sum:12:2);
  writeln(ouf);
  writeln(ouf,'Fs BB : ',(BB1s/poly_deposition[3]):12:5);
  writeln(ouf,'Fs bb : ',(bb2s/poly_deposition[4]):12:5);
  close(ouf);


end.