# Import Statements
import pylab as pl	# PyLab is Numpy, Scipy and Matlplotlib
import numpy as np	# Numpy explicitly used for arrays
from drug_dosage_sensitivity import *	# This imports the function made from the base code
# This is the code for running a sensitivity analysis on the
# excretion rate from the base code.  Changing this parameter
# has the effect of pulling the steady state up or down due to the
# changes in the frequency of excretion of the drug from the body 
# This code utilizes numpy arrays and plotting to show the results
# of the sensitivity analysis.  The first loop generates a 2D array
# from plasma_concentration values from base code, which was turned
# into a large function for the purposes of the sensitivity analysis

## Set up the ranges for sensitivity
er = np.arange(1,5,0.5)	# An array of excretion rates
t = np.arange(0,48,0.25)	# An array holding the time steps
sens_conc = np.zeros([len(er),len(t)+1])	# A 2D array for plasma_concentrations at each excretion rate
half_life = 6	# The half-life of the drug in the body

## Run the Sensitivity Analysis
# This loops over the values in the excretion_rate array
for i in range(len(er)):
	# Pulls plasma_concentration for each excretion rate and adds them to a row in the 2D array
	sens_conc[i,:] = drug_dosage_sensitivity(er[i])

# This loops through the excretion rates and plots it in the time domain against the corresponding
# plamsa_concentration.  There is a legend for excretion rate.
for i in range(len(er)):
	pl.plot(t, sens_conc[i,0:192], label='Excretion Rate: %.2f' % ((math.log(er[i]) / half_life)))
	pl.legend(loc=4)	# Sets the legend at the bottom left hand corner


N = len(t)-1		   # Ensure that all arrays are the same length
medicinal_level = 800  # level at which drug becomes effective mg/l
toxic_level = 1000     # level at which drug becomes toxic  mg/l
m_level = [medicinal_level]*(N+1)	# Create an array of medicinal levels
t_level = [toxic_level]*(N+1)		# Create an array of toxic levels
# Plot the medicinal and toxic levels
plt.plot(t,m_level, 'og', label = "Medicinal Level")
plt.plot(t,t_level, 'or', label = "Toxic Level")
# Set title and axis labels
plt.title('Sensitivity Plot')
plt.xlabel('Time in Hours')
plt.ylabel('Concentration')
# Set plotting parameters
plt.minorticks_on()
plt.ylim(0,1100)
# Only show once to ensure all plots are on the same graph
pl.show()