In [1]:
import numpy as np
from scipy.optimize import leastsq
from matplotlib.pylab import *
from matplotlib.colors import LogNorm
%matplotlib inline
Define all functions necessary for the fit
In [2]:
LogDDM = lambda p, dts: np.log(p[0] * (1-np.exp(-dts/p[2])) + p[1])
Read the tables of $\mathcal{D}(\Delta t, q)$ and $\Delta t$ from files. Set the pixel size and thus the scale of wavenumbers $q$. Set the longest time index fit.
In [38]:
DDM = np.load('DDM_Colloid.npy')
dts = np.load('dt_Colloid.npy')
pixelSize = 6.450/10. #in micrometre
qs = 2*np.pi/(2*DDM.shape[-1]*pixelSize) * np.arange(DDM.shape[-1])
tmax = -6
Perform the fit at each $q$. Generate the values of the fitted function
In [39]:
params = np.zeros((DDM.shape[-1], 3))
matrixfit = np.zeros(DDM[:tmax].T.shape)
for iq, ddm in enumerate(DDM[:tmax].T):
params[iq] = leastsq(
#function to minimize
lambda p, dts, logd: LogDDM(p, dts) - logd,
#initial parameters
[ddm.ptp(), ddm.min(), 1],
#data on which to perform minimization
args=(dts[:tmax], np.log(ddm))
)[0]
matrixfit[iq] = np.exp(LogDDM(params[iq], dts[:tmax]))
Show parameters $A(q)$ and $B(q)$ and deduce the range of $q$ where the fit is correct
In [118]:
figure(figsize=(8,6))
plot(qs, params[:,0], 'o')
plot(qs, params[:,1], 'o')
xscale('log')
yscale('log')
xlabel(r'$q\,(\mu m^{-1})$')
ylabel(r'$A(q),\, B(q)$')
iqmin = np.argmax(params[:,0])
iqmax = np.where(params[:,0]<10*params[:,1])[0][0]
axvspan(qs[iqmin], qs[iqmax], color=(0.9,0.9,0.9))
print iqmin, iqmax
Show the results of the fit function of $q$ for various $\Delta t$ and function of $\Delta t$ for various $q$.
In [119]:
figure(figsize=(16,6))
ax1 = subplot(1,2,1)
for i in range(0,len(dts[:tmax]),5):
plot(qs, DDM[i]/512**2, 'o', color=cm.autumn(i/float(len(dts[:tmax]))))
plot(qs, matrixfit[:,i]/512**2, '-k')
xscale('log')
yscale('log')
ylabel(r'$\mathcal{D}$')
xlabel(r'$q\,(\mu m^{-1})$')
axvspan(qs[iqmin], qs[iqmax], color=(0.9,0.9,0.9))
ax2 = subplot(1,2,2, sharey=ax1)
for i,iq in enumerate(np.logspace(np.log10(iqmin),np.log10(iqmax),10, base=10).astype(int)):
plot(dts[:tmax], DDM[:tmax,iq]/512**2, 'o', color=cm.autumn(i/10.))
plot(dts[:tmax], matrixfit[iq]/512**2, '-k')
xscale('log')
xlabel(r'$\Delta t\,(s)$')
setp(ax2.get_yticklabels(), visible=False)
ylim(1, 1e4)
Out[119]:
The corresponding ISF for various $q$ function of $\Delta t$ or, to show the Brownian rescaling, function of $q^2\Delta t$.
In [120]:
figure(figsize=(16,6))
ax1 = subplot(1,2,1)
for i,iq in enumerate(np.logspace(np.log10(iqmin),np.log10(iqmax),10, base=10).astype(int)):
plot(dts[:tmax], 1- (DDM[:tmax,iq]-params[iq,1])/params[iq,0], 'o', color=cm.autumn(i/10.))
plot(dts[:tmax], np.exp(-dts[:tmax]/params[iq,2]), '-k')
xscale('log')
ylabel(r'$f(q,\Delta t)$')
xlabel(r'$\Delta t\,(s)$')
ax2 = subplot(1,2,2, sharey=ax1)
for i,iq in enumerate(np.logspace(np.log10(iqmin),np.log10(iqmax),10, base=10).astype(int)):
plot(qs[iq]**2*dts[:tmax], 1- (DDM[:tmax,iq]-params[iq,1])/params[iq,0], 'o', color=cm.autumn(i/10.))
xscale('log')
xlabel(r'$q^2\Delta t\,(s/\mu m^2)$')
setp(ax2.get_yticklabels(), visible=False)
ylim(-0.1,1.1)
Out[120]:
Show the fit parameter $\tau_d(q)$ and fit it to obtain the diffusion coefficient $D$.
In [121]:
figure(figsize=(8,6))
plot(qs, params[:,2], 'o')
xscale('log')
yscale('log')
xlabel(r'$q\,(\mu m^{-1})$')
ylabel(r'$A(q),\, B(q)$')
axvspan(qs[iqmin], qs[iqmax], color=(0.9,0.9,0.9))
D = np.exp(-leastsq(
lambda p, q, td: p[0] - 2*np.log(q) - np.log(td),
[-1.],
args=(qs[iqmin:iqmax], params[iqmin:iqmax,2])
)[0][0])
plot([0.1, 4], 1/(D*np.array([0.1, 4])**2), '-r')
print u'D = {:.02f} µm²/s'.format(D)
In [ ]: