(* 
Mathematica code to implement the calibration method of
Sader et al., Rev. Sci. Instr. Vol. 70, 3967-3969 (1999).

To run:
(1) first evaluate this file;
(2) then type in k[L, b, f, Q], where 
    L is the length in um,
    b is the width in um, 
    f is the resonant frequency in kHz,
    Q is the quality factor.

E.g.: To evaluate k for 397x29um cantilever in Table 1: Type k[397, 29, 17.36, 55.5]
*)

(* Density and viscosity of air *)
(* ---------------------------------------- *)

rho = 1.18;
eta = 1.86 10^(-5);


(* Hydrodynamic functions *)
(* ---------------------------------------- *)

Ren[f_, b_] := N[rho f b^2 / 4 / eta];

omre[x_] := (0.91324 - 0.48274*x + 0.46842*x^2 - 0.12886*x^3 +
     0.044055*x^4 - 0.0035117*x^5 + 0.00069085*x^6)/
   (1 - 0.56964*x + 0.48690*x^2 - 0.13444*x^3 + 0.045155*x^4 -
     0.0035862*x^5 + 0.00069085*x^6);

omimag[x_] := (-0.024134- 0.029256*x + 0.016294*x^2 - 0.00010961*x^3 +
     0.000064577*x^4 - 0.000044510*x^5)/ (1 - 0.59702*x + 0.55182*x^2 -
     0.18357*x^3 + 0.079156*x^4 - 0.014369*x^5 + 0.0028361*x^6);

om[x_] := omre[x] + I omimag[x];

omega[x_] := om[Log[10, x]];

LAMBDA[f_, b_] := omega[Ren[f, b]] ( 1 + 4 I BesselK[1, - I Sqrt[I Ren[f, b]]] / ( Sqrt[I Ren[f, b]] BesselK[0, - I Sqrt[I Ren[f, b]]] ) );


(* Calculation of spring constant *)
(* ---------------------------------------- *)
(* L, b in um. Resonant frequency f in kHz. Quality factor Q *)

k[L_, b_, f_, Q_] := N[0.1906 rho b^2 L 10^(-18) Q Im[LAMBDA[2 Pi 1000 f, b 10^(-6)]] (2 Pi 1000 f)^2];