import matplotlib.pyplot as plt
import math
# drug_dosage_start.py
#
# VARIABLE DICTIONARY DRUG DOSAGE MODEL:
# --------------------------------------------------------------------
#     Drug levels in the blood, for certain drugs, can be modeled by a pair of
#     coupled first-order linear differential equations.  For compactness, let
#     A= medicine_in_intestines and B = plasma_level.  Variable names are defined
#     below.
#     dA/dt = -absorption_rate*A + intake
#     dB/dt = -excretion_rate*B + absorption_rate*A

# Model parameters
time_step = 0.25       # simulation time step in hours
start_time = 0         # in hours
end_time = 48          # in hours
doses_per_day = 4.0
dosage_per_day = 11800  #mg of the drug
absorption_rate = 0.25 # absorption rate constant
half_life = 6          # in hours
blood_volume = 4.6     # in liters
medicinal_level = 800  # level at which drug becomes effective mg/l
toxic_level = 1000     # level at which drug becomes toxic  mg/l

# Derived constants
N = int((end_time - start_time) / time_step)   # number of simulation steps
dosage = dosage_per_day / doses_per_day        # amount of each dose
dosage_interval = 24 / doses_per_day           # in hours
excretion_rate = math.log(2) / half_life
steps_between_doses = dosage_interval / time_step

# Time-varying quantities, arrays with one value per time step
# The syntax, "[0]*(N+1)," creates a one-dimensional array of length N+1 whose values are all zero.

t = [0]*(N+1)                       # time in hours
intake = [0]*(N+1)                  # dose at this time step
medicine_in_intestines = [0]*(N+1)  # level of drug in intestines
plasma_level = [0]*(N+1)            # level of drug in the blood
plasma_concentration = [0]*(N+1)    # concentration of drug in blood
absorption = [0]*(N+1)              # amount absorbed at this time step
excretion = [0]*(N+1)               # amount excreted at this time step

# Initialize variables - assume all others start at 0
t[0] = start_time
step_next_dose = 1         # give first dose at start time

# Loop to calculate intake of drug over time 
# 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. Number of doses must be an
# integer divisor of 24 for this to work.  Otherwise change the time step
for i in range(N):
    t[i+1] = t[i] + time_step
    # Is it time for a dose?
    if i == step_next_dose:
        intake[i] = dosage
        step_next_dose = i + steps_between_doses
    else:
        intake[i] = 0

# Loop over all of the timesteps
for i in range(N):
	# Compute the absorption and excretion levels at current timestep
	absorption[i] = time_step*absorption_rate*medicine_in_intestines[i]
	excretion[i] = time_step * excretion_rate * plasma_level[i]
	# Compute the amount of the drug in the intestines at the next timestep
	medicine_in_intestines[i+1] = medicine_in_intestines[i] - absorption[i] + intake[i]
	# Compute the plasma_level and plasma_concentration levels at next timestep
	plasma_level[i+1] = plasma_level[i] - excretion[i] + absorption[i]
	plasma_concentration[i+1] = plasma_level[i+1]/blood_volume
        
#TODO: Create another loop that calculates the absorption, excretion,
#plasma level, and plasma concentration over time

  
    
    
# Now Plot the results
m_level = [medicinal_level]*(N+1)
t_level = [toxic_level]*(N+1)

plt.plot(t,plasma_concentration, color = 'b', label= "Plasma Concentration")
plt.plot(t,m_level, color = 'g', label = "Medicinal Level")
plt.plot(t,t_level, color = 'r', label = "Toxic Level")
plt.title('Pharmacokinetic Model')
plt.xlabel('Time in Hours')
plt.ylabel('Concentration')


plt.legend(loc=4)

plt.minorticks_on()
plt.ylim(0,1100)
plt.show()