function drug_dosage ( )

%*****************************************************************************80
%
%% DRUG_DOSAGE models the variation in blood levels of a medicinal drug.
%
%  Discussion:
%
%    Drug levels in the blood, for certain drugs, can be modeled by
%    coupled first-order linear differential equations.
%
%    If we define:
%
%      A = medicine_in_intestines,
%      B = plasma_level. 
%      C = excreted medicine 
%
%    then the governing differential equations are:
%
%      dA/dt = intake - absorption_rate * A
%      dB/dt =        + absorption_rate * A - excretion_rate * B
%      dC/dt =                              + excretion_rate * B
%
%  Modified:
%
%    29 July 2019
%
%  Reference:
%
%    Soraya Dhillon, Andrzej Kostrzewski,
%    Clinical Pharmacokinetics,
%    Pharamceutical Press, 2006,
%    ISBN: 978-0-85369-571-4.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'DRUG_DOSAGE:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Simulate the variation over time of the concentration\n' );
  fprintf ( 1, '  of a medicinal drug in the bloodstream.\n' );

  dt = 0.25;
  end_time = 48.0;
  doses_per_day = 1;
  dosage_per_day = 7500;

% doses_per_day = 4;
% dosage_per_day = 12000;

  absorption_rate = 0.25;
  blood_volume = 4.6;
  medicinal_level = 800;
  toxic_level = 1000;
  N = floor ( ( end_time ) / dt );
  dosage = dosage_per_day / doses_per_day;
  dosage_interval = 24 / doses_per_day;
  excretion_rate = log ( 2.0 ) / 6.0;
  steps_between_doses = dosage_interval / dt;

  t = zeros ( N + 1, 1 );
  A = zeros ( N + 1, 1 );
  B = zeros ( N + 1, 1 );
  C = zeros ( N + 1, 1 );
%
%  Initialize variables;
%
  step_next_dose = 1;

  for j = 1 : N

    t(j+1) = t(j) + dt;

    if ( j == step_next_dose )
      intake = dosage;
      step_next_dose = j + steps_between_doses;
    else
      intake = 0;
    end

    AtoB = dt * absorption_rate * A(j);
    BtoC = dt * excretion_rate * B(j);

    A(j+1) = A(j) + intake - AtoB;
    B(j+1) = B(j)          + AtoB - BtoC;
    C(j+1) = C(j)                 + BtoC;

  end
%   
%  Plot the results
%
  clf ( );
  hold ( 'on' );

    plot ( t, B / blood_volume, 'b-', 'linewidth', 2 );

    times = [ t(1), t(N+1) ];
    m_level = [ medicinal_level, medicinal_level ];
    plot ( times, m_level, 'g-', 'linewidth', 4 );

    t_level = [ toxic_level, toxic_level ];
    plot ( times, t_level, 'r-', 'linewidth', 4 );

  hold ( 'off' );

  legend ( 'Bloodstream Level', 'Medicinal Level', 'Toxic Level', ...
    'Location', 'southeast' );

  grid ( 'on' )
  title ( 'Drug concentration', 'Fontsize', 24 );
  xlabel ( '<-- Hours -->', 'Fontsize', 16 );
  ylabel ( '<-- Concentration (mg/l) -->', 'Fontsize', 16 );
  ylim ( [ 0.0, 1100.0] );
  print ( '-dpng', 'dosage7500_p1.png' );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'DRUG_DOSAGE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
