<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">'''
 TrapAnalysis.py

 Functions for analyzing position data obtained from the video of
 a trapped bead. Produces information about the optical trap.

 Author: Christopher Dydula, University of Toronto    -    May 23, 2011
'''

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D

def my_show():
    '''Display a figure. Same as show() in backend_tkagg.py but without
    Tk.mainloop.
    '''
    for manager in plt._pylab_helpers.Gcf.get_all_fig_managers():
        manager.show()
    
def positions_plot(x, y, frames):
    '''Create a 3d plot of the x and y positions vs frame from the position data
       of a trapped bead video.

       x: array_like
          Array containing the x positions
       y: array_like
          Array containing the y positions
       frames: array_like
          Array containing the frame numbers
    '''

    fig = plt.figure(1)
    ax = Axes3D(fig)
    ax.scatter(x, y, frames, label='Positions', c=frames)
    #ax.legend(loc=9) # legend is placed at top centre, acts as title

    ax.set_xlabel('x (m)')
    ax.set_ylabel('y (m)')
    ax.set_zlabel('Frame')
    ax.set_title('Position vs. frame')

    return fig
   # my_show()

def disp_v_frame(direc, frames, varp, name):
    '''Create a plot of a vector of positions in the direction direc vs a
       vector of frames.

       direc: array_like
              Array of positions in the desired direction (x, y or radial)
       frames: array_like
              Array of frame numbers
       varp: array_like
              An array with n variances for the positions in direc, where
              n is the length of direc and frames
       name: string
              Either 'x', 'y' or 'r' for the three respective directions
    '''
    plt.plot(frames, direc, frames, varp, linewidth=0.5)
    plt.xlabel('Frame')
    if name == 'x':
        plt.ylabel('x (m)')
        plt.title('x displacement vs. frame')
        plt.legend(('x for individual frames', '$&lt;x^2&gt;$ = %9.4e' % varp[0]))
    elif name == 'y':
        plt.ylabel('y (m)')
        plt.title('y displacement vs. frame')
        plt.legend(('y for individual frames', '$&lt;y^2&gt;$ = %9.4e' % varp[0]))
    elif name == 'r':
        plt.ylabel('$r^2 (m^2)$')
        plt.title('Radial displacement vs. frame')
        plt.legend(('$r^2$ for individual frames', '$&lt;r^2&gt;$ = %9.4e' % varp[0]))
    plt.grid()

def disp_distr(direc, bins, name):
    '''Create a histogram showing the distribution of the displacements in the
       direction direc (x, y or radial)
       
       direc: array_like
              Array of positions in the desired direction (x, y or radial)
        bins: positive integer
              The number of bins for the histogram
        name: string
              Either 'x', 'y' or 'r' for the three respective directions
    '''
    plt.hist(direc, bins)
    plt.ylabel('Count')
    if name == 'x':
        plt.title('Distribution of x displacements')
        plt.xlabel('Displacement (m)')
    elif name == 'y':
        plt.title('Distribution of y displacements')
        plt.xlabel('Displacement (m)')
    elif name == 'r':
        plt.title('Distribution of radial displacements')
        plt.xlabel('Displacement ($m^2$)')

def rms(x):
    '''Return the root mean square of the values in the array vector x.

       x: array_like
          Array of a numbers
    '''
    return (np.sum(x**2) / np.size(x)) ** (0.5)

def analyze(x, y, temp, psize_um, deltemp, delpsize_um):
    '''Given x and y pixel position data of a bead from a video, convert
       pixel units to metric. Produce 7 plots in two figure windows, a 3d
       position plot in the first and 3 position plots and 3 displacement
       distrubution histograms in the second. Calculate and return trap
       stiffness in the x, y and radial direction.

       x: array_like
          Array containing the x positions (pixels)
       y: array_like
          Array containing the y positions (pixels)

       x and y must be of the same length
    '''
    kb = 1.38065e-23 # [m^2*kg*s^-2*K^-1], boltzmann constant
    T = temp # [K], temperature in room at time of recording of video
    delT = deltemp
    degfree = 2 # number of degrees of freedom

    # pixel sizes and uncertainties in metres
    # by default, every 14.5 pixels in the image is one micrometre, approximately
    # (1 micrometre is 14.5 +/- 0.3 px)
    psize_m = psize_um * 1e-6
    delpsize_m = delpsize_um * 1e-6

    # zero-mean postions in metres
    x -= np.mean(x)
    y -= np.mean(y)
    x_m = x * psize_m
    y_m = y * psize_m
    # x_m -= np.mean(x_m)
    # y_m -= np.mean(y_m)

    # use the standard deviations of x, y as proxies for the pixel uncertainties
    errxpx = np.std(x)
    errypx = np.std(y)

    # uncertainties in x, y values (metres)
    delx_m = ( (psize_m * errxpx)**2 + (x * delpsize_m)**2 )**(0.5)
    dely_m = ( (psize_m * errypx)**2 + (y * delpsize_m)**2 )**(0.5)

    # uncertainty in &lt;r^2&gt;
    delrvar_m = (2/(np.size(x))) * ( np.sum((delx_m*x_m)**2) + np.sum((dely_m*y_m)**2) )**(0.5)

    # calculate &lt;r^2&gt; and trap stiffness
    r2 = x_m**2 + y_m**2
    rvar = np.sum(r2) / (np.size(r2)) # - 1)
    delk = degfree * kb * np.sqrt( (delT/rvar)**2 + (T*delrvar_m/(rvar**2))**2 )
    ktrap = degfree * kb * T / rvar

    # create a frame vector
    n = np.size(x)
    frames = np.linspace(1,n,n)

    # create a 3d plot of the x and y positions vs frame
    mpl.rcParams['legend.fontsize'] = 8
    mpl.rcParams['figure.subplot.left'] = 0.07
    mpl.rcParams['figure.subplot.right'] = .95
    mpl.rcParams['figure.subplot.hspace'] = .2
    mpl.rcParams['ytick.labelsize'] = 8
    mpl.rcParams['xtick.labelsize'] = 8
    mpl.rcParams['font.size'] = 10
    fig1 = positions_plot(x_m, y_m, frames)

    # plot the x, y and radial displacements vs frame as well as histograms
    # of the x, y and radial distributions in the same figure
    fig2 = plt.figure(2, figsize=(12,8))

    xvar = np.var(x_m)
    yvar = np.var(y_m)
    rvar = rvar
    xvarp = xvar * np.ones((n,1)) # a vector with n xvar's
    yvarp = yvar * np.ones((n,1)) # a vector with n yvar's
    rvarp = rvar * np.ones((n,1)) # a vector with n rvar's

    plt.subplot(2,3,1)
    disp_v_frame(x_m, frames, xvarp, 'x')
    plt.subplot(2,3,2)
    disp_v_frame(y_m, frames, yvarp, 'y')
    plt.subplot(2,3,3)
    disp_v_frame(r2, frames, rvarp, 'r')

    bins = n // 10 # the number of bins for the histograms; may be changed

    plt.subplot(2,3,4)
    disp_distr(x_m, bins, 'x')
    plt.subplot(2,3,5)
    disp_distr(y_m, bins, 'y')
    plt.subplot(2,3,6)
    disp_distr(r2, bins, 'r')
    plt.tight_layout()

    figs = [fig1, fig2]

    return ktrap, delk, figs
</pre></body></html>