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

def simulation( doses_per_day, dosage_per_day, blood_volume):
    # 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 = 2
    # dosage_per_day = 12000  #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 = round(i + steps_between_doses, 0)
        else:
            intake[i] = 0
            
    #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]

        plasma_level[i+1] = plasma_level[i] - excretion[i] + absorption[i]

        plasma_concentration[i+1] = plasma_level[i+1] / blood_volume

        medicine_in_intestines[i+1] = medicine_in_intestines[i] - absorption[i] + intake[i]
        

    # After the simulation has been run check to see that the last doseage
    # interval is in the safety region
    if all(test > medicinal_level  and test < toxic_level 
            for test in plasma_concentration[-1*int(steps_between_doses):-1]):
        return True
    else:
        return False

def ask_Weight():

    # Ask User for weight
    Weight = input('Enter in your weight (kg): ')

    # Calculate blood volume from weight
    # BloodVol = Weight * AvgBloodVolperKg
    # from http://reference.medscape.com/calculator/estimated-blood-volume
    AvgBloodVolperKg = 5.75
    BloodVol = 1e-2*Weight * AvgBloodVolperKg

    mg = []

    # Loop through the dose levels for a single blood volume
    for j in np.linspace(1000, 20000, 100):
        check = simulation(4, j, BloodVol)
        if check:
            mg.append(j)
    print 'For your weight of {}kg the total daily dose {:.1f}mg assuming 4 doses a day'.format(Weight, np.mean(mg))

ask_Weight()