ESBS

Biophysics practicals 2012

This SAGE notebook contains some hints and programs that will allow you to analyse the data recorded during the experiments. Exemples are provided to show you the more appropriate way to analyse and display your data. The graphics could be exported directly to fit within your reports as png files.

1. TP 2 Mass transport : Non-linear fitting of the time dependant evolution of the gel's temperature

The time evolution of the gel's temperature could be modelled using the following equation:

$T(t)=A+B\left(1-e^{-Ct}\right)$

where $A=T_{0}, B=\frac{P}{aC} and C=a$

First : Get an idea of the effect of each of the parameters on the curve using the following interactive plot

In particular, check which parameter influences 1) the value of the plateau 2) the initial slope

@interact
def _(A=(10..100),B=(1..10),CI=(1..20)):
    C=CI*0.001
    f(t) = A + B*(1-exp(-C*t))
    pf = plot(f,0,1000,axes_labels=['Time (s)','Temp (Deg C)'])
    show(pf)

Second : We should define the function to be fitted, in a way that will be recognized by NUMPY

from numpy import *
def temp(time,pA,pB,pC):
    """ Equation describing heat exchange within an electrophoresis gel"""
    y= pA+pB*(1 - exp(-time*pC))
    return y

Third : The measured values have to be entered in the program. First option : enter the values manually as shown below

Times = [10, 30, 100, 300, 400, 600, 800, 1200, 2000]
Texp =  [20.2, 20.5, 21.3, 22.9, 23.3, 23.7, 23.9, 24.0, 24.0]

Second option: download a text file named arduino.dat, that contains the recorded temperatures according the following format:

Fri Mar 16 17:31:34 2012 23.56
Fri Mar 16 17:31:40 2012 23.5

And run the following script

...

# arduino prep
# Read Arduino data 
# B. Kieffer Mar 2012

import time 

# 1. Read data into a list all_line
fname = DATA+'arduino.dat'
myfile = open(fname,'r')
all_lines = myfile.readlines()
myfile.close()

txtime = []
Texp = []
Times = []
# 2. Extract relevent informations from all_line
for line in all_lines:
	line_field = line.split()
	time_field = line_field[3]
	val = line_field[5]
	txtime.append(time_field)
	Texp.append(float(val))

# 3. Build a vector of time delays starting from the first point
first = True	
for t in txtime :
   t1 = time.strptime(t, "%H:%M:%S")
   t_sec = t1.tm_hour*3600+t1.tm_min*60+t1.tm_sec
   if first :
      Times.append(0)
      t0 = t_sec # First point is the time 0
      first = False
   else :
      Times.append(t_sec-t0)
print " %i lines read from file arduino.dat " % len(Times)

61 lines read from file arduino.dat

print Times

[0, 31, 60, 91, 120, 150, 180, 210, 240, 270, 301, 330, 361, 390, 421, 450, 481, 510, 541, 570, 601, 630, 661, 690, 721, 750, 780, 810, 840, 870, 900, 931, 960, 991, 1020, 1051, 1080, 1111, 1140, 1171, 1200, 1231, 1260, 1291, 1320, 1351, 1380, 1410, 1440, 1470, 1500, 1530, 1560, 1590, 1621, 1650, 1681, 1710, 1741, 1770, 1801]

Finally : Perform the non-linear fit using the program curve_fit that is included in the optimize library of SCIPY

from scipy import optimize
# 1. Initialize parameters
Pini = [18., 2., 0.001] # Initial guess of the parameters
Pow = 10.0 # Power in the gel

# 2. Perform the fit
PP, pcov = optimize.curve_fit(temp, array(Times), array(Texp), Pini)

# 3. Report and plot the results
Cap = Pow/(PP[1]*PP[2])
print "Optimized parameters:"
print "T0: %f (Deg C) " % PP[0]
print "Calorific capacity: %f (J.K-1)  Rate: %f (s-1)"%(Cap,PP[2])
a = scatter_plot(zip(Times,Texp))
a += plot(temp(x,PP[0],PP[1],PP[2]), 0,max(Times),axes_labels=['Time (s)','Temp (Deg C)'])
show(a)

Optimized parameters:
T0: 20.059797 (Deg C) 
Calorific capacity: 1141.509034 (J.K-1)  Rate: 0.001589 (s-1)