# -*- coding: utf-8 -*-
# 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
import math
import numpy as np
import matplotlib.pyplot as plt
import pylab

# Model parameters
time_step = 0.25       # simulation time step in hours
start_time = 0         # in hours
end_time = 48          # in hours
doses_per_day = 2
#Start dosage_per_day below medicinal level
dosage_per_day = 10000  #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

plasma_concentration = [0]*(N+1)    # concentration of drug in blood
min_not_found = True
#Loop until we reach one below the toxic level
while (max(plasma_concentration) < (toxic_level - 1)):
	
	if max(plasma_concentration) > medicinal_level and min_not_found:
		min_plasma_concentration = plasma_concentration
		min_dosage_per_day = dosage_per_day
		min_not_found = False

	#Increase the dosage_per_day each time through the loop
	dosage_per_day = dosage_per_day + 1  #mg of the drug

	# 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
	    #For time released dosage, we make intake constant
	    intake[i] = dosage_per_day*time_step/24
		
	#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
	    
print 'Minimum concentration found at a ', min_dosage_per_day, ' mg dosage a day.'
print 'Maximum concentration found at a ', dosage_per_day, ' mg dosage a day.'
# Now Plot the results

plt.plot(t,plasma_concentration, color = 'b', label= "Maximum Plasma Concentration")
plt.plot(t,min_plasma_concentration, color = 'y', label= "Minimum 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()