In [8]:
%pylab inline
import sys
sys.path.append(PATH_TO_CONTRACT_PY) ## Put the path to contract.py here !!!!
import contract as contract
from scipy.integrate import odeint
from scipy.optimize import curve_fit
##set parameters of contract simulation + initialize grid
ngp=64
contract.ngp=ngp
contract.setup_grid()
t=linspace(0,80.,2000)
In [9]:
##redefine compute function with more flexible parameter input
def compute(a,b,r0,w_init):
l=sqrt(a/b)
s0=-1./sqrt(a*b)
#set parameters
contract.l=l
contract.s0=s0
contract.rho0=1.
contract.rho_init=r0*ones(ngp)
contract.w_init=w_init
return odeint(contract.dt_packed,contract.init_packed(),t)
##define fit function for time profile of contractions
def fit_func(t,aa,tau):
return aa*(1.-exp(-t/tau))
##define functional form of time scaling
def ftaus(_ww,_aa,_bb):
return _aa+_bb*_ww**2
In [13]:
## choose a= eta/s rho0^2 and b= gamma/s rho0^2 around which to determine alpha beta
## (taken from fit to a typical contraction profile)
a,b,r0=array([ 2.26413429e+00, 7.09505347e-05, 3.26194452e-01]) #
a=a/3.
b=b/3.
b=b/4.
## define interval over which to sample
_w=(arange(5)+1)*400
aa=linspace(a/2.,a*2.,10.)
bb=linspace(b/2.,b*2.,10.)
afacs=[]
bfacs=[]
## this loop determines output time scales for many parameters
for _a in aa:
for _b in bb:
taus=[]
t=linspace(0,80.,2000)
for w_init in _w:
o=compute(_a,_b,r0,w_init)
w=o[:,0]
w0=-(w-w_init)/w_init
cf=curve_fit(fit_func,t,w0)
a,tau=cf[0]
taus.append(tau)
cftau=curve_fit(ftaus,_w,taus)
afacs.append(cftau[0][0]/_a)
bfacs.append(cftau[0][1]/_b)
In [30]:
##print the alpha and beta determined from the sampling
alpha=mean(afacs)
beta=mean(bfacs)
alpha_err=sqrt(var(afacs))
beta_err=sqrt(var(bfacs))
print("alpha=%0.2e +- %0.2e"%(alpha,alpha_err))
print("beta=%0.2e +- %0.2e"%(beta,beta_err))