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 a pair of
%    coupled first-order linear differential equations.
%
%    If we define:
%
%      A = medicine_in_intestines,
%      B = plasma_level.  
%
%    then the governing differential equations are:
%
%      dA/dt = - absorption_rate * A + intake
%      dB/dt = - excretion_rate * B + absorption_rate * A
%
%  Reference:
%
%    Soraya Dhillon, Andrzej Kostrzewski,
%    Clinical Pharmacokinetics,
%    Pharamceutical Press, 2006,
%    ISBN: 978-0-85369-571-4.
%
%  Local parameters:
%
%    time_step: simulation time step in hours
%    start_time: starting time in hours
%    end_time: ending time in hours
%    doses_per_day: number of daily doses
%    dosage_per_day: mg of the drug
%    absorption_rate: absorption rate constant
%    half_life: half life in hours
%    blood_volume: blood volume in liters
%    medicinal_level: level at which drug becomes effective mg/l
%    toxic_level: level at which drug becomes toxic  mg/l
%    N: number of simulation steps
%    dosage: amount of each dose
%    dosage_interval: interval between doses in hours
%    excretion_rate: excretion rate from plasma
%    steps_between_doses: number of time steps between doses
%    t: time in hours
%    intake: dose at this time step
%    medicine_in_intestines: level of drug in intestines
%    plasma_level: level of drug in the blood
%    plasma_concentration: concentration of drug in blood
%    absorption: amount absorbed at this time step
%    excretion: amount excreted at this time step
%
  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' );

  time_step = 0.25;
  start_time = 0;
  end_time = 48;
  doses_per_day = 1;
  dosage_per_day = 7500;
  absorption_rate = 0.25;
  half_life = 6;
  blood_volume = 4.6;
  medicinal_level = 800;
  toxic_level = 1000;
  N = floor ( ( end_time - start_time ) / time_step );
  dosage = dosage_per_day / doses_per_day;
  dosage_interval = 24 / doses_per_day;
  excretion_rate = log ( 2.0 ) / half_life;
  steps_between_doses = dosage_interval / time_step;

  t = zeros ( N + 1, 1 );
  intake = zeros ( N + 1, 1 );
  medicine_in_intestines = zeros ( N + 1, 1 );
  plasma_level = zeros ( N + 1, 1 );
  plasma_concentration = zeros ( N + 1, 1 );
  absorption = zeros ( N + 1, 1 );
  excretion = zeros ( N + 1, 1 );
%
%  Initialize variables;
%
  t(1) = start_time;
%
%  Give the first dose at the start time.
%
  step_next_dose = 1;
%
%  Loop to calculate intake of drug over time, using Euler's method.
%
%  This part calculates when each dose of the drug is added, based on
%  the number of doses over a 24 hour period. The number of doses must be an
%  integer divisor of 24 for this to work.
%
  for j = 1 : N

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

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

    absorption(j) = time_step * absorption_rate * medicine_in_intestines(j);
    excretion(j) = time_step * excretion_rate * plasma_level(j);

    plasma_level(j+1) = plasma_level(j) - excretion(j) + absorption(j);
    plasma_concentration(j+1) = plasma_level(j+1) / blood_volume;
    medicine_in_intestines(j+1) = medicine_in_intestines(j) - absorption(j) + intake(j);

  end
%
%  Write plasma level to a file.
%
  filename = 'drug_dosage.txt';
  output = fopen ( filename, 'wt' );
  for j = 1 : N + 1
    fprintf ( output, '%g  %g\n', t(j), plasma_concentration(j) );
  end
  fclose ( output );
%   
%  Plot the results
%
  clf
  hold on

    plot ( t, plasma_concentration, '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 ( 'Plasma Concentration',  'Medicinal Level', 'Toxic Level' );

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

  return
end
