Added latex framework and cantilever_calib appendix.
authorW. Trevor King <wking@drexel.edu>
Wed, 2 Dec 2009 20:22:15 +0000 (15:22 -0500)
committerW. Trevor King <wking@drexel.edu>
Wed, 2 Dec 2009 20:22:15 +0000 (15:22 -0500)
24 files changed:
tex/Makefile [new file with mode: 0644]
tex/src/Makefile [new file with mode: 0644]
tex/src/cantilever_calib/Makefile [new file with mode: 0644]
tex/src/cantilever_calib/README [new file with mode: 0644]
tex/src/cantilever_calib/cantilever_calib.tex [new file with mode: 0644]
tex/src/cantilever_calib/contour_integration.tex [new file with mode: 0644]
tex/src/cantilever_calib/dot/concept_map.dot [new file with mode: 0644]
tex/src/cantilever_calib/integrals.tex [new file with mode: 0644]
tex/src/cantilever_calib/mp/Makefile [new file with mode: 0644]
tex/src/cantilever_calib/mp/contour.mp [new file with mode: 0644]
tex/src/cantilever_calib/mp/sample.tex [new file with mode: 0644]
tex/src/cantilever_calib/overview.tex [new file with mode: 0644]
tex/src/cantilever_calib/setup_general.tex [new file with mode: 0644]
tex/src/cantilever_calib/solve_general.tex [new file with mode: 0644]
tex/src/cantilever_calib/solve_highly_damped.tex [new file with mode: 0644]
tex/src/cantilever_calib/test/Lorentz_root.py [new file with mode: 0644]
tex/src/cantilever_calib/test/Makefile [new file with mode: 0644]
tex/src/cantilever_calib/test/locate_roots.py [new file with mode: 0644]
tex/src/cantilever_calib/test/test_integral.py [new file with mode: 0644]
tex/src/cantilever_calib/test/test_roots.py [new file with mode: 0644]
tex/src/local_cmmds.tex [new file with mode: 0644]
tex/src/packages.tex [new file with mode: 0644]
tex/src/root.tex [new file with mode: 0644]
tex/src/wtk.bib [new symlink]

diff --git a/tex/Makefile b/tex/Makefile
new file mode 100644 (file)
index 0000000..2dbb26c
--- /dev/null
@@ -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 (file)
index 0000000..8e0b458
--- /dev/null
@@ -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 (file)
index 0000000..f1b083b
--- /dev/null
@@ -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 (file)
index 0000000..368f355
--- /dev/null
@@ -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 (file)
index 0000000..ea6b2fb
--- /dev/null
@@ -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 (file)
index 0000000..c8b5de9
--- /dev/null
@@ -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 (file)
index 0000000..7bf5d60
--- /dev/null
@@ -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 (file)
index 0000000..d9b6752
--- /dev/null
@@ -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 (file)
index 0000000..be5e052
--- /dev/null
@@ -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 (file)
index 0000000..ce60879
--- /dev/null
@@ -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 (file)
index 0000000..e8fca4f
--- /dev/null
@@ -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 (file)
index 0000000..2948376
--- /dev/null
@@ -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 (file)
index 0000000..d391a05
--- /dev/null
@@ -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 (file)
index 0000000..bb4076d
--- /dev/null
@@ -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 (file)
index 0000000..2855b14
--- /dev/null
@@ -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 (file)
index 0000000..624bfbd
--- /dev/null
@@ -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 (file)
index 0000000..386a300
--- /dev/null
@@ -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 (file)
index 0000000..c4b6d39
--- /dev/null
@@ -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 (file)
index 0000000..93ce2da
--- /dev/null
@@ -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 (file)
index 0000000..88c1026
--- /dev/null
@@ -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 (file)
index 0000000..e713b13
--- /dev/null
@@ -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 (file)
index 0000000..7cda711
--- /dev/null
@@ -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 (file)
index 0000000..d7925ee
--- /dev/null
@@ -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 (symlink)
index 0000000..7cfefbc
--- /dev/null
@@ -0,0 +1 @@
+/home/wking/rsrch/bib/wtk.bib
\ No newline at end of file