```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

data = np.fromfile("fit_result_sedimenting_n_1p595r_1p45fix_0p0535_wav_532.dat", dtype=np.float64)
data = data.reshape(len(data)//10,10)

x = data[:,0]*0.0535
y = data[:,1]*0.0535
z = data[:,2]*0.0535
```


```python
import numpy as np
import matplotlib.pyplot as plt

### Just some matplotlib tweaks
import matplotlib as mpl

mpl.rcParams["xtick.direction"] = "in"
mpl.rcParams["ytick.direction"] = "in"
mpl.rcParams["lines.markeredgecolor"] = "k"
mpl.rcParams["lines.markeredgewidth"] = 1
mpl.rcParams["figure.dpi"] = 200
from matplotlib import rc
rc('font', family='serif')
rc('text', usetex=True)
rc('xtick', labelsize='medium')
rc('ytick', labelsize='medium')
rc("axes", labelsize = "large")
def cm2inch(value):
    return value/2.54
```


```python
def traj_sedim(x, a, b):
    return a * x + b
    
```


```python
# We take the data where we are sure that the particle is far enough from the wall - and from the 
z = z[0:6000]
time = np.arange(0,len(z)) / 100

```


```python
figsize=( cm2inch(16),cm2inch(8))

plt.plot(time, z)
plt.xlabel("$t$ (s)")
plt.ylabel("$z$ ($\\mu$m)")
plt.tight_layout(pad=0.1)
plt.savefig("traj_sedimenting.pdf")
```


```python
popt, pcov = curve_fit(traj_sedim, time, z)
```


```python
plt.plot(time, z)
plt.plot(time, traj_sedim(time, *popt))
plt.xlabel("$t$ (s)")
plt.ylabel("$z$ ($\\mu$m)")
```


```python
popt
```




    array([-0.26669228, 30.0528722 ])



The sedimentation drift velocity for a spherical particle in bulk writes:

\begin{equation}
v_\mathrm{sed} = \frac{2}{9} \frac{\Delta \rho a^2 g}{\eta}
\end{equation}


```python
a = 1.45e-6
v_th = 2/9 * 50 * a ** 2 * 9.81 / 0.001 * 1e6
```


```python
v_th
```




    0.22917250000000006




```python
def msd(x, Dt):
	"""Function that returns the MSD for a list of time indices Dt for a trajectory x"""
	_msd = lambda x, Dt: np.mean((x[:-Dt] - x[Dt:]) ** 2)
	return [_msd(x, i) for i in Dt]
```


```python
Dt = np.arange(1,1000,5)
```


```python
MSD = msd(z, Dt)

```


```python
def msd_th(x, d, v, n):
    return 2*d*x + v**2 * x**2 + n**2

def msd_log(x, d, v, n):
    return np.log( 2*d*x + v**2 * x**2 + n**2)
```


```python
plt.plot(np.log(Dt/100),np.log(MSD), "o")
```


```python
popt, pcov = curve_fit(msd_log, Dt/100, np.log(MSD), p0 = [0.13, 0.22, 4e-2], maxfev=2000, ftol=1e-2)
```


```python
popt
```




    array([0.14092578, 0.22348519, 0.04612557])




```python
perr = np.sqrt(np.diag(pcov))
```


```python
perr
```




    array([0.0007985 , 0.0008066 , 0.00149602])




```python
Dt = np.array([*np.arange(1,10,1), *np.arange(10,100,10), *np.arange(100,1000,100)])
MSD = msd(z, Dt)

```


```python
figsize=( cm2inch(16),cm2inch(8))
plt.loglog(Dt/100,MSD, "o")
t_th=np.linspace(0.7e-2, 20,10000)
plt.plot(t_th, msd_th(t_th, *popt), color = "k", linewidth=2)
plt.xlim((0.8e-2, 15))
plt.ylim((None, 10))
plt.xlabel("$\Delta t$ (s)")
plt.ylabel("MSD ($\mathrm{\mu m^2}$)")
plt.tight_layout(pad=0.1)
plt.savefig("msd_sedimentation.pdf")
```


```python
relat_err_D = abs(0.144 - 0.146)/0.146
relat_err_D*100
```




    1.3698630136986314




```python
relat_err_v = abs(222 - 229)/229
relat_err_v*100
```




    3.056768558951965




```python
popt[1]
```




    0.2217820942226398




```python
popt[0]
```




    0.14377799719977521




```python
9/2 * 222e-9 * 0.001 / (1.45e-6 ** 2) / 9.81 
```




    48.43513074212656




```python
(50 - 48.5) / 50 * 100
```




    3.0




```python

```

fae45a4835298df4aa2ed78a1c67458b859177fd
ipynb
Jupyter Notebook
02_body/chapter2/images/sedimenting/Sedimenting_particle.ipynb
eXpensia/Confined-Brownian-Motion
bd0eb6dea929727ea081dae060a7d1aa32efafd1
MIT
02_body/chapter2/images/sedimenting/Sedimenting_particle.ipynb
eXpensia/Confined-Brownian-Motion
bd0eb6dea929727ea081dae060a7d1aa32efafd1
MIT
02_body/chapter2/images/sedimenting/Sedimenting_particle.ipynb
eXpensia/Confined-Brownian-Motion
bd0eb6dea929727ea081dae060a7d1aa32efafd1
MIT
Qwen/Qwen-72B
1. YES
2. YES