function dosage_euler ( )

%*****************************************************************************80
%
%% dosage_euler() uses euler() to model medicinal blood levels.
%
%  Discussion:
%
%    Drug levels in the blood, for certain drugs, can be modeled by a set of
%    coupled first-order linear differential equations.
%
%    If we define:
%
%      A = medicine level in intestines,
%      B = medicinal level in bloodstream.  
%
%    then the governing differential equations are:
%
%      dA/dt = intake - absorption_rate * A
%      dB/dt =        + absorption_rate * A - excretion_rate * B
%
%    or, in discretized, forward Euler form:
%
%      Anew = Aold + dt * ( intake - absorption_rate * Aold )
%      Bnew = Bold + dt * (        + absorption_rate * Aold - excretion_rate * Bold )
%
%    After computation, the medicinal bloostream level is divided by the blood volume
%    to yield the medicinal bloodstream concentration.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    23 April 2021
%
%  Reference:
%
%    Soraya Dhillon, Andrzej Kostrzewski,
%    Clinical Pharmacokinetics,
%    Pharamceutical Press, 2006,
%    ISBN: 978-0-85369-571-4.
%
%  Local:
%
%    real dt: simulation time step 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
%    intake: dose at this time step
%
%    t: time in hours
%    medicine_in_intestines: level of drug in intestines
%    plasma: level of drug in blood, then converted to concentration.
%
%    real t0: the initial time, in hours.
%
%    real tstop: the final time, in hours.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'dosage_euler():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Use the Euler method to simulate the variation over time\n' );
  fprintf ( 1, '  of the level of a medicinal drug in the intestines\n' );
  fprintf ( 1, '  and the bloodstream.\n' );

  [ t0, y0, tstop ] = dosage_parameters ( );

  dt = 0.25;
  n = floor ( ( tstop - t0 ) / dt );

  doses_per_day = 1;
  dosage_per_day = 7500;
  absorption_rate = 0.25;
  half_life = 6.0;
  dosage = dosage_per_day / doses_per_day;
  dosage_interval = 24.0 / doses_per_day;
  excretion_rate = log ( 2.0 ) / half_life;
  steps_between_doses = dosage_interval / dt;
%
%  Allocate space.
%
  t = zeros ( n + 1, 1 );
  intake = zeros ( n + 1, 1 );
  medicine_in_intestines = zeros ( n + 1, 1 );
  plasma = zeros ( n + 1, 1 );
%
%  Initialize variables;
%
  t(1) = t0;
  intake(1) = y0(1);
  medicine_in_intestines(1) = y0(2);
  plasma(1) = y0(3);
%
%  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) + dt;

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

    medicine_in_intestines(j+1) = medicine_in_intestines(j) ...
      + intake(j) ...
      - dt * absorption_rate * medicine_in_intestines(j);

    plasma(j+1) = plasma(j) ...
      + dt * absorption_rate * medicine_in_intestines(j) ...
      - dt * excretion_rate * plasma(j);

  end
%   
%  Plot the intake.
%
  clf ( );
  hold ( 'on' );
    plot ( t, intake, 'b-', 'linewidth', 3 );
    plot ( [ t(1),t(n+1)], [0.0, 0.0], 'k-', 'linewidth', 3 );
  hold ( 'off' )
  grid ( 'on' )
  title ( 'Intake', 'Fontsize', 24 );
  xlabel ( '<-- Hours -->', 'Fontsize', 16 );
  ylabel ( '<-- Level (mg) -->', 'Fontsize', 16 );
  filename = 'dosage_intake.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Convert plasma level to plasma concentration.
%
  blood_volume = 4.6;
  plasma = plasma / blood_volume;
%   
%  Plot the intestinal level.
%
  clf ( );
  hold ( 'on' );
    plot ( t, medicine_in_intestines, 'b-', 'linewidth', 3 );
    plot ( [ t(1),t(n+1)], [0.0, 0.0], 'k-', 'linewidth', 3 );
  hold ( 'off' )
  grid ( 'on' )
  title ( 'Intestinal drug level.', 'Fontsize', 24 );
  xlabel ( '<-- Hours -->', 'Fontsize', 16 );
  ylabel ( '<-- Level (mg) -->', 'Fontsize', 16 );
  filename = 'dosage_intestines.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%   
%  Plot the bloodstream concentration.
%
  medicinal_level = 800.0;
  toxic_level = 1000.0;
  clf ( );
  hold ( 'on' );
    plot ( t, plasma, '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 ( 'Concentration',  'Therapeutic Concentration', 'Toxic Concentration', ...
    'Location', 'South' );
  grid ( 'on' )
  title ( 'Bloodstream medicinal concentration.', 'Fontsize', 24 );
  xlabel ( '<-- Hours -->', 'Fontsize', 16 );
  ylabel ( '<-- Concentration (mg/l) -->', 'Fontsize', 16 );
  ylim ( [ 0.0, 1100.0 ] );
  filename = 'dosage_bloodstream.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end
