From b6dbee3b4b05496f89ecbd3282164f80daefb612 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Wed, 2 Dec 2009 15:22:15 -0500 Subject: [PATCH] Added latex framework and cantilever_calib appendix. --- tex/Makefile | 14 ++ tex/src/Makefile | 11 ++ tex/src/cantilever_calib/Makefile | 17 ++ tex/src/cantilever_calib/README | 8 + tex/src/cantilever_calib/cantilever_calib.tex | 58 ++++++ .../cantilever_calib/contour_integration.tex | 61 ++++++ tex/src/cantilever_calib/dot/concept_map.dot | 60 ++++++ tex/src/cantilever_calib/integrals.tex | 145 +++++++++++++++ tex/src/cantilever_calib/mp/Makefile | 31 +++ tex/src/cantilever_calib/mp/contour.mp | 55 ++++++ tex/src/cantilever_calib/mp/sample.tex | 12 ++ tex/src/cantilever_calib/overview.tex | 176 ++++++++++++++++++ tex/src/cantilever_calib/setup_general.tex | 64 +++++++ tex/src/cantilever_calib/solve_general.tex | 60 ++++++ .../cantilever_calib/solve_highly_damped.tex | 64 +++++++ tex/src/cantilever_calib/test/Lorentz_root.py | 44 +++++ tex/src/cantilever_calib/test/Makefile | 9 + tex/src/cantilever_calib/test/locate_roots.py | 51 +++++ .../cantilever_calib/test/test_integral.py | 80 ++++++++ tex/src/cantilever_calib/test/test_roots.py | 65 +++++++ tex/src/local_cmmds.tex | 59 ++++++ tex/src/packages.tex | 25 +++ tex/src/root.tex | 54 ++++++ tex/src/wtk.bib | 1 + 24 files changed, 1224 insertions(+) create mode 100644 tex/Makefile create mode 100644 tex/src/Makefile create mode 100644 tex/src/cantilever_calib/Makefile create mode 100644 tex/src/cantilever_calib/README create mode 100644 tex/src/cantilever_calib/cantilever_calib.tex create mode 100644 tex/src/cantilever_calib/contour_integration.tex create mode 100644 tex/src/cantilever_calib/dot/concept_map.dot create mode 100644 tex/src/cantilever_calib/integrals.tex create mode 100644 tex/src/cantilever_calib/mp/Makefile create mode 100644 tex/src/cantilever_calib/mp/contour.mp create mode 100644 tex/src/cantilever_calib/mp/sample.tex create mode 100644 tex/src/cantilever_calib/overview.tex create mode 100644 tex/src/cantilever_calib/setup_general.tex create mode 100644 tex/src/cantilever_calib/solve_general.tex create mode 100644 tex/src/cantilever_calib/solve_highly_damped.tex create mode 100644 tex/src/cantilever_calib/test/Lorentz_root.py create mode 100644 tex/src/cantilever_calib/test/Makefile create mode 100644 tex/src/cantilever_calib/test/locate_roots.py create mode 100644 tex/src/cantilever_calib/test/test_integral.py create mode 100644 tex/src/cantilever_calib/test/test_roots.py create mode 100644 tex/src/local_cmmds.tex create mode 100644 tex/src/packages.tex create mode 100644 tex/src/root.tex create mode 120000 tex/src/wtk.bib diff --git a/tex/Makefile b/tex/Makefile new file mode 100644 index 0000000..2dbb26c --- /dev/null +++ b/tex/Makefile @@ -0,0 +1,14 @@ +thesis.pdf : build + (cd ./build && xelatex root) + (cd ./build && bibtex root) + (cd ./build && xelatex root) + (cd ./build && xelatex root) + mv build/root.pdf $@ + +.PHONY: build +build : + rsync -a src/ build + (cd ./build && make) + +clean : + rm -rf build thesis.pdf diff --git a/tex/src/Makefile b/tex/src/Makefile new file mode 100644 index 0000000..8e0b458 --- /dev/null +++ b/tex/src/Makefile @@ -0,0 +1,11 @@ +SUBDIRS = cantilever_calib + +all : + @for i in $(SUBDIRS); do \ + echo "make all in $$i..."; \ + (cd $$i; $(MAKE) $(MFLAGS) all); done + +clean : + @for i in $(SUBDIRS); do \ + echo "make clean in $$i..."; \ + (cd $$i; $(MAKE) $(MFLAGS) clean); done diff --git a/tex/src/cantilever_calib/Makefile b/tex/src/cantilever_calib/Makefile new file mode 100644 index 0000000..f1b083b --- /dev/null +++ b/tex/src/cantilever_calib/Makefile @@ -0,0 +1,17 @@ +GEN_FILES = mp/contour.1 dot/concept_map.png + +all : $(GEN_FILES) + +clean : + rm -f $(GEN_FILES) + $(MAKE) -C mp clean + $(MAKE) -C test clean + +check : + $(MAKE) -C test all + +mp/contour.1 : + $(MAKE) -C mp images + +dot/concept_map.png : + dot -Tpng dot/concept_map.dot > $@ diff --git a/tex/src/cantilever_calib/README b/tex/src/cantilever_calib/README new file mode 100644 index 0000000..368f355 --- /dev/null +++ b/tex/src/cantilever_calib/README @@ -0,0 +1,8 @@ +I didn't have a good understanding of the theory behind thermally +calibrating an AFM cantilever, so I work it out here with all the +gorey details :p. + +The testq subdirectory contains some python scripts I used to test my +algebra and get a better feel for what was going on. The dot +subdirectory is a preliminary data-flow map of the cantilever +calibration procedure. diff --git a/tex/src/cantilever_calib/cantilever_calib.tex b/tex/src/cantilever_calib/cantilever_calib.tex new file mode 100644 index 0000000..ea6b2fb --- /dev/null +++ b/tex/src/cantilever_calib/cantilever_calib.tex @@ -0,0 +1,58 @@ +\input{cantilever_calib/overview} +\input{cantilever_calib/setup_general} +\input{cantilever_calib/solve_highly_damped} +\input{cantilever_calib/solve_general} +\input{cantilever_calib/contour_integration} +\input{cantilever_calib/integrals} + +%\begin{thebibliography}{77} %77 is width of longest number of your reference? +%% Defines the standard Lorentzian function +%\bibitem{mathworld_lorentzian} +% Eric W.\ Weisstein. +% {\it Lorentzian Function.} +% MathWorld--A Wolfram Web Resource. +% \url{http://mathworld.wolfram.com/LorentzianFunction.html} +% +%\bibitem{cos_halfangle} +% S.\ Thornton and J.\ Marion. +% {Classical Dynamics of Particles and Systems}. +% 5\sth Edition, Brooks/Cole, 2004. \\ +% See Eq. D.16 on page 609. +% +%\bibitem{four_deriv} +% Hmm, it is suprisingly difficult to find an `official' reference for this. +% I obviously need to get a spectral analysis book :p. +% See Wikipedia's currently excellent page (Feb 15th, 2008) \\ +% \url{http://en.wikipedia.org/wiki/Fourier_transform#Functional_relationships} ,\\ +% or derive it for yourself in about three lines :p. +% +%\bibitem{parseval} +% W.\ Press, S.\ Teukolsky, W.\ Vetterling, and B.\ Flannery. +% {\it Numerical Recipies in C: The Art of Scientific Computing}. +% 2\nd Edition, Cambridge University Press, 1992. \\ +% \url{http://www.nrbook.com/a/bookcpdf.php}\\ +% See Eq. 12.0.13 on page 498. +% +%\bibitem{PSD} +% Press et al. See Eq. 12.0.14 on page 498. +% +% +%\bibitem{wiener_khinchin} +% Press et al. See Eq. 12.0.12 on page 498. +% +%\bibitem{wikipedia_wiener_khinchin} +% Wikipedia. +% {\it Wiener-Khinchin theorem}. \\ +% \url{http://en.wikipedia.org/wiki/Wiener\%E2\%80\%93Khinchin_theorem} +% +%\bibitem{tweezer_lab_notes} +% C.\ Grossman and A.\ Stout. +% {\it Optical Tweezers Advanced Lab} +% Fall 2005. See section 4.3. +% \url{http://chirality.swarthmore.edu/PHYS81/OpticalTweezers.pdf} +% +% \url{http://nlolab.swarthmore.edu/webstuff/Phys81_82/OpticalTweezersLabTheory.pdf} +% +%\end{thebibliography} +% +%\end{document} diff --git a/tex/src/cantilever_calib/contour_integration.tex b/tex/src/cantilever_calib/contour_integration.tex new file mode 100644 index 0000000..c8b5de9 --- /dev/null +++ b/tex/src/cantilever_calib/contour_integration.tex @@ -0,0 +1,61 @@ +\section{Contour integration} + +\begin{figure} + \begin{centering} + %\includegraphics[width=0.3\textwidth]{mp/contour.1} \\ + \caption{Integral contour \C\ enclosing the upper half plane.} + \label{f.UHP_contour} + \end{centering} +\end{figure} + +As a brief review, some definite integrals from $-\infty$ to $\infty$ +can be evaluated by integrating along the contour \C\ +shown in Figure \ref{f.UHP_contour}. + +A sufficient condition on the function $f(z)$ to be integrated, is +that $\lim_{|z|\rightarrow\infty}|f(z)|$ falls off at least as fast as +$\frac{1}{z^2}$. +When this is the case, the integral around the outer semicircle of \C\ is 0, +so the $\iC{f(z)} = \iInfInf{z}{f(z)}$. + +We can evaluate the integral using residue theorem +\begin{equation} + \iC{f(x)} = \sum_{z_p \in \text{poles in \C}} 2\pi i \Res{z_p}{f(z)}, + \label{eq.res_thm} +\end{equation} +where for simple poles (single roots) +\begin{equation} + \Res{z_p}{f(z)} = \limZp(z-z_p) f(z), \label{eq.res_simple} +\end{equation} +and in general for a pole of order $n$ +\begin{equation} + \Res{z_p}{f(z)} = \frac{1}{(n-1)!} \cdot\limZp + \nderiv{n-1}{z}{}\left[ (z-z_p)^n \cdot f(z) \right] + \label{eq.res_general} +\end{equation} + +%We also use the following theorem +%\begin{thm} \label{thm.res_even_fn} +%If $f(z) = f(-z)$, then $\Res{z_p}{f(z)} = \Res{-z_p}{f(z)}$.\\ +%\end{thm} +%\begin{proof} +%If $g(z)$ is even, then $\deriv{z}{g}$ must be odd +% (consider the Lorentz expansion). +%By induction, $\nderiv{n}{z}{g}$ is even when $n$ is even, +%and odd when $n$ is odd. +%By inspection, $(z-z_p)^n$ is also even iff $n$ is even. +%Therefore, +%\begin{equation} +% (z-z_p)^n \cdot f(z) +%\end{equation} +%is even iff $n$ is even, so +%\begin{equation} +% \nderiv{n-1}{z}{}\left( (z-z_p)^n \cdot f(z) \right) +%\end{equation} +%is always even, so +%\begin{equation} +% \Res{z_p}{f(z)} = \frac{1}{(n-1)!} \cdot\limZp +% \nderiv{n-1}{z}{}\left[ (z-z_p)^n \cdot f(z) \right] +%\end{equation} +%must be even with respect to $z_p$ $\qed$. +%\end{proof} diff --git a/tex/src/cantilever_calib/dot/concept_map.dot b/tex/src/cantilever_calib/dot/concept_map.dot new file mode 100644 index 0000000..7bf5d60 --- /dev/null +++ b/tex/src/cantilever_calib/dot/concept_map.dot @@ -0,0 +1,60 @@ +graph message_handler { + + epsilon=0.01; + nodesep=0.1; /* inches */ + + edge [ +// arrowhead=normal, +// dir=both, +// len=1 + ]; + + /* concept nodes */ + node [ + shape=box, +// fontsize=12, + color=black, + fillcolor=white, + style="rounded,filled" + ]; + + k_spring [label="Spring constant"]; + equipartition [label="Equipartition theorem"]; + + temperature [label="Temperature"]; + + photo_sense [label="Photodiode sensitivity"]; + + voltage_var_spring [label="Voltage variance\ndue to the spring"]; + psd_fit [label="Fit power spectral density"]; + psd_measured [label="Measured power spectral density"]; + + + /* data nodes */ + node [ + fillcolor="#ffd83d", /* drexel yellow */ + ]; + + thermocouple [label="Thermocouple voltage"]; + contact_slope [label="Surface contact"]; + voltage_noise [label="Voltage noise"]; + + /* Deeper calibration nodes */ + node [ + fillcolor="#ff9999", /* reddish */ + ]; + zp_calib [label="Z-piezo calibration"]; + volt_calib [label="Voltmeter calibration"]; + thermocouple_calib [label="Thermocouple calibration"]; + + k_spring -- equipartition; + + equipartition -- temperature -- thermocouple; + equipartition -- photo_sense -- contact_slope; + equipartition -- voltage_var_spring -- psd_fit -- psd_measured -- voltage_noise; + + temperature -- thermocouple_calib; + photo_sense -- zp_calib; + psd_measured -- volt_calib; + +} diff --git a/tex/src/cantilever_calib/integrals.tex b/tex/src/cantilever_calib/integrals.tex new file mode 100644 index 0000000..d9b6752 --- /dev/null +++ b/tex/src/cantilever_calib/integrals.tex @@ -0,0 +1,145 @@ +\section{Integrals} +\label{sec.integrals} + +\subsection{Highly damped integral} + +\begin{align} + I &= \iOInf{z}{\frac{1}{k^2 + z^2}} \\ + &= \frac{1}{2} \iInfInf{z}{\frac{1}{k^2 + z^2}} \\ + &= \frac{1}{2k} \iInfInf{u}{\frac{1}{u^2+1}} \\ +\end{align} +where $u \equiv z/k$, $du = dz/k$. +There are simple poles at $u = \pm i$ +\begin{align} + I &= \frac{1}{2k} \cdot 2 \pi i \Res{i}{f(u)} \\ + &= \frac{1}{2k} \cdot \frac{2 \pi i}{i+i} \\ + &= \frac{1}{2k} \pi \\ + &= \frac{\pi}{2 k}, +\end{align} + +\subsection{General case integral} + + +We will show that for any $(a,b > 0) \in \Reals$ +\begin{equation} + I = \iInfInf{z}{\frac{1}{(a^2-z^2) + b^2 z^2}} = \frac{\pi}{b a^2}. +\end{equation} + +First we note that $|f(z)| \rightarrow 0$ like $|z^{-4}|$ for $|z| \gg 1$, +and that $f(z)$ is even, so +\begin{equation} + I = \iC{\frac{1}{(a^2-z^2)^2 + b^2 z^2}}, +\end{equation} +where \C\ is the contour shown in Figure \ref{f.UHP_contour}. + +Because the denominator is of the form $A^2 + B^2$, we can factor it +into $(A+iB)(A-iB)$ like so % thanks Prof. Yuan +\begin{equation} + (a^2-z^2)^2 + b^2 z^2 + = (a^2-z^2 \colA{+} ibz)(a^2-z^2 \colA{-} ibz) +\end{equation} +And the roots of $z^2 \colA{\pm} ibz - a^2$ +\begin{align} + z_{r\colB{\pm}} + &= \colA{\pm}\frac{ib}{2} \left( + 1 \colB{\pm} \sqrt{1-4\frac{-a^2}{(ib)^2}} + \right) \\ + &= \pm\frac{ib}{2} \left( + 1 \pm \sqrt{1-4\frac{a^2}{b^2}} + \right) \\ + &= \pm\frac{ib}{2} \left( + 1 \pm S + \right) +\end{align} +Where $S \equiv \sqrt{1-4\frac{a^2}{b^2}}$. + +%critical damping when $\omega_0^2 = \beta'^2$ % TM +%where our $a = \omega_0$ and $b = \beta$, +%and $\beta = \gamma/m = 2 \beta'$ +%Critical damping when $a^2 = b^2/4$, so $S = 0$ +To determine the nature and locations of the roots, consider the following +cases (in order of increasing $a$). +\begin{itemize} + \item $a < b/2$, overdamped. + \item $a = b/2$, critically damped. + \item $a > b/2$, underdamped. +\end{itemize} + +In the overdamped case $S \in \Reals$ and $S > 0$, +so $z_{r\pm}$ is purely imaginary, and $z_{r+} != z_{r-}$. +For any $a < b/2$, we have $0 < S < 1$, so $\Imag(z_{r\pm}) > 0$. +Thus, there are two single poles in the upper half plane ($z_{r\pm}$), +and two single poles in the lower half plane ($-z_{r\pm}$). + +In the critically damped case $S = 0$, so $z_{r+} = z_{r-}$, +and we have double poles at $\pm z_{r+} = \frac{ib}{2}$. + +In the underdamped case $S$ is purely imaginary, +so $z_{r\pm}$ is complex, with $z_{r+}$ in the 2\nd quarter, +and $z_{r-}$ in the 1\st quarter. +The other two simple poles, $-z_{r-}$ and $-z_{r+}$, are in the +3\rd and 4\sth quarters respectively. + +Our contour \C\ always encloses the poles $z_{r\pm}$. +We will deal with the simple pole cases first, +and then return to the critically damped case. + +\subsubsection{Over- and under-damped} + +Our factored function $f(z)$ is +\begin{equation} + f(z) = \frac{1}{(z-z_{r+})(z+z_{r+})(z+z_{r-})(z-z_{r-})} +\end{equation} + +Applying eq. \ref{eq.res_thm} and \ref{eq.res_simple} we have +\begin{align} + I &= 2\pi i \left( \Res{z_{r+}}{f(z)} + \Res{z_{r-}}{f(z)} \right) \\ + &= 2\pi i \left( + \frac{1}{ (z_{r+}+z_{r+}) + \colA{(z_{r+}+z_{r-}) + (z_{r+}-z_{r-})} } + + \frac{1}{ \colA{(z_{r-}-z_{r+}) + (z_{r-}+z_{r+})} + (z_{r-}+z_{r-}) } + \right) \\ + &= \frac{\pi i}{\colA{z_{r+}^2-z_{r-}^2}} \left( + \frac{1}{z_{r+}} + \colA{-} \frac{1}{z_{r-}} + \right) \\ + &= \frac{\pi i}{ \left( \colB{\frac{ib}{2}} (1+S) \right)^2 + - \left( \colB{\frac{ib}{2}} (1-S) \right)^2 } + \cdot \frac{z_{r-}-z_{r+}}{z_{r+}z_{r-}} \\ + &= \frac{\colB{-4}\pi i / \colB{b^2}}{ (1+2S+S^2) - (1-2S+S^2) } + \cdot \frac{ \colA{\frac{ib}{2}} [(1-S) - (1+S)] } + { \left(\frac{ib}{2}\right)^{\colA{2}} (1+S)(1-S) } \\ + &= \frac{-8\pi / b^3}{ 4S } + \cdot \frac{-2S} + {(1 - S^2)} \\ + &= \frac{ 4\pi }{ b^3 (1 - S^2)} \\ + &= \frac{ 4\pi }{ b^3 [1 - (1-4\frac{a^2}{b^2})]} \\ + &= \frac{ 4\pi }{ b^3 \cdot 4\frac{a^2}{b^2}} \\ + &= \frac{ \pi }{ b a^2 } \label{eq.gen_int_noncrit} +\end{align} +Hooray! + +\subsubsection{Critically damped} + +Our factored function $f(z)$ is +\begin{equation} + f(z) = \frac{1}{(z-z_{r+})^2(z-z_{r-})^2} +\end{equation} + +Applying eq. \ref{eq.res_thm} and \ref{eq.res_general} we have +\begin{align} + I &= 2\pi i \Res{z_{r+}}{f(z)} \\ + &= \colA{2}\pi i \left( \colA{\frac{1}{2!}} + \limZ{z_{r+}} + \deriv{z}{} \frac{1}{(z + z_{r+})^2} + \right) \\ + &= \pi i \limZ{z_{r+}} -2 \cdot \frac{1}{(z_{r+} + z_{r+})^3} \\ + &= - 2 \pi i \frac{1}{z_{r+}^3} \\ + &= \colA{-} 2 \pi \colA{i} \frac{1}{(\frac{\colA{i}b}{2})^3} \\ + &= \frac{\pi}{b (\frac{b}{2})^2} \\ + &= \frac{\pi}{b a^2}, \label{eq.gen_int_crit} +\end{align} +which matches eq. \ref{eq.gen_int_noncrit} diff --git a/tex/src/cantilever_calib/mp/Makefile b/tex/src/cantilever_calib/mp/Makefile new file mode 100644 index 0000000..be5e052 --- /dev/null +++ b/tex/src/cantilever_calib/mp/Makefile @@ -0,0 +1,31 @@ +MPS = contour + + +TEMP_FILES = $(MPS:%=%.mp[a-z]) $(MPS:%=%.log) mp* texnum.mpx +GENERATED_FILES = $(MPS:%=%.[0-9]*) sample.log sample.aux sample.pdf + +#PDFVIEWER = gv +PDFVIEWER = xpdf +#PDFVIEWER = evince + +images : $(MPS:%=%.1) + +all : view + +view : sample.pdf + $(PDFVIEWER) sample.pdf & + +clean : semi-clean + rm -f $(GENERATED_FILES) + +semi-clean : + rm -f $(TEMP_FILES) + +# generate a pdf containing all mp images for previewing-troubleshooting +# depend on the first image from each mp file +sample.pdf : sample.tex images + pdflatex sample.tex + +# if we call for the first image from an mp file, make them all +%.1 : %.mp + mpost $< diff --git a/tex/src/cantilever_calib/mp/contour.mp b/tex/src/cantilever_calib/mp/contour.mp new file mode 100644 index 0000000..ce60879 --- /dev/null +++ b/tex/src/cantilever_calib/mp/contour.mp @@ -0,0 +1,55 @@ +%% Sarcomere graphic + +boolean labels; labels := false; + +% coloring in rgb +color axis_color, contour_color, text_color; +%acolor = (0,.1736,.62); % drexel blue +%mcolor = (1,0.85,0.24); % drexel yellow +axis_color = (0,0,0); % black +contour_color = (0,0,1); % blue +text_color = (0,0,0); % black + +% sizing +numeric radius; +radius := 0.7; +u := 2cm; % define units + +% Generate and draw axes with labels. +% Origin at CENTER, +% LL (lower left) and UR (upper right) define the length of the axis markings +% XLABEL and YLABEL contain text to label the axes +def draw_axis(expr center, ll, ur, xlabel, ylabel) = + save start, stop; + pair start, stop; + + % draw X + start := (ll dotprod (1,0),0); % get the x portion + stop := (ur dotprod (1,0),0); + drawarrow start .. stop withcolor axis_color; + label.bot(xlabel, stop) withcolor text_color; + + % draw Y + start := (0, ll dotprod (0,1)); % get the y portion + stop := (0, ur dotprod (0,1)); + drawarrow start .. stop withcolor axis_color; + label.lft(ylabel, stop) withcolor text_color; +enddef; + +% Generate and draw a semicircle in the upper half plane +% Centered at CENTER with radius RAD. +def draw_UHP_semicircle(expr center, rad) = + % straight line + draw (center-(rad,0)) .. (center+(rad,0)) withcolor contour_color; + % semicircle + draw (center-(rad,0)){up} .. tension(1) .. {down}(center+(rad,0)) + withcolor contour_color; + drawarrow center .. (center+(rad/2,0)) withcolor contour_color; +enddef; + +beginfig(1) + draw_axis(origin, -(u,0.2u), (u,u), "Re", "Im"); + draw_UHP_semicircle(origin, radius * u); +endfig; + +end diff --git a/tex/src/cantilever_calib/mp/sample.tex b/tex/src/cantilever_calib/mp/sample.tex new file mode 100644 index 0000000..e8fca4f --- /dev/null +++ b/tex/src/cantilever_calib/mp/sample.tex @@ -0,0 +1,12 @@ +\documentclass{article} +\usepackage{graphicx} +\usepackage{fullpage} + +\DeclareGraphicsRule{*}{mps}{*}{} % if you don't recognize the type, it's mps + +\newcommand{\cmpi}[1]{contour #1 \includegraphics{contour.#1}\vspace{1cm} \\} + +\begin{document} +\centering +\cmpi{1} +\end{document} diff --git a/tex/src/cantilever_calib/overview.tex b/tex/src/cantilever_calib/overview.tex new file mode 100644 index 0000000..2948376 --- /dev/null +++ b/tex/src/cantilever_calib/overview.tex @@ -0,0 +1,176 @@ +\section{Overview} + +In order to measure forces accurately with an Atomic Force Microscope (AFM), +it is important to measure the cantilever spring constant. +The force exerted on the cantilever can then be deduced from it's deflection +via Hooke's law $F = -kx$. + +The basic idea is to use the equipartition theorem\cite{hutter93}, +\begin{equation} + \frac{1}{2} k \avg{x^2} = \frac{1}{2} k_BT \label{eq.equipart}, +\end{equation} +where $k_B$ is Boltzmann's constant, + $T$ is the absolute temperature, and + $\avg{x^2}$ denotes the expectation value of $x^2$ as measured over a + very long interval $t_T$, +\begin{equation} + \avg{A} \equiv \iLimT{A}. +\end{equation} +Solving the equipartition theorem for $k$ yields +\begin{equation} + k = \frac{k_BT}{\avg{x^2}}, \label{eq.equipart_k} +\end{equation} +so we need to measure (or estimate) the temperature $T$ and variance +of the cantilever position $\avg{x^2}$ in order to estimate $k$. + +\subsection{Related papers} + +Various corrections taking into acount higher order modes +\cite{butt95,stark01}, and cantilever tilt \cite{hutter05} have been +proposed and reviewed \cite{florin95,levy02,ohler07}, but we will +focus here on the derivation of Lorentzian noise in damped simple +harmonic oscillators that underlies all frequency-space methods for +improving the basic $k\avg{x^2} = k_BT$ method. + +Roters and Johannsmann describe a similar approach to deriving the Lorentizian +power spectral density\cite{roters96}. %, +%as do +% see Gittes 1998 for more thermal noise details +% see Berg-Sorenson for excellent overdamped treament. + +\emph{WARNING}: It is popular to refer to the power spectral density +as a ``Lorentzian''\cite{hutter93,roters96,levy02,florin95} even +though eq.~\ref{eq.model_psd} differs from the classic +Lorentzian\cite{mathworld_lorentzian}. +\begin{equation} + L(x) = \frac{1}{\pi}\frac{\frac{1}{2}\Gamma} + {(x-x_0)^2 + \p({\frac{1}{2}\Gamma})^2} +\end{equation} +It is unclear whether the references are due to uncertainty about the +definition of the Lorentzian or to the fact that +eq.~\ref{eq.model_psd} is also peaked. In order to avoid any +uncertainty, we will leave eq.~\ref{eq.model_psd} unnamed. + +\section{Methods} + +To find $\avg{x^2}$, the raw photodiode voltages $V_p(t)$ are +converted to distances $x(t)$ using the photodiode sensitivity +$\sigma_p$ (the slope of the voltage vs.~distance curve of data taken +while the tip is in contact with the surface) via +\begin{equation} + x(t) = \frac{V_p(t)}{\sigma_p} +\end{equation} +Rather than computing the variance of $x(t)$ directly, we attempt to +filter out noise by fitting the spectral power density (\PSD) of +$x(t)$ to the theoretically predicted \PSD\ for a damped harmonic +oscillator (eq.~\ref{eq.model_psd}) +\begin{align} + \ddt{x} + \beta\dt{x} + \omega_0^2 x &= \frac{F_\text{thermal}}{m} \\ + \PSD(x, \omega) &= \frac{G_1}{(\omega_0^2-\omega^2)^2 + \beta^2\omega^2}, +\end{align} +where $G_1\equiv G_0/m^2$, $\omega_0$, and $\beta$ are used as the +fitting parameters (see eqn.s \ref{eq.model_psd}). The variance of +$x(t)$ is then given by eq.~\ref{eq.DHO_var} +\begin{equation} + \avg{x(t)^2} = \frac{\pi G_1}{2\beta\omega_0^2}, +\end{equation} +which we can plug into the equipartition theorem +(eqn.~\ref{eq.equipart}) yielding +\begin{align} + k = \frac{2 \beta \omega_0^2 k_BT}{\pi G_1}. +\end{align} + +From eqn. \ref{eq.GO}, we find the expected value of $G_1$ to be +\begin{equation} + G_1 \equiv G_0/m^2 = \frac{2}{\pi m} k_BT \beta. \label{eq.Gone} +\end{equation} + + +\subsection{Fitting deflection voltage directly} + +In order to keep our errors in measuring $\sigma_p$ seperate from +other errors in measuring $\avg{x(t)^2}$, we can fit the voltage +spectrum before converting to distance. +\begin{align} + \ddt{V_p}/\sigma_p + \beta\dt{V_p}/\sigma_p + \omega_0^2 V_p/\sigma_p + &= F_\text{thermal} \\ + \ddt{V_p} + \beta\dt{V_p} + \omega_0^2 V_p + &= \sigma_p\frac{F_\text{thermal}}{m} \\ + \ddt{V_p} + \beta\dt{V_p} + \omega_0^2 V_p + &= \frac{F_\text{thermal}}{m_p} \\ + \PSD(V_p, \omega) &= \frac{G_{1p}} + { (\omega_0^2-\omega^2)^2 + \beta^2\omega^2 } \\ + \avg{V_p(t)^2} &= \frac{\pi G_{1p}}{2\beta\omega_0^2} + = \frac{\pi \sigma_p^2 G_{1}}{2\beta\omega_0^2} + = \sigma_p^2 \avg{x(t)^2}, +\end{align} +where $m_p\equiv m/\sigma_p$, $G_{1p}\equiv G_0/m_p^2=\sigma_p^2 G_1$. +Plugging into the equipartition theorem yeilds +\begin{align} + k &= \frac{\sigma_p^2 k_BT}{\avg{V_p(t)^2}} + = \frac{2 \beta\omega_0^2 \sigma_p^2 k_BT}{\pi G_{1p}}. +\end{align} + +From eqn. \ref{eq.Gone}, we find the expected value of $G_{1p}$ to be +\begin{equation} + G_{1p} \equiv \sigma_p^2 G_1 = \frac{2}{\pi m} \sigma_p^2 k_BT \beta. + \label{eq.Gone-p} +\end{equation} + + +\subsection{Fitting deflection voltage in frequency space} + +Note: the math in this section depends on some definitions from +section \ref{sec.setup}. + +As yet another alternative, you could fit in frequency +$f\equiv\omega/2\pi$ instead of angular frequency $\omega$. But we +must be careful with normalization. Comparing the angular frequency +and normal frequency unitary Fourier transforms +\begin{align} + \Four{x(t)}(\omega) + &\equiv \frac{1}{\sqrt{2\pi}} \iInfInf{t}{x(t) e^{-i \omega t}} \\ + \Fourf{x(t)}(f) &\equiv \iInfInf{t}{x(t) e^{-2\pi i f t}} + = \iInfInf{t}{x(t) e^{-i \omega t}} + = \sqrt{2\pi}\cdot\Four{x(t)}(\omega=2\pi f), +\end{align} +from which we can translate the \PSD +\begin{align} + \PSD(x, \omega) &\equiv \normLimT 2 \magSq{ \Four{x(t)}(\omega) } \\ + \PSD_f(x, f) &\equiv \normLimT 2 \magSq{ \Fourf{x(t)}(f) } + = 2\pi \cdot \normLimT 2 \magSq{ \Four{x(t)}(\omega=2\pi f) } + = 2\pi \PSD(x, \omega=2\pi f). +\end{align} +The variance of the function $x(t)$ is then given by plugging into +eqn.~\ref{eq.parseval_var} (our corollary to Parseval's theorem) +\begin{align} + \avg{x(t)^2} &= \iOInf{\omega}{\PSD(x,\omega)} + = \iOInf{f}{\frac{1}{2\pi}\PSD_f(x,f)2\pi\cdot} + = \iOInf{f}{\PSD_f(x,f)}. +\end{align} +Therefore +\begin{align} + \PSD_f(V_p, f) &= 2\pi\PSD(V_p,\omega) + = \frac{2\pi G_{1p}}{(4\pi f_0^2-4\pi^2f^2)^2 + \beta^2 4\pi^2f^2} + = \frac{2\pi G_{1p}}{16\pi^4(f_0^2-f^2)^2 + \beta^2 4\pi^2f^2} + = \frac{G_{1p}/8\pi^3}{(f_0^2-f^2)^2 + \frac{\beta^2 f^2}{4\pi^2}} \\ + &= \frac{G_{1f}}{(f_0^2-f^2)^2 + \beta_f^2 f^2} \\ + \avg{V_p(t)^2} &= \frac{\pi G_{1f}}{2\beta_f f_0^2}. +% = \frac{\pi G_{1p} / (2\pi)^3}{2\beta/(2\pi) \omega_0^2/(2\pi)^2} +% = \frac{\pi G_{1p}}{2\beta\omega_0^2} = \avg{V_p(t)^2} % check! +\end{align} +where $f_0\equiv\omega_0/2\pi$, $\beta_f\equiv\beta/2\pi$, and +$G_{1f}\equiv G_{1p}/8\pi^3$. Finally +\begin{align} + k &= \frac{\sigma_p^2 k_BT}{\avg{V_p(t)^2}} + = \frac{2 \beta_f f_0^2 \sigma_p^2 k_BT}{\pi G_{1f}}. +\end{align} + +From eqn. \ref{eq.Gone}, we expect $G_{1f}$ to be +\begin{equation} + G_{1f} = \frac{G_{1p}}{8\pi^3} + = \frac{\sigma_p^2 G_1}{8\pi^3} + = \frac{\frac{2}{\pi m} \sigma_p^2 k_BT \beta}{8\pi^3} + = \frac{\sigma_p^2 k_BT \beta}{4\pi^4 m}. + \label{eq.Gone-f} +\end{equation} diff --git a/tex/src/cantilever_calib/setup_general.tex b/tex/src/cantilever_calib/setup_general.tex new file mode 100644 index 0000000..d391a05 --- /dev/null +++ b/tex/src/cantilever_calib/setup_general.tex @@ -0,0 +1,64 @@ +\section{Theoretical power spectral density for a damped harmonic oscillator} +\label{sec.setup} + +Our cantilever can be approximated as a damped harmonic oscillator +\begin{equation} + m\ddt{x} + \gamma \dt{x} + k x = F(t), \label{eq.DHO} + % DHO for Damped Harmonic Oscillator +\end{equation} +where $x$ is the displacement from equilibrium, + $m$ is the effective mass, + $\gamma$ is the effective drag coefficient, + $k$ is the spring constant, and + $F(t)$ is the external driving force. +During the non-contact phase of calibration, + $F(t)$ comes from random thermal noise. + +In the following analysis, we use the unitary, angular frequency Fourier transform normalization +\begin{equation} + \Four{x(t)} \equiv \frac{1}{\sqrt{2\pi}} \iInfInf{t}{x(t) e^{-i \omega t}} +\end{equation} + +We also use the following theorems (proved elsewhere): +\begin{align} + \cos\left(\frac{\theta}{2}\right) &= \pm\sqrt{\frac{1}{2}[1+\cos(\theta)]} + &\text{\cite{cos_halfangle},} \label{eq.cos_halfangle} \\ + \Four{\nderiv{n}{t}{x(t)}} &= (i \omega)^n x(\omega) + &\text{\cite{four_deriv},} \label{eq.four_deriv} \\ +% \Four{x*y} &= x(\omega) y(\omega), \label{eq.four_conv} +% & \text{and} \\ + \iInfInf{t}{\magSq{x(t)}} &= \iInfInf{\omega}{\magSq{x(w)}} + &\text{(Parseval's)\cite{parseval}.} \label{eq.parseval} +\end{align} +%where $x*y$ denotes the convolution of $x$ and $y$, +%\begin{equation} +% x*y \equiv \iInfInf{\tau}{x(t-\tau)y(\tau)}. +%\end{equation} +As a corollary to Parseval's theorem, we note that the one sided power spectral density per unit time (\PSD) defined by +\begin{align} + \PSD(x, \omega) &\equiv \normLimT 2 \left| x(\omega) \right|^2 + &\text{\cite{PSD}} \label{eq.psd_def} +\end{align} +relates to the variance by +\begin{align} + \avg{x(t)^2} + &= \iLimT{\magSq{x(t)}} + = \normLimT \iInfInf{\omega}{\magSq{x(\omega)}} + = \iOInf{\omega}{\PSD(x,\omega)}, \label{eq.parseval_var} +\end{align} +where $t_T$ is the total time over which data has been aquired. + + +We also use the Wiener-Khinchin theorem, +which relates the two sided power spectral density $S_{xx}(\omega)$ +to the autocorrelation function $r_{xx}(t)$ via +\begin{align} + S_{xx}(\omega) &= \Four{ r_{xx}(t) } + &\text{(Wiener-Khinchin)\cite{wiener_khinchin},} \label{eq.wiener_khinchin} +\end{align} +where $r_{xx}(t)$ is defined in terms of the expectation value +\begin{align} + r_{xx}(t) &\equiv \avg{x(\tau)\conj{x}(\tau-t)} + &\text{\cite{wikipedia_wiener_khinchin}} +\end{align} +and $\conj{x}$ represents the complex conjugate of $x$. diff --git a/tex/src/cantilever_calib/solve_general.tex b/tex/src/cantilever_calib/solve_general.tex new file mode 100644 index 0000000..bb4076d --- /dev/null +++ b/tex/src/cantilever_calib/solve_general.tex @@ -0,0 +1,60 @@ +\subsection{General form} + +The procedure here is exactly the same as the previous section. +The integral normalizing $G_0$ just become a little more complicated\ldots + +Fourier transforming eq. \ref{eq.DHO} and applying \ref{eq.four_deriv} we have +\begin{align} + (-m\omega^2 + i \gamma \omega + k) x(\omega) &= F(\omega) + \label{eq.DHO_freq} \\ + (\omega_0^2-\omega^2 + i \beta \omega) x(\omega) &= \frac{F(\omega)}{m} \\ + |x(\omega)|^2 &= \frac{|F(\omega)|^2/m^2} + {(\omega_0^2-\omega^2)^2 + \beta^2\omega^2} + \label{eq.DHO_xmag}, +\end{align} +where $\omega_0 \equiv \sqrt{k/m}$ is the resonant angular frequency + and $\beta \equiv \gamma / m$ is the drag-acceleration coefficient. + +We compute the \PSD\ by plugging eq. \ref{eq.DHO_xmag} into \ref{eq.psd_def} +\begin{equation} + \PSD(x, \omega) + = \normLimT \frac{2 |F(\omega)|^2/m^2} + {(\omega_0^2-\omega^2)^2 + \beta^2\omega^2}. + \label{eq.DHO_psd} +\end{equation} + +Plugging eq. \ref {eq.GOdef} into \ref{eq.DHO_psd} we have +\begin{equation} + \PSD(x, \omega) = \frac{G_0/m^2}{(\omega_0^2-\omega^2)^2 + \beta^2\omega^2}. + \label{eq.model_psd} +\end{equation} +Integrating over positive $\omega$ to find the total power per unit time yields +\begin{align} + \iOInf{\omega}{\PSD(x, \omega)} + &= \frac{G_0}{2m^2} + \iInfInf{\omega}{\frac{1}{(\omega_0^2-\omega^2)^2 + \beta^2\omega^2}} \\ + &= \frac{G_0}{2m^2} \cdot \frac{\pi}{\beta\omega_0^2} \\ + &= \frac{G_0 \pi}{2m^2\beta\omega_0^2} \\ + &= \frac{G_0 \pi}{2m^2\beta \frac{k}{m}} \\ + &= \frac{G_0 \pi}{2m \beta k} +\end{align} +The integration is detailed in Section \ref{sec.integrals}. +By our corollary to Parseval's theorem (eq. \ref{eq.parseval_var}), we have +\begin{equation} + \avg{x(t)^2} = \frac{G_0 \pi}{2m^2\beta\omega_0^2} \label{eq.DHO_var} +\end{equation} + +Plugging eq. \ref{eq.DHO_var} into the equipartition theorem +(eqn. \ref{eq.equipart}) we have +\begin{align} + k \frac{G_0 \pi}{2m \beta k} &= k_BT \\ + G_0 &= \frac{2}{\pi} k_BT m \beta. \label{eq.GO} +\end{align} + +So we expect $x(t)$ to have a power spectral density per unit time given by +\begin{equation} + \PSD(x, \omega) = \frac{2 k_BT \beta} + {\pi m \left[ + (\omega_0^2-\omega^2)^2 + \beta^2\omega^2 + \right] } +\end{equation} diff --git a/tex/src/cantilever_calib/solve_highly_damped.tex b/tex/src/cantilever_calib/solve_highly_damped.tex new file mode 100644 index 0000000..2855b14 --- /dev/null +++ b/tex/src/cantilever_calib/solve_highly_damped.tex @@ -0,0 +1,64 @@ +\subsection{Highly damped case} + +For highly damped systems, the inertial term becomes insignificant + ($m \rightarrow 0$). +This model is commonly used for optically trapped beads. % \cite{} +Because it is simpler and solutions are more easily available, + %cite{grossman05}{}{}{} +we'll use it outline the general approach before diving into the general case. + +Fourier transforming eq. \ref{eq.DHO} with $m=0$ and applying \ref{eq.four_deriv} we have +% ODHO stands for very Over Damped Harmonic oscillator +\begin{align} + (i \gamma \omega + k) x(\omega) &= F(\omega) \label{eq.ODHO_freq} \\ + |x(\omega)|^2 &= \frac{|F(\omega)|^2}{k^2 + \gamma^2 \omega^2}. + \label{eq.ODHO_xmag} +\end{align} +We compute the \PSD\ by plugging eq. \ref{eq.ODHO_xmag} into \ref{eq.psd_def} +\begin{equation} + \PSD(x, \omega) + = \normLimT \frac{2\magSq{F(\omega)}}{k^2 + \gamma^2\omega^2}. + \label{eq.ODHO_psd} +\end{equation} + +Because thermal noise is white (not autocorrelated + Wiener-Khinchin Theorem), +we can denote the one sided thermal power spectral density per unit time by +\begin{equation} + \PSD(F, \omega) = G_0 + = \normLimT 2 \magSq{F(\omega)} \label{eq.GOdef} % label O != zero +\end{equation} + +Plugging eq. \ref {eq.GOdef} into \ref{eq.ODHO_psd} we have +\begin{equation} + \PSD(x, \omega) = \frac{G_0}{k^2 + \gamma^2\omega^2}. +\end{equation} +This is the formula we would use to fit our measured \PSD, but let us go a +bit farther to find the expected \PSD\ and thermal noise + given $m$, $\gamma$ and $k$. + +Integrating over positive $\omega$ to find the total power per unit time yields +\begin{align} + \iOInf{\omega}{\PSD(x, \omega)} + &= \iOInf{\omega}{\frac{G_0}{k^2 + \gamma^2\omega^2}} \\ + &= \frac{G_0}{\gamma}\iOInf{z}{\frac{1}{k^2 + z^2}} \\ + &= \frac{G_0 \pi}{2 \gamma k}, +\end{align} +where the integral is solved in Section \ref{sec.integrals}. + +Plugging into our corollary to Parseval's theorem (eq. \ref{eq.parseval_var}), +\begin{equation} + \avg{x(t)^2} = \frac{G_0 \pi}{2 \gamma k} \label{eq.ODHO_var} +\end{equation} + +Plugging eq. \ref{eq.ODHO_var} into eqn. \ref{eq.equipart} we have +\begin{align} + k \frac{G_0 \pi}{2 \gamma k} &= k_BT \\ + G_0 &= \frac{2 \gamma k_BT}{\pi}. +\end{align} + +So we expect $X(t)$ to have a power spectral density per unit time given by +\begin{equation} + \PSD(x, \omega) = \frac{2}{\pi} + \cdot + \frac{\gamma k_BT}{k^2 + \gamma^2\omega^2}. +\end{equation} diff --git a/tex/src/cantilever_calib/test/Lorentz_root.py b/tex/src/cantilever_calib/test/Lorentz_root.py new file mode 100644 index 0000000..624bfbd --- /dev/null +++ b/tex/src/cantilever_calib/test/Lorentz_root.py @@ -0,0 +1,44 @@ +from numpy import sqrt, complex + +w = 1 + +def f(g) : + g = complex(g) + if g == 0 : + groupAp = 0 + groupAm = 0 + else : + rootA = sqrt(1.0-4.0*w**2/complex(g**2)) + groupAp = g**2 /2.0 * (1+ rootA) + groupAm = g**2 /2.0 * (1- rootA) + rootBp = sqrt(w**2 - groupAp) + rootBm = sqrt(w**2 - groupAm) + return (rootBp, rootBm) + +if __name__ == "__main__" : + import matplotlib + matplotlib.use('Agg') # select backend that doesn't require X Windows + import pylab + from numpy import linspace, array + + gammas = linspace(0,4*w,100) + roots = [] + for g in gammas : + Bp, Bm = f(g) + roots.append( [Bp, Bm, -Bp, -Bm] ) + roots = array(roots) + print roots.shape + pylab.figure() + pylab.subplot(311) + pylab.plot(roots[:,0].real, roots[:,0].imag, 'r+') + pylab.plot(roots[:,1].real, roots[:,1].imag, 'bx') + pylab.plot(roots[:,2].real, roots[:,2].imag, 'k.') + pylab.plot(roots[:,3].real, roots[:,3].imag, 'k.') + pylab.subplot(312) + pylab.plot(gammas, roots[:,0].real, 'r+') + pylab.plot(gammas, roots[:,1].real, 'bx') + pylab.subplot(313) + pylab.plot(gammas, roots[:,0].imag, 'r+') + pylab.plot(gammas, roots[:,1].imag, 'bx') + pylab.show() + diff --git a/tex/src/cantilever_calib/test/Makefile b/tex/src/cantilever_calib/test/Makefile new file mode 100644 index 0000000..386a300 --- /dev/null +++ b/tex/src/cantilever_calib/test/Makefile @@ -0,0 +1,9 @@ +PYTHON_FILES = Lorentz_root locate_roots test_integral test_roots + +all : $(PYTHON_FILES:%=%.png) + +%.png : %.py + python $< + +clean : + rm -f $(PYTHON_FILES:%=%.png) diff --git a/tex/src/cantilever_calib/test/locate_roots.py b/tex/src/cantilever_calib/test/locate_roots.py new file mode 100644 index 0000000..c4b6d39 --- /dev/null +++ b/tex/src/cantilever_calib/test/locate_roots.py @@ -0,0 +1,51 @@ +import matplotlib +matplotlib.use('Agg') # select backend that doesn't require X Windows +from pylab import * +from numpy import * + +N=100 +root_height = 1 +i = complex(0,1) + +def zr(a,b,p,q) : # p, q \in {-1,1} + return p*i*b/2.0 * (1 + q*sqrt(complex(1,0) - (2*a/b)**2)) +def zrc(c,b,p,q) : # p, q \in {-1,1} + return p*i*b/2.0 * (1 + q*sqrt(complex(1,0) - (2*c)**2)) +def loc_roots(c,b) : + figure() + subplot(311) + hold(True) + dpp = zrc(c,b,+1,+1) + dmp = zrc(c,b,-1,+1) + dpm = zrc(c,b,+1,-1) + dmm = zrc(c,b,-1,-1) + plot(dpp.real,dpp.imag,'ro', label="++") + plot(dmp.real,dmp.imag,'b+', label="-+") + plot(dpm.real,dpm.imag,'gx', label="+-") + plot(dmm.real,dmm.imag,'k.', label="--") + title('plane') + legend(loc='upper right') + subplot(312) + plot(c, dpp.real,'ro', + c, dmp.real,'b+', + c, dpm.real,'gx', + c, dmm.real,'k.') + title('real') + subplot(313) + plot(c, dpp.imag,'ro', + c, dmp.imag,'b+', + c, dpm.imag,'gx', + c, dmm.imag,'k.') + title('imag') + +ioff() + +# critically damped for b^2 = 4a^2 +# so for c = a/b = 0.5 +#for b in [0.1, 1, 2, 4, 6] : +# testroots(a,b,x) +c = linspace(0,4,N-1) +loc_roots(c,1) + + FIGURE.savefig("figure.png") +show() diff --git a/tex/src/cantilever_calib/test/test_integral.py b/tex/src/cantilever_calib/test/test_integral.py new file mode 100644 index 0000000..93ce2da --- /dev/null +++ b/tex/src/cantilever_calib/test/test_integral.py @@ -0,0 +1,80 @@ +import matplotlib +matplotlib.use('Agg') # select backend that doesn't require X Windows +from pylab import * +from numpy import * + +N=1000 +root_height = 1 +i = complex(0,1) + +def d(a,b,x) : # know denominator + return (x, ((a**2 - x**2)**2 + (b*x)**2)) +def rm(r, height) : # root marker + return ([r, r], float(height) * array([0.0001,1])) +def zr(a,b,p,q) : # p, q \in {-1,1} + return p*i*b/2.0 * (1 + q*sqrt(complex(1,0) - (2*a/b)**2)) +def zrc(c,b,p,q) : # p, q \in {-1,1} + return p*i*b/2.0 * (1 + q*sqrt(complex(1,0) - (2*c)**2)) +def dr(a,b,x) : # attempt at factoring + c = float(a)/b + return (x, (x-zrc(c,b,+1,+1)) \ + *(x-zrc(c,b,-1,+1)) \ + *(x-zrc(c,b,+1,-1)) \ + *(x-zrc(c,b,-1,-1)) ) + +def numeric(a,b,x) : + (x, dt) = d(a,b,x) + dx = x[1]-x[0] + return (1/dt).mean()*(x[-1]-x[0])*2 # 2 because we want int from -inf to inf +def theory_residue_thm(a,b) : + c = complex(a,0)/b + zpp = zrc(c,b,+1,+1) + zmp = zrc(c,b,-1,+1) + zpm = zrc(c,b,+1,-1) + zmm = zrc(c,b,-1,-1) + return 2*pi*i* ( 1/((zpp-zpm)*(zpp-zmp)*(zpp-zmm)) + 1/((zpm-zpp)*(zpm-zmp)*(zpm-zmm)) ) +def theory_before_plugging_in_roots(a,b) : + c = complex(a,0)/b + zpp = zrc(c,b,+1,+1) + zmp = zrc(c,b,-1,+1) + zpm = zrc(c,b,+1,-1) + zmm = zrc(c,b,-1,-1) + return pi*i/(zpp**2 - zpm**2) * ( 1/zpp - 1/zpm ) +def theory_before_mc(a,b) : + c = complex(a,0)/b + zpp = zrc(c,b,+1,+1) + zmp = zrc(c,b,-1,+1) + zpm = zrc(c,b,+1,-1) + zmm = zrc(c,b,-1,-1) + sq = sqrt(1-4*c**2) + return (-4*pi/b**2) \ + / ((1+sq)**2 - (1-sq)**2) * ((1-sq) - (1+sq)) \ + /(b/2*(1+sq)*(1-sq)) +def theory_final(a,b) : + return pi/(complex(b)* a**2) + +theory = theory_final + +PLOT_COMPARISONS = False + +def test(a) : + x = linspace(0,8*a,N-1) + nA = [] ; tA = [] + if PLOT_COMPARISONS == True : + figure() + for b in linspace(a/10.0, 8*a, 500) : + if PLOT_COMPARISONS == True : + (x1, dt) = d(a,b,x) + plot(x1, 1/dt, 'r.') + (x1, dt) = dr(a,b,x) + plot(x1, 1/dt, 'b-') + nA.append(numeric(a,b,x)) + tA.append(theory(a,b)) + figure() + plot(nA, 'r.', tA, 'b-') + +test(10) +test(1) +test(0.1) + +show() diff --git a/tex/src/cantilever_calib/test/test_roots.py b/tex/src/cantilever_calib/test/test_roots.py new file mode 100644 index 0000000..88c1026 --- /dev/null +++ b/tex/src/cantilever_calib/test/test_roots.py @@ -0,0 +1,65 @@ +import matplotlib +matplotlib.use('Agg') # select backend that doesn't require X Windows +from pylab import * +from numpy import * + +N=100 +root_height = 1 +i = complex(0,1) + +def d(a,b,x) : # know denominator + return (x, ((a**2 - x**2)**2 + (b*x)**2)) +def rm(r, height) : # root marker + return ([r, r], float(height) * array([0.0001,1])) +def zr(a,b,p,q) : # p, q \in {-1,1} + return p*i*b/2.0 * (1 + q*sqrt(complex(1,0) - (2*a/b)**2)) +def rp(a,b,p,q,h) : # plot root + (x,d) = rm(zr(a,b,p,q), h) + if x[0].imag == 0 and x[0] > 0 : plot(x,d,'-x') + print "Root at ", x[0] +def dr(a,b,x) : # attempt at factoring + return (x, (x-zr(a,b,+1,+1)) \ + *(x-zr(a,b,-1,+1)) \ + *(x-zr(a,b,+1,-1)) \ + *(x-zr(a,b,-1,-1)) ) +def testroots(a,b,x) : + figure() + hold(True) + (x,dt) = d(a, b, x) + plot(x,dt,'.') + print "For a = %g, b = %g" % (a, b) + rp(a,b,+1,+1, root_height) + rp(a,b,-1,+1, root_height*2) + rp(a,b,+1,-1, root_height*3) + rp(a,b,-1,-1, root_height*4) + title("For a = %g, b = %g" % (a, b)) +def testfn(a,b,x) : + figure() + hold(True) + (x,dt) = d(a, b, x) + plot(x,1/dt,'-') + (x,dt) = dr(a, b, x) + plot(x,1/dt,'.') + title("For a = %g, b = %g" % (a, b)) + +ioff() + +plot = semilogy + +# critically damped for b^2 = 4a^2 +#for b in [0.1, 1, 2, 4, 6] : +# testroots(a,b,x) +a = 10.0 +x = linspace(0,4*a,N-1) +for b in linspace(a/10, 4*a, 5) : + testfn(a,b,x) +a = 1.0 +x = linspace(0,4*a,N-1) +for b in linspace(a/10, 4*a, 5) : + testfn(a,b,x) +a = .1 +x = linspace(0,4*a,N-1) +for b in linspace(a/10, 4*a, 5) : + testfn(a,b,x) + +show() diff --git a/tex/src/local_cmmds.tex b/tex/src/local_cmmds.tex new file mode 100644 index 0000000..e713b13 --- /dev/null +++ b/tex/src/local_cmmds.tex @@ -0,0 +1,59 @@ +% Fourier Transform to angular momentum space +\newcommand{\Four}[1]{\ensuremath{{\mathcal F}\left\{ {#1} \right\}}} +% Fourier Transform to frequency space +\newcommand{\Fourf}[1]{\ensuremath{{\mathcal F}_f\left\{ {#1} \right\}}} + +% Integral from #1 to #2 of #4 with respect to #3 +\newcommand{\integral}[4]{\ensuremath{\int_{#1}^{#2} {#4} \dd{#3}}} +% Integral from -infty to +infty of #2 with respect to #1 +\newcommand{\iInfInf}[2]{\integral{-\infty}{\infty}{#1}{#2}} +% Integral from 0 to +infty of #2 with respect to #1 +% note that the second char in the name is a capital o, not a 0 +% because ?Latex doesn't allow digits in command names +\newcommand{\iOInf}[2]{\integral{0}{\infty}{#1}{#2}} +% #3 evaluated at #2 and #1 +\newcommand{\evaluated}[3]{\ensuremath{\left.{#3}\right|_{#1}^{#2}}} +% #2 evaluated at #1 +\newcommand{\eval}[2]{\ensuremath{\left.{#2}\right|_{#1}}} +% evaluated at infty and -infty +\newcommand{\eInfInf}[1]{\evaluated{-\infty}{\infty}{#1}} +% evaluated at infty and 0 +\newcommand{\eOInf}[1]{\evaluated{0}{\infty}{#1}} + +% we do a lot of lim t_T -> infty, so macro that out +% lim as #1 -> infty +\newcommand{\limInf}[1]{\ensuremath{\lim_{{#1}\rightarrow \infty}}} +\newcommand{\limT}{\limInf{t_T}} +\newcommand{\normLimT}{\limT \frac{1}{t_T}\ } % '\ ' for a normal space +% lim as t_T -> infty of integral from -t_T/2 to t_T/2 of #1 with respect to t +\newcommand{\iLimT}[1]{\normLimT \integral{-t_T/2}{t_T/2}{t}{#1}} + +\newcommand{\magSq}[1]{\ensuremath{\left|{#1}\right|^2}} +\newcommand{\PSD}{\ensuremath{\operatorname{PSD}}} + +% complex conjugate +\newcommand{\conj}[1]{\ensuremath{\overline{#1}}} + +% Symbol denoting the contour C +\newcommand{\C}{\ensuremath{\mathcal{C}}} +% Integral about contour C of #1 with respect to z +\newcommand{\iC}[1]{\integral{\C}{}{z}{#1}} +\newcommand{\Reals}{\ensuremath{\mathbb{R}}} +\newcommand{\Imags}{\ensuremath{\mathbb{I}}} +\newcommand{\Real}{\ensuremath{\operatorname{Re}}} +\newcommand{\Imag}{\ensuremath{\operatorname{Im}}} +\newcommand{\Res}[2]{\operatorname{Res}\left({z=#1},{#2}\right)} +\newcommand{\limZ}[1]{\lim_{z \rightarrow {#1}}} +\newcommand{\limZp}{\limZ{z_p}} +\newcommand{\CPV}{\ensuremath{\mathbb{P}}} + +\newcommand{\supScript}[1]{\ensuremath{^{\mathrm{#1}}}} +\newcommand{\st}{\supScript{st}} +\newcommand{\nd}{\supScript{nd}} +\newcommand{\rd}{\supScript{rd}} +\newcommand{\sth}{\supScript{th}} % th, TH already taken + +\newcommand{\colA}[1]{\textcolor{Red}{#1}} +\newcommand{\colB}[1]{\textcolor{Blue}{#1}} +\newcommand{\colAB}[1]{\textcolor{Purple}{#1}} +\newcommand{\colC}[1]{\textcolor{Green}{#1}} diff --git a/tex/src/packages.tex b/tex/src/packages.tex new file mode 100644 index 0000000..7cda711 --- /dev/null +++ b/tex/src/packages.tex @@ -0,0 +1,25 @@ +\usepackage[super,sort&compress]{natbib} % fancy citation extensions +% super selects citations in superscript mode +% sort&compress automatically sorts and compresses compound citations (\cite{a,b,...}) + +%\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet +%\bibliographystyle{plain} % pick the bibliography style, includes dates +\bibliographystyle{plainnat} + + +% print labels in margins for debugging +%\usepackage{showkeys} + +% environments for multiline displayed equations, and other enhancements +\usepackage{amsmath} +% proof environment and extensions for the \newtheorem command +\usepackage{amsthm} +% fonts like mathbb (blackboard) +\usepackage{amsfonts} + +\usepackage{graphicx} % to include images +% allow colored text (usenames -> with 68 dvips named colors) +\usepackage[dvipsnames, usenames]{color} + +\usepackage{wtk_cmmds} % common personal macros, in ~/texmf/tex/latex/ +\input{local_cmmds} diff --git a/tex/src/root.tex b/tex/src/root.tex new file mode 100644 index 0000000..d7925ee --- /dev/null +++ b/tex/src/root.tex @@ -0,0 +1,54 @@ +\documentclass[% +% draft% (remove the double spacing) + ]{thesis} +% See thesis.cls for more options. + +%\includeonly{% +% } + +\author{William Trevor King} % Fullname +\title{Temperature dependent mechanical protein unfolding} % Title of thesis +\month{July} % Name of the month of you defense +\year{2010} % Year you are defending (YYYY) +\degree{Doctor of Philosophy} % Should be fine +\advisor{Guoliang Yang, Ph.D.} % Put advisor's full name, degree + +\input{packages} + +\begin{document} +\begin{preamble} + +\begin{dedications} +% Type dedications in here - OPTIONAL +\end{dedications} + +\begin{acknowledgments} +% Type acknowledgments in here - OPTIONAL +\end{acknowledgments} + +\tableofcontents +% Include these following commands only if you have tables or figures. +% Tables should come before figures! +\listoftables +\listoffigures + +\begin{abstract} +% Type abstract in here +\end{abstract} +\end{preamble} + +\begin{thesis} +% Use include statements to include your main thesis text +% from seperate files. \include{chapter1} etc... +\end{thesis} + +\bibliography{wtk} + +\appendix +\include{cantilever_calib/cantilever_calib} + +\begin{vita} +% Ph.D. only. See manual for details. +\end{vita} + +\end{document} diff --git a/tex/src/wtk.bib b/tex/src/wtk.bib new file mode 120000 index 0000000..7cfefbc --- /dev/null +++ b/tex/src/wtk.bib @@ -0,0 +1 @@ +/home/wking/rsrch/bib/wtk.bib \ No newline at end of file -- 2.26.2