Code examples¶
Applications¶
Creep compliance analysis¶
This example uses the contour data of an HL60 cell in the OS to compute its GGF and creep compliance. The contour data were determined from this phase-contrast video (prior to video compression). During stretching, the total laser power was increased from 0.2W to 1.3W (reflexes due to second harmonic effects appear as white spots).
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 | import ggf
import h5py
import lmfit
import matplotlib.pylab as plt
import numpy as np
def ellipse_fit(radius, theta):
'''Fit a centered ellipse to data in polar coordinates
Parameters
----------
radius: 1d ndarray
radial coordinates
theta: 1d ndarray
angular coordinates [rad]
Returns
-------
a, b: floats
semi-axes of the ellipse; `a` is aligned with theta=0.
'''
def residuals(params, radius, theta):
a = params["a"].value
b = params["b"].value
r = a*b / np.sqrt(a**2 * np.sin(theta)**2 + b**2 * np.cos(theta)**2)
return r - radius
parms = lmfit.Parameters()
parms.add(name="a", value=radius.mean())
parms.add(name="b", value=radius.mean())
res = lmfit.minimize(residuals, parms, args=(radius, theta))
return res.params["a"].value, res.params["b"].value
# load the contour data (stored in polar coordinates)
with h5py.File("data/creep_compliance_data.h5", "r") as h5:
radius = h5["radius"].value * 1e-6 # [µm] to [m]
theta = h5["theta"].value
time = h5["time"].value
meta = dict(h5.attrs)
factors = np.zeros(len(radius), dtype=float)
semimaj = np.zeros(len(radius), dtype=float)
semimin = np.zeros(len(radius), dtype=float)
strains = np.zeros(len(radius), dtype=float)
complnc = np.zeros(len(radius), dtype=float)
for ii in range(len(radius)):
# determine semi-major and semi-minor axes
smaj, smin = ellipse_fit(radius[ii], theta[ii])
semimaj[ii] = smaj
semimin[ii] = smin
# compute GGF
if (time[ii] > meta["time_stretch_begin [s]"]
and time[ii] < meta["time_stretch_end [s]"]):
power_per_fiber = meta["power_per_fiber_stretch [W]"]
f = ggf.get_ggf(
model="boyde2009",
semi_major=smaj,
semi_minor=smin,
object_index=meta["object_index"],
medium_index=meta["medium_index"],
effective_fiber_distance=meta["effective_fiber_distance [m]"],
mode_field_diameter=meta["mode_field_diameter [m]"],
power_per_fiber=power_per_fiber,
wavelength=meta["wavelength [m]"],
poisson_ratio=.5)
else:
power_per_fiber = meta["power_per_fiber_trap [W]"]
f = np.nan
factors[ii] = f
# compute compliance
strains = (semimaj-semimaj[0]) / semimaj[0]
complnc = strains / factors
compl_ival = (time > meta["time_stretch_begin [s]"]) * \
(time < meta["time_stretch_end [s]"])
stretch_index = np.where(compl_ival)[0][0]
complnc_1 = strains/factors[stretch_index]
# plots
plt.figure(figsize=(8, 7))
ax1 = plt.subplot(221, title="ellipse fit semi-axes")
ax1.plot(time, semimaj*1e6, label="semi-major axis")
ax1.plot(time, semimin*1e6, label="semi-minor axis")
ax1.legend()
ax1.set_xlabel("time [s]")
ax1.set_ylabel("axis radius [µm]")
ax2 = plt.subplot(222, title="GGF")
ax2.plot(time, factors)
ax2.set_xlabel("time [s]")
ax2.set_ylabel("global geometric factor [Pa]")
ax3 = plt.subplot(223, title="strain")
ax3.plot(time, (strains)*100)
ax3.set_xlabel("time [s]")
ax3.set_ylabel("deformation $w(t)/r_0$ [%]")
ax4 = plt.subplot(224, title="creep compliance")
ax4.plot(time[compl_ival], complnc[compl_ival])
ax4.set_xlabel("time [s]")
ax4.set_ylabel("compliance $J(t)$ [Pa⁻¹]")
for ax in [ax1, ax2, ax3, ax4]:
ax.set_xlim(0, np.round(time.max()))
ax.axvline(x=meta["time_stretch_begin [s]"], c="r")
ax.axvline(x=meta["time_stretch_end [s]"], c="r")
plt.tight_layout()
plt.show()
|
Reproduction tests¶
Radial stresses of a prolate spheroid¶
This examples computes radial stress profiles for spheroidal objects in the optical stretcher, reproducing figures (9) and (10) of [BCG09].
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 | import matplotlib.pylab as plt
import numpy as np
import percache
from ggf.stress.boyde2009.core import stress
@percache.Cache("stress_reproduced.cache", livesync=True)
def compute(**kwargs):
"Locally cached version of ggf.core.stress"
return stress(**kwargs)
# variables from the publication
alpha = 47
wavelength = 1064e-9
radius = alpha * wavelength / (2 * np.pi)
kwargs = {"stretch_ratio": .1,
"object_index": 1.375,
"medium_index": 1.335,
"wavelength": wavelength,
"beam_waist": 3,
"semi_minor": radius,
"power_left": 1,
"power_right": 1,
"poisson_ratio": 0,
"n_points": 200,
}
kwargs1 = kwargs.copy()
kwargs1["power_right"] = 0
kwargs1["stretch_ratio"] = 0
kwargs1["dist"] = 90e-6
kwargs2 = kwargs.copy()
kwargs2["power_right"] = 0
kwargs2["stretch_ratio"] = .05
kwargs2["dist"] = 90e-6
kwargs3 = kwargs.copy()
kwargs3["power_right"] = 0
kwargs3["stretch_ratio"] = .1
kwargs3["dist"] = 90e-6
kwargs4 = kwargs.copy()
kwargs4["dist"] = 60e-6
kwargs5 = kwargs.copy()
kwargs5["dist"] = 120e-6
kwargs6 = kwargs.copy()
kwargs6["dist"] = 200e-6
# polar plots
plt.figure(figsize=(8, 5))
th1, sigma1 = compute(**kwargs1)
ax1 = plt.subplot(231, projection='polar')
ax1.plot(th1, sigma1, "k")
ax1.plot(th1 + np.pi, sigma1[::-1], "k")
th2, sigma2 = compute(**kwargs2)
ax2 = plt.subplot(232, projection='polar')
ax2.plot(th2, sigma2, "k")
ax2.plot(th2 + np.pi, sigma2[::-1], "k")
th3, sigma3 = compute(**kwargs3)
ax3 = plt.subplot(233, projection='polar')
ax3.plot(th3, sigma3, "k")
ax3.plot(th3 + np.pi, sigma3[::-1], "k")
for ax in [ax1, ax2, ax3]:
ax.set_rticks([0, 1.5, 3, 4.5])
ax.set_rlim(0, 4.5)
th4, sigma4 = compute(**kwargs4)
ax4 = plt.subplot(234, projection='polar')
ax4.plot(th4, sigma4, "k")
ax4.plot(th4 + np.pi, sigma4[::-1], "k")
ax4.set_rticks([0, 4, 8, 12])
ax4.set_rlim(0, 12)
th5, sigma5 = compute(**kwargs5)
ax5 = plt.subplot(235, projection='polar')
ax5.plot(th5, sigma5, "k")
ax5.plot(th5 + np.pi, sigma5[::-1], "k")
ax5.set_rticks([0, 1.5, 3, 4.5])
ax5.set_rlim(0, 4.5)
th6, sigma6 = compute(**kwargs6)
ax6 = plt.subplot(236, projection='polar')
ax6.plot(th6, sigma6, "k")
ax6.plot(th6 + np.pi, sigma6[::-1], "k")
ax6.set_rticks([0, 0.6, 1.2, 1.8])
ax6.set_rlim(0, 1.8)
for ax in [ax1, ax2, ax3, ax4, ax5, ax6]:
ax.set_thetagrids(np.linspace(0, 360, 12, endpoint=False))
plt.tight_layout()
plt.show()
|
Decomposition of stress in Legendre polynomials¶
To compute the GGF, ggf.globgeomfact.coeff2ggf()
uses the
coefficients of the decomposition of the stress into Legendre
polynomials \(P_n(\text{cos}(\theta))\). This example visualizes
the small differences between the original stress and the stress
computed from the Legendre coefficients. This plot is automatically
produced by the original Matlab script StretcherNStress.m.
Note that the original Matlab yields different results for the
same set of parameters, because the Poisson’s ratio (keyword
argument poisson_ratio
) is non-zero;
see issue #1.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 | import matplotlib.pylab as plt
import numpy as np
import percache
from ggf.sci_funcs import legendrePlm
from ggf.stress.boyde2009.core import stress
@percache.Cache("stress_decomposition.cache", livesync=True)
def compute(**kwargs):
"Locally cached version of ggf.core.stress"
return stress(**kwargs)
# compute default stress
theta, sigmarr, coeff = compute(ret_legendre_decomp=True,
n_points=300)
# compute stress from coefficients
numpoints = theta.size
sigmarr_c = np.zeros((numpoints, 1), dtype=float)
for ii in range(numpoints):
for jj, cc in enumerate(coeff):
sigmarr_c[ii] += coeff[jj] * \
np.real_if_close(legendrePlm(0, jj, np.cos(theta[ii])))
# polar plot
plt.figure(figsize=(8, 8))
ax = plt.subplot(111, projection="polar")
plt.plot(theta, sigmarr, '-b', label="computed stress")
plt.plot(theta + np.pi, sigmarr[::-1], '-b')
plt.plot(theta, sigmarr_c, ':r', label="from Legendre coefficients")
plt.plot(theta + np.pi, sigmarr_c[::-1], ':r')
plt.legend()
plt.tight_layout()
plt.show()
|
Object boundary: stretching and Poisson’s ratio¶
This example illustrates how the parameters Poisson’s ratio
\(\nu\) and stretch ratio \(\epsilon\) influence
the object boundary used in ggf.core.stress()
and
defined in ggf.core.boundary()
.
Note that the boundary function was not defined correctly prior to version 0.3.0 (issue #1). Since version 0.3.0, the semi-minor axis is equivalent to the keyword argument a (1 by default).
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 | import numpy as np
import matplotlib.pylab as plt
from ggf.stress.boyde2009.core import boundary
theta = np.linspace(0, 2*np.pi, 300)
costheta = np.cos(theta)
# change epsilon
eps = [.0, .05, .10, .15, .20]
b1s = []
for ep in eps:
b1s.append(boundary(costheta=costheta,
epsilon=ep,
nu=.0))
# change Poisson's ratio
nus = [.0, .25, .5]
b2s = []
for nu in nus:
b2s.append(boundary(costheta=costheta,
epsilon=.1,
nu=nu))
# plot
plt.figure(figsize=(8, 4))
ax1 = plt.subplot(121, projection="polar")
for ep, bi in zip(eps, b1s):
ax1.plot(theta, bi, label="ϵ={:.2f}, ν=0".format(ep))
ax1.legend()
ax2 = plt.subplot(122, projection="polar")
for nu, bi in zip(nus, b2s):
ax2.plot(theta, bi, label="ϵ=.1, ν={:.1f}".format(nu))
ax2.legend()
plt.tight_layout()
plt.show()
|