--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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 > $@
--- /dev/null
+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.
--- /dev/null
+\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}
--- /dev/null
+\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}
--- /dev/null
+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;
+
+}
--- /dev/null
+\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}
--- /dev/null
+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 $<
--- /dev/null
+%% 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
--- /dev/null
+\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}
--- /dev/null
+\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}
--- /dev/null
+\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$.
--- /dev/null
+\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}
--- /dev/null
+\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}
--- /dev/null
+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()
+
--- /dev/null
+PYTHON_FILES = Lorentz_root locate_roots test_integral test_roots
+
+all : $(PYTHON_FILES:%=%.png)
+
+%.png : %.py
+ python $<
+
+clean :
+ rm -f $(PYTHON_FILES:%=%.png)
--- /dev/null
+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()
--- /dev/null
+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()
--- /dev/null
+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()
--- /dev/null
+% 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}}
--- /dev/null
+\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}
--- /dev/null
+\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}
--- /dev/null
+/home/wking/rsrch/bib/wtk.bib
\ No newline at end of file