4619.sagews

# This calculates the level of coffee effecting you in units of 1 coffee. Here is an example of a coffee at 9am, 10am and 1pm. This uses an 1 hour exponential increase until peak blood level, then a 5.7 hour half life. Scroll down to see the graph. This results in about 0.9 coffees still in the system at 10pm.

# change these values
weight=80 #in kg
dose=[1,1,1] # Doses of coffee. A shot of coffee has (~55mg caffeine)
time=[9,10,13] # array of times after t0 of doses. Time in hours since 1st dose
t_half=5.7 # biological half life 
t_absorb=1 # time to peak blood level

# this makes a exponential (exp) model of blood levels
# 1-exp(-x) model of absorption and a exp(-x) model elimination
# this is the same as the nuci number of a duaghter radionuclide in radioactive decay:
# Nd(t)=Np(0)* lambdap/(lambdad-lambdap)* (exp(-lambdap*t)-exp(-lambdad*t)

# dont change these values
#dose=[n/0.3*8/weight for n in dose] # 0.3mg/kg > 8ng/mLdose
lambd=log(2)/t_half # the mean life -the deacy constant
absorb=log(2)/(t_absorb/2) # as it takes 2 half lives to ~peak, so decay constant is 4.5*times for absorbtion

# loop over each dose, adding a function for each one independantly.
var('t')
f=0
for x in range(len(dose)):
    f=f+ heaviside(t-time[x])* dose[x] * absorb/(absorb-lambd) * ( exp(-lambd*(t-time[x]))-exp(-absorb*(t-time[x])) )
    
# Now plot the graph
a=plot(f,time[0],ceil(time[len(time)-1]*2),gridlines=true) # plot function from 0 to twice highest dosing time
a.axes_labels(['hours','coffees'])
a.show(gridlines="minor",figsize=[10,5])

# References
# Caffeine half life is 5.7 hours 
# Source :http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&list_uids=7361718&dopt=Abstract
# Peak blood concentration is at 1 hour
# Source: http://www.fass.se/LIF/produktfakta/artikel_produkt.jsp?NplID=20030808000012&DocTypeID=6#pharmacokinetic
# A typical expresso coffee has 40-75 mg