Code examples

Applications

Creep compliance analysis

This example uses the contour data of an 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).

creep_compliance.py

  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
117
118
119
120
121
122
123
124
125
import ggf
import h5py
import lmfit
import matplotlib.pylab as plt
import numpy as np
import percache

mycache = percache.Cache("creep_compliance.cache", livesync=True)


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

@mycache
def get_ggf(**kw):
    f = ggf.get_ggf(use_lut=True, **kw)
    return f

# 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
    print("compute ggf smaj={:.3e}, smin={:.3e}".format(smaj, smin))
    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]"]
    else:
        power_per_fiber=meta["power_per_fiber_trap [W]"]
    f = 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)
    print("... ", ii, f)
    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="compliance")
ax4.plot(time[compl_ival], complnc_1[compl_ival], label="single GGF")
ax4.plot(time[compl_ival], complnc[compl_ival], label="series GGF")
ax4.legend()
ax4.set_xlabel("time [s]")
ax4.set_ylabel("creep 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].

_images/stress_reproduced.jpg

stress_reproduced.py

  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,
          "radius": 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.

_images/stress_decomposition.jpg

stress_decomposition.py

 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,
                                numpoints=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().

_images/boundary.jpg

boundary.py

 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 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()