Use relative imports (PEP 328) for calibcant sibling imports.
[calibcant.git] / calibcant / calibrate.py
1 #!/usr/bin/python
2 #
3 # calibcant - tools for thermally calibrating AFM cantilevers
4 #
5 # Copyright (C) 2007-2009 William Trevor King
6 #
7 # This program is free software; you can redistribute it and/or
8 # modify it under the terms of the GNU General Public License as
9 # published by the Free Software Foundation; either version 3 of the
10 # License, or (at your option) any later version.
11 #
12 # This program is distributed in the hope that it will be useful, but
13 # WITHOUT ANY WARRANTY; without even the implied warranty of
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
15 # See the GNU General Public License for more details.
16 #
17 # You should have received a copy of the GNU General Public License
18 # along with this program; if not, write to the Free Software
19 # Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
20 # 02111-1307, USA.
21 #
22 # The author may be contacted at <wking@drexel.edu> on the Internet, or
23 # write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
24 # Philadelphia PA 19104, USA.
25
26 """
27 Aquire and analyze cantilever calibration data.
28
29 W. Trevor King Dec. 2007-Jan. 2008
30
31 GPL BOILERPLATE
32
33
34 The relevent physical quantities are :
35  Vzp_out  Output z-piezo voltage (what we generate)
36  Vzp      Applied z-piezo voltage (after external ZPGAIN)
37  Zp       The z-piezo position
38  Zcant    The cantilever vertical deflection
39  Vphoto   The photodiode vertical deflection voltage (what we measure)
40  Fcant    The force on the cantilever
41  T        The temperature of the cantilever and surrounding solution
42           (another thing we measure or guess)
43  k_b      Boltzmann's constant
44
45 Which are related by the parameters :
46  zpGain           Vzp_out / Vzp
47  zpSensitivity    Zp / Vzp
48  photoSensitivity Vphoto / Zcant
49  k_cant           Fcant / Zcant
50
51 Cantilever calibration assumes a pre-calibrated z-piezo
52 (zpSensitivity) and a amplifier (zpGain).  In our lab, the z-piezo is
53 calibrated by imaging a calibration sample, which has features with
54 well defined sizes, and the gain is set with a knob on the Nanoscope.
55
56 photoSensitivity is measured by bumping the cantilever against the
57 surface, where Zp = Zcant (see bump_aquire() and the bump_analyze
58 submodule).
59
60 k_cant is measured by watching the cantilever vibrate in free solution
61 (see the vib_aquire() and the vib_analyze submodule).  The average
62 energy of the cantilever in the vertical direction is given by the
63 equipartition theorem.
64     1/2 k_b T   =   1/2 k_cant <Zcant**2>
65  so     k_cant  = k_b T / Zcant**2
66  but    Zcant   = Vphoto / photoSensitivity
67  so     k_cant  = k_b T * photoSensitivty**2 / <Vphoto**2>
68
69 We measured photoSensitivity with the surface bumps.  We can either
70 measure T using an external function (see temperature.py), or just
71 estimate it (see T_aquire() and the T_analyze submodule).  Guessing
72 room temp ~22 deg C is actually fairly reasonable.  Assuming the
73 actual fluid temperature is within +/- 5 deg, the error in the spring
74 constant k_cant is within 5/273.15 ~= 2%.  A time series of Vphoto
75 while we're far from the surface and not changing Vzp_out will give us
76 the average variance <Vphoto**2>.
77
78 We do all these measurements a few times to estimate statistical
79 errors.
80
81 The functions are layed out in the families:
82  bump_*(), vib_*(), T_*(), and calib_*()
83 where calib_{save|load|analyze}() deal with derived data, not
84 real-world data.
85
86 For each family, * can be any of :
87  aquire       get real-world data
88  save         store real-world data to disk
89  load         get real-world data from disk
90  analyze      interperate the real-world data.
91  plot         show a nice graphic to convince people we're working :p
92  load_analyze_tweaked
93               read a file with a list of paths to previously saved
94               real world data load each file using *_load(), analyze
95               using *_analyze(), and optionally plot using *_plot().
96               Intended for re-processing old data.
97 A family name without any _* extension (e.g. bump()), runs *_aquire(),
98  *_save(), *_analyze(), *_plot().
99
100 We also define the two positioning functions:
101  move_just_onto_surface() and move_far_from_surface()
102 which make automating the calibration procedure more straightforward.
103 """
104
105 import numpy
106 import time 
107
108 import FFT_tools
109 import piezo.z_piezo_utils as z_piezo_utils
110 from splittable_kwargs import splittableKwargsFunction, \
111     make_splittable_kwargs_function
112
113 import .common
114 import .config
115 import .bump_analyze
116 import .T_analyze
117 import .vib_analyze
118 import .analyze
119
120
121 # bump family
122
123 @splittableKwargsFunction()
124 def bump_aquire(zpiezo, push_depth=200, npoints=1024, push_speed=1000) :
125     """
126     Ramps closer push_depth and returns to the original position.
127     Inputs:
128      zpiezo     an opened zpiezo.zpiezo instance
129      push_depth distance to approach, in nm
130      npoints    number points during the approach and during the retreat
131      push_speed piezo speed during approach and retreat, in nm/s
132     Returns the aquired ramp data dictionary, with data in DAC/ADC bits.
133     """
134     # generate the bump output
135     nm_per_step = float(push_depth) / npoints
136     freq = push_speed / nm_per_step # freq is sample frequency in Hz
137     start_pos = zpiezo.curPos()
138     pos_dist = zpiezo.pos_nm2out(push_depth) - zpiezo.pos_nm2out(0)
139     close_pos = start_pos + pos_dist
140     appr = numpy.linspace(start_pos, close_pos, npoints)
141     retr = numpy.linspace(close_pos, start_pos, npoints)
142     out = numpy.concatenate((appr, retr))
143     # run the bump, and measure deflection
144     if config.TEXT_VERBOSE :
145         print "Bump %g nm at %g nm/s" % (push_depth, push_speed)
146     data = zpiezo.ramp(out, freq)
147     return data
148
149 @splittableKwargsFunction(bump_aquire,
150                           (bump_analyze.bump_save, 'data'),
151                           (bump_analyze.bump_analyze, 'data'))
152 def bump(**kwargs):
153     """
154     Wrapper around bump_aquire(), bump_save(), bump_analyze()
155     """
156     bump_aquire_kwargs,bump_save_kwargs,bump_analyze_kwargs = \
157         bump._splitargs(bump, kwargs)
158     data = bump_aquire(**bump_aquire_kwargs)
159     bump_analyze.bump_save(data, **bump_save_kwargs)
160     photoSensitivity = bump_analyze.bump_analyze(data, **bump_analyze_kwargs)
161     return photoSensitivity
162
163 # T family.
164 # Fairly stubby, since a one shot Temp measurement is a common thing.
165 # We just wrap that to provide a consistent interface.
166
167 @splittableKwargsFunction()
168 def T_aquire(get_T=None) :
169     """
170     Measure the current temperature of the sample, 
171     or, if get_T == None, fake it by returning config.DEFAULT_TEMP
172     """
173     if get_T == None :
174         if config.TEXT_VERBOSE :
175             print "Fake temperature %g" % config.DEFAULT_TEMP
176         return config.DEFAULT_TEMP
177     else :
178         if config.TEXT_VERBOSE :
179             print "Measure temperature"
180         return get_T()
181
182 @splittableKwargsFunction(T_aquire,
183                           (T_analyze.T_save, 'T'),
184                           (T_analyze.T_analyze, 'T'))
185 def T(**kwargs):
186     """
187     Wrapper around T_aquire(), T_save(), T_analyze(), T_plot()
188     """
189     T_aquire_kwargs,T_save_kwargs,T_analyze_kwargs = \
190         T._splitargs(T, kwargs)
191     T_raw = T_aquire(**T_aquire_kwargs)
192     T_analyze.T_save(T_raw, **T_save_kwargs)
193     T_ret = T_analyze.T_analyze(T_raw, **T_analyze_kwargs) # returns array
194     return T_ret[0]
195
196 # vib family
197
198 @splittableKwargsFunction()
199 def vib_aquire(zpiezo, time=1, freq=50e3) :
200     """
201     Record data for TIME seconds at FREQ Hz from ZPIEZO at it's current position.
202     """
203     # round up to the nearest power of two, for efficient FFT-ing
204     nsamps = FFT_tools.ceil_pow_of_two(time*freq)
205     time = nsamps / freq
206     # take some data, keeping the position voltage at it's current value
207     out = numpy.ones((nsamps,), dtype=numpy.uint16) * zpiezo.curPos()
208     if config.TEXT_VERBOSE :
209         print "get %g seconds of data" % time
210     data = zpiezo.ramp(out, freq)
211     data['sample frequency Hz'] = numpy.array([freq])
212     return data
213
214 @splittableKwargsFunction(vib_aquire,
215                           (vib_analyze.vib_save, 'data'),
216                           (vib_analyze.vib_analyze, 'deflection_bits', 'freq'))
217 def vib(**kwargs) :
218     """
219     Wrapper around vib_aquire(), vib_save(), vib_analyze()
220     """
221     vib_aquire_kwargs,vib_save_kwargs,vib_analyze_kwargs = \
222         vib._splitargs(vib, kwargs)
223     data = vib_aquire(**vib_aquire_kwargs)
224     vib_analyze.vib_save(data, **vib_save_kwargs)
225     freq = data['sample frequency Hz']
226     deflection_bits = data['Deflection input']
227     Vphoto_var = vib_analyze.vib_analyze(deflection_bits=deflection_bits,
228                                          freq=freq, **vib_analyze_kwargs)
229     return Vphoto_var
230
231 # A few positioning functions, so we can run bump_aquire() and vib_aquire()
232 # with proper spacing relative to the surface.
233
234 @splittableKwargsFunction()
235 def move_just_onto_surface(stepper, zpiezo, Depth_nm=-50, setpoint=2) :
236     """
237     Uses z_piezo_utils.getSurfPos() to pinpoint the position of the
238     surface.  Adjusts the stepper position as required to get within
239     stepper_tol nm of the surface.  Then set Vzp to place the
240     cantilever Depth_nm onto the surface.  Negative Depth_nm values
241     will place the cantilever that many nm _off_ the surface.
242     
243     If getSurfPos() fails to find the surface, backs off (for safety)
244     and steps in (without moving the zpiezo) until Vphoto > setpoint.
245     """
246     stepper_tol = 250 # nm, generous estimate of the fullstep stepsize
247
248     if config.TEXT_VERBOSE :
249         print "moving just onto surface"
250     # Zero the piezo
251     if config.TEXT_VERBOSE :
252         print "zero the z piezo output"
253     zpiezo.jumpToPos(zpiezo.pos_nm2out(0))
254     # See if we're near the surface already
255     if config.TEXT_VERBOSE :
256         print "See if we're starting near the surface"
257     try :
258         dist = zpiezo.pos_out2nm( \
259             z_piezo_utils.getSurfPos(zpiezo, zpiezo.def_V2in(setpoint))
260                                 )
261     except (z_piezo_utils.tooClose, z_piezo_utils.poorFit), string :
262         if config.TEXT_VERBOSE :
263             print "distance failed with: ", string
264             print "Back off 200 half steps"
265         # Back away 200 steps
266         stepper.step_rel(-400)
267         stepper.step_rel(200)
268         sp = zpiezo.def_V2in(setpoint) # sp = setpoint in bits
269         zpiezo.updateInputs()
270         cd = zpiezo.curDef()           # cd = current deflection in bits
271         if config.TEXT_VERBOSE :
272             print "Single stepping approach"
273         while cd < sp :
274             if config.TEXT_VERBOSE :
275                 print "deflection %g < setpoint %g.  step closer" % (cd, sp)
276             stepper.step_rel(2) # Full step in
277             zpiezo.updateInputs()
278             cd = zpiezo.curDef()
279         # Back off two steps (protecting against backlash)
280         if config.TEXT_VERBOSE :
281             print "Step back 4 half steps to get off the setpoint"
282         stepper.step_rel(-200)
283         stepper.step_rel(196)
284         # get the distance to the surface
285         zpiezo.updateInputs()
286         if config.TEXT_VERBOSE :
287             print "get surf pos, with setpoint %g (%d)" % (setpoint, zpiezo.def_V2in(setpoint))
288         for i in range(20) : # HACK, keep stepping back until we get a distance
289             try :
290                 dist = zpiezo.pos_out2nm( \
291                     z_piezo_utils.getSurfPos(zpiezo,zpiezo.def_V2in(setpoint)))
292             except (z_piezo_utils.tooClose, z_piezo_utils.poorFit), string :
293                 stepper.step_rel(-200)
294                 stepper.step_rel(198)
295                 continue
296             break
297         if i >= 19 :
298             print "tried %d times, still too close! bailing" % i
299             print "probably an invalid setpoint."
300             raise Exception, "weirdness"
301     if config.TEXT_VERBOSE :
302         print 'distance to surface ', dist, ' nm'
303     # fine tune the stepper position
304     while dist < -stepper_tol : # step back if we need to
305         stepper.step_rel(-200)
306         stepper.step_rel(198)
307         dist = zpiezo.pos_out2nm( \
308             z_piezo_utils.getSurfPos(zpiezo, zpiezo.def_V2in(setpoint)))
309         if config.TEXT_VERBOSE :
310             print 'distance to surface ', dist, ' nm, step back'
311     while dist > stepper_tol : # and step forward if we need to
312         stepper.step_rel(2)
313         dist = zpiezo.pos_out2nm( \
314             z_piezo_utils.getSurfPos(zpiezo, zpiezo.def_V2in(setpoint)))
315         if config.TEXT_VERBOSE :
316             print 'distance to surface ', dist, ' nm, step closer'
317     # now adjust the zpiezo to place us just onto the surface
318     target = dist + Depth_nm
319     zpiezo.jumpToPos(zpiezo.pos_nm2out(target))
320     # and we're there :)
321     if config.TEXT_VERBOSE :
322         print "We're %g nm into the surface" % Depth_nm
323
324 @splittableKwargsFunction()
325 def move_far_from_surface(stepper, um_back=50) :
326     """
327     Step back a specified number of microns.
328     (uses very rough estimate of step distance at the moment)
329     """
330     step_nm = 100
331     steps = int(um_back*1000/step_nm)
332     print "step back %d steps" % steps
333     stepper.step_rel(-steps)
334
335
336 # and finally, the calib family
337
338 @splittableKwargsFunction((move_just_onto_surface, 'stepper', 'zpiezo'),
339                           (bump, 'zpiezo', 'log_dir', 'Vphoto_in2V'),
340                           (move_far_from_surface, 'stepper'),
341                           (T, 'log_dir'),
342                           (vib, 'zpiezo', 'log_dir', 'Vphoto_in2V'),
343                           (analyze.calib_save, 'bumps','Ts','vibs','log_dir'))
344 def calib_aquire(stepper, zpiezo, num_bumps=10, num_Ts=10, num_vibs=20,
345                  log_dir=config.LOG_DIR, Vphoto_in2V=config.Vphoto_in2V,
346                  **kwargs):
347     """
348     Aquire data for calibrating a cantilever in one function.
349     return (bump, T, vib), each of which is an array.
350     Inputs :
351      stepper       a stepper.stepper_obj for coarse Z positioning
352      zpiezo        a z_piezo.z_piezo for fine positioning and deflection readin
353      num_bumps     number of 'bumps' (see Outputs)
354      num_temps     number of 'Ts' (see Outputs)
355      num_vibs      number of 'vib's (see Outputs)
356      log_dir       directory to log data to.  Default 'None' disables logging.
357      Vphoto_in2V   function to convert photodiode input bits to Volts
358
359      + other kwargs.  Run calib_aquire._kwargs(calib_aquire) to see
360      all options.  Run calib_aquire._childSplittables to see a list
361      of kwarg functions that this function calls.
362
363     Outputs (all are arrays of recorded data) :
364      bumps measured (V_photodiode / nm_tip) proportionality constant
365      Ts    measured temperature (K)
366      vibs  measured V_photodiode variance in free solution
367     """
368     move_just_onto_surface_kwargs,bump_kwargs,move_far_from_surface_kwargs, \
369         T_kwargs,vib_kwargs,calib_save_kwargs = \
370         calib_aquire._splitargs(calib_aquire, kwargs)
371     # get bumps
372     bumps = numpy.zeros((num_bumps,), dtype=numpy.float)
373     for i in range(num_bumps) :
374         move_just_onto_surface(stepper, zpiezo, **move_just_onto_surface_kwargs)
375         bumps[i] = bump(zpiezo=zpiezo, log_dir=log_dir,
376                         Vphoto_in2V=Vphoto_in2V, **bump_kwargs)
377     if config.TEXT_VERBOSE :
378         print bumps
379
380     move_far_from_surface(stepper, **move_far_from_surface_kwargs)
381
382     # get Ts
383     Ts = numpy.zeros((num_Ts,), dtype=numpy.float)
384     for i in range(num_Ts) :
385         Ts[i] = T(**T_kwargs)
386         time.sleep(1) # wait a bit to get an independent temperature measure
387     print Ts
388
389     # get vibs
390     vibs = numpy.zeros((num_vibs,), dtype=numpy.float)
391     for i in range(num_vibs) :
392         vibs[i] = vib(zpiezo=zpiezo, log_dir=log_dir, Vphoto_in2V=Vphoto_in2V,
393                       **vib_kwargs)
394     print vibs
395     
396     analyze.calib_save(bumps, Ts, vibs, log_dir, **calib_save_kwargs)
397     
398     return (bumps, Ts, vibs)
399
400
401 @splittableKwargsFunction( \
402     (calib_aquire, 'log_dir'),
403     (analyze.calib_analyze, 'bumps','Ts','vibs'))
404 def calib(log_dir=config.LOG_DIR, **kwargs) :
405     """
406     Calibrate a cantilever in one function.
407     The I-don't-care-about-the-details black box version :p.
408     return (k, k_s)
409     Inputs:
410      (see calib_aquire()) 
411     Outputs :
412      k    cantilever spring constant (in N/m, or equivalently nN/nm)
413      k_s  standard deviation in our estimate of k
414     Notes :
415      See get_calibration_data() for the data aquisition code
416      See analyze_calibration_data() for the analysis code
417     """
418     calib_aquire_kwargs,calib_analyze_kwargs = \
419         calib._splitargs(calib, kwargs)
420     a, T, vib = calib_aquire(**calib_aquire_kwargs)
421     k,k_s,ps2_m, ps2_s,T_m,T_s,one_o_Vp2_m,one_o_Vp2_s = \
422         analyze.calib_analyze(a, T, vib, **calib_analyze_kwargs)
423     analyze.calib_save_analysis(k, k_s, ps2_m, ps2_s, T_m, T_s,
424                                 one_o_Vp2_m, one_o_Vp2_s, log_dir)
425     return (k, k_s)
426
427     
428