import math
import numpy as np
import matplotlib.pyplot as plt
import pylab

# 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*2          # in hours
doses_per_day = 128
dosage_per_day = 10000  #mg of the drug
absorption_rate = 0.25 # absorption rate constant
half_life = 12        # 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
#ASSUMPTION CORRECTION
#dosage_interval was changed to account for 16 hour awake day
dosage_interval = 16 / 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 = 0         # 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(0,N):
    t[i+1] = t[i] + time_step
    # Is it time for a dose?
    if i >= step_next_dose and t[i]%24 <= 16:
        intake[i] = dosage
	#the step_next_dose is rounded in order to allow for a doses_per_day that is not a divisor of 24
        step_next_dose = round(i + steps_between_doses,0)
    else:
        intake[i] = 0
        
#TODO: Create another loop that calculates the absorption, excretion,
#plasma level, and plasma concentration over time

for i in range(N):
	absorption[i] = time_step * absorption_rate *medicine_in_intestines[i] 
	excretion[i] = time_step * excretion_rate * plasma_level[i]
	medicine_in_intestines[i+1] = medicine_in_intestines[i] - absorption[i] + intake[i]
	plasma_level[i+1] = plasma_level[i] - excretion[i] +absorption[i]
	plasma_concentration[i+1] = plasma_level[i+1] /blood_volume
  
    
    
# Now Plot the results

plt.plot(t,plasma_concentration, color = 'b', label= "Plasma Concentration")

m_level = [medicinal_level]*(N+1)
plt.plot(t, m_level, color= 'r', label= "Medicinal Level",)

t_level = [toxic_level]*(N+1)
plt.plot(t, t_level, color = 'g', label = "Toxic Level")

plt.legend(loc=4)

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