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

def drug_dosage_sensitivity(e_rate):
	# 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
	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(float(e_rate)) / 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

	for i in range(N):
		# Otherwise, run the model
		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
		    
	#TODO: Create another loop that calculates the absorption, excretion,
	#plasma level, and plasma concentration over time
	return plasma_concentration