1 # calibcant - tools for thermally calibrating AFM cantilevers
3 # Copyright (C) 2008-2011 W. Trevor King <wking@drexel.edu>
5 # This file is part of calibcant.
7 # calibcant is free software: you can redistribute it and/or
8 # modify it under the terms of the GNU Lesser General Public
9 # License as published by the Free Software Foundation, either
10 # version 3 of the License, or (at your option) any later version.
12 # calibcant is distributed in the hope that it will be useful,
13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 # GNU Lesser General Public License for more details.
17 # You should have received a copy of the GNU Lesser General Public
18 # License along with calibcant. If not, see
19 # <http://www.gnu.org/licenses/>.
21 """Calculate `k` from arrays of bumps, temperatures, and vibrations.
23 Separate the more general `calib_analyze()` from the other `calib_*()`
24 functions in calibcant.
26 The relevent physical quantities are :
27 Vzp_out Output z-piezo voltage (what we generate)
28 Vzp Applied z-piezo voltage (after external ZPGAIN)
29 Zp The z-piezo position
30 Zcant The cantilever vertical deflection
31 Vphoto The photodiode vertical deflection voltage (what we measure)
32 Fcant The force on the cantilever
33 T The temperature of the cantilever and surrounding solution
34 (another thing we measure)
35 k_b Boltzmann's constant
37 Which are related by the parameters:
39 zp_sensitivity Zp / Vzp
40 photo_sensitivity Vphoto / Zcant
47 >>> from h5config.storage.hdf5 import HDF5_Storage, pprint_HDF5
48 >>> from .config import CalibrationConfig
50 >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
53 >>> calibration_config = CalibrationConfig(storage=HDF5_Storage(
54 ... filename=filename, group='/calib/config/'))
55 >>> bumps = numpy.array((15.9e6, 16.9e6, 16.3e6))
56 >>> temperatures = numpy.array((295, 295.2, 294.8))
57 >>> vibrations = numpy.array((2.20e-5, 2.22e-5, 2.21e-5))
59 >>> k,k_s = calib_analyze(bumps=bumps, temperatures=temperatures,
60 ... vibrations=vibrations)
61 >>> (k, k_s) # doctest: +ELLIPSIS
62 (0.0493..., 0.00248...)
64 Most of the error in this example comes from uncertainty in the
65 photodiode sensitivity (bumps).
67 >>> k_s/k # doctest: +ELLIPSIS
69 >>> bumps.std()/bumps.mean() # doctest: +ELLIPSIS
71 >>> temperatures.std()/temperatures.mean() # doctest: +ELLIPSIS
73 >>> vibrations.std()/vibrations.mean() # doctest: +ELLIPSIS
76 >>> calib_save(filename=filename, group='/calib/',
77 ... bumps=bumps, temperatures=temperatures, vibrations=vibrations,
78 ... calibration_config=calibration_config, k=k, k_s=k_s)
79 >>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
83 <HDF5 dataset "bump": shape (), type "|S1">
85 <HDF5 dataset "num-bumps": shape (), type "<i4">
87 <HDF5 dataset "num-temperatures": shape (), type "<i4">
89 <HDF5 dataset "num-vibrations": shape (), type "<i4">
91 <HDF5 dataset "temperature": shape (), type "|S1">
93 <HDF5 dataset "temperature-sleep": shape (), type "<i4">
95 <HDF5 dataset "vibration": shape (), type "|S1">
97 <HDF5 dataset "vibration-spacing": shape (), type "<f8">
100 /calib/processed/spring-constant
101 <HDF5 dataset "data": shape (), type "<f8">
103 <HDF5 dataset "standard-deviation": shape (), type "<f8">
105 <HDF5 dataset "units": shape (), type "|S3">
108 /calib/raw/photodiode-sensitivity
109 <HDF5 dataset "data": shape (3,), type "<f8">
110 [ 15900000. 16900000. 16300000.]
111 <HDF5 dataset "units": shape (), type "|S3">
113 /calib/raw/temperature
114 <HDF5 dataset "data": shape (3,), type "<f8">
116 <HDF5 dataset "units": shape (), type "|S1">
118 /calib/raw/thermal-vibration-variance
119 <HDF5 dataset "data": shape (3,), type "<f8">
120 [ 2.20...e-05 2.220...e-05 2.210...e-05]
121 <HDF5 dataset "units": shape (), type "|S3">
124 >>> bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load(
125 ... filename=filename, group='/calib/')
126 >>> (k, k_s) # doctest: +ELLIPSIS
127 (0.0493..., 0.00248...)
129 >>> os.remove(filename)
133 import numpy as _numpy
135 from scipy.constants import Boltzmann as _kB # in J/K
137 from scipy.constants import Bolzmann as _kB # in J/K
138 # Bolzmann -> Boltzmann patch submitted:
139 # http://projects.scipy.org/scipy/ticket/1417
140 # Fixed in scipy commit 4716d91, Apr 2, 2011, during work after v0.9.0rc5.
143 import matplotlib as _matplotlib
144 import matplotlib.pyplot as _matplotlib_pyplot
145 import time as _time # for timestamping lines on plots
146 except (ImportError, RuntimeError), e:
148 _matplotlib_import_error = e
150 from h5config.storage.hdf5 import HDF5_Storage as _HDF5_Storage
151 from h5config.storage.hdf5 import h5_create_group as _h5_create_group
153 from . import LOG as _LOG
154 from . import package_config as _package_config
155 from .bump_analyze import bump_analyze as _bump_analyze
156 from .bump_analyze import bump_load as _bump_load
157 from .bump_analyze import bump_save as _bump_save
158 from .config import CalibrationConfig as _CalibrationConfig
159 from .T_analyze import T_analyze as _temperature_analyze
160 from .T_analyze import T_load as _temperature_load
161 from .T_analyze import T_save as _temperature_save
162 from .vib_analyze import vib_analyze as _vibration_analyze
163 from .vib_analyze import vib_load as _vibration_load
164 from .vib_analyze import vib_save as _vibration_save
167 def calib_analyze(bumps, temperatures, vibrations):
168 """Analyze data from `get_calibration_data()`
170 Inputs (all are arrays of recorded data):
171 bumps measured (V_photodiode / nm_tip) proportionality constant
172 temperatures measured temperature (K)
173 vibrations measured V_photodiode variance in free solution (V**2)
175 k cantilever spring constant (in N/m, or equivalently nN/nm)
176 k_s standard deviation in our estimate of k
180 We're assuming vib is mostly from thermal cantilever vibrations
181 (and then only from vibrations in the single vertical degree of
182 freedom), and not from other noise sources.
184 If the error is large, check the relative errors
185 (`x.std()/x.mean()`)of your input arrays. If one of them is
186 small, don't bother repeating that measurment too often. If one
187 is large, try repeating that measurement more. Remember that you
188 need enough samples to have a valid error estimate in the first
189 place, and that none of this addresses any systematic errors.
191 ps_m = bumps.mean() # ps for photo-sensitivity
193 T_m = temperatures.mean()
194 T_s = temperatures.std()
195 v2_m = vibrations.mean() # average voltage variance
196 v2_s = vibrations.std()
198 # Vphoto / photo_sensitivity = x
199 # k = kB T / <x**2> = kB T photo_sensitivity**2 / Vphoto_var
201 # units, photo_sensitivity = Vphoto/(Zcant in m),
202 # so Vphoto/photo_sensitivity = Zcant in m
203 # so k = J/K * K / m^2 = J / m^2 = N/m
204 k = _kB * T_m * ps_m**2 / v2_m
206 # propogation of errors
210 dk_ps = 2*k/ps_m * ps_s
212 dk_v = -k/v2_m * v2_s
214 k_s = _numpy.sqrt(dk_T**2 + dk_ps**2 + dk_v**2)
216 _LOG.info('variable (units) : '
217 'mean +/- std. dev. (relative error)')
218 _LOG.info('cantilever k (N/m) : %g +/- %g (%g)' % (k, k_s, k_s/k))
219 _LOG.info('photo sensitivity (V/m) : %g +/- %g (%g)'
220 % (ps_m, ps_s, ps_s/ps_m))
221 _LOG.info('T (K) : %g +/- %g (%g)'
222 % (T_m, T_s, T_s/T_m))
223 _LOG.info('vibration variance (V^2) : %g +/- %g (%g)'
224 % (v2_m, v2_s, v2_s/v2_m))
226 if _package_config['matplotlib']:
227 calib_plot(bumps, temperatures, vibrations)
231 def calib_save(filename, group='/', bumps=None, temperatures=None,
232 vibrations=None, calibration_config=None, k=None, k_s=None):
233 with _h5py.File(filename, 'a') as f:
234 cwg = _h5_create_group(f, group)
235 if calibration_config is not None:
236 config_cwg = _h5_create_group(cwg, 'config')
237 calibration_config.save(group=config_cwg)
238 if bumps is not None:
240 del cwg['raw/photodiode-sensitivity/data']
244 del cwg['raw/photodiode-sensitivity/units']
247 cwg['raw/photodiode-sensitivity/data'] = bumps
248 cwg['raw/photodiode-sensitivity/units'] = 'V/m'
249 if temperatures is not None:
251 del cwg['raw/temperature/data']
255 del cwg['raw/temperature/units']
258 cwg['raw/temperature/data'] = temperatures
259 cwg['raw/temperature/units'] = 'K'
260 if vibrations is not None:
262 del cwg['raw/thermal-vibration-variance/data']
266 del cwg['raw/thermal-vibration-variance/units']
269 cwg['raw/thermal-vibration-variance/data'] = vibrations
270 cwg['raw/thermal-vibration-variance/units'] = 'V^2'
273 del cwg['processed/spring-constant/data']
277 del cwg['processed/spring-constant/units']
280 cwg['processed/spring-constant/data'] = k
281 cwg['processed/spring-constant/units'] = 'N/m'
284 del cwg['processed/spring-constant/standard-deviation']
287 cwg['processed/spring-constant/standard-deviation'] = k_s
289 def calib_load(filename, group='/'):
290 assert group.endswith('/')
291 bumps = temperatures = vibrations = k = k_s = None
293 with _h5py.File(filename, 'a') as f:
295 bumps = f[group+'raw/photodiode-sensitivity/data'][...]
299 temperatures = f[group+'raw/temperature/data'][...]
303 vibrations = f[group+'raw/thermal-vibration-variance/data'][...]
307 k = f[group+'processed/spring-constant/data'][...]
311 k_s = f[group+'processed/spring-constant/standard-deviation'][...]
314 calibration_config = _CalibrationConfig(storage=_HDF5_Storage(
315 filename=filename, group=group+'config/'))
316 calibration_config.load()
317 return (bumps, temperatures, vibrations, calibration_config, k, k_s)
319 def calib_plot(bumps, temperatures, vibrations):
321 raise _matplotlib_import_error
322 figure = _matplotlib_pyplot.figure()
324 bump_axes = figure.add_subplot(3, 1, 1)
325 T_axes = figure.add_subplot(3, 1, 2)
326 vib_axes = figure.add_subplot(3, 1, 3)
328 timestamp = _time.strftime('%H%M%S')
329 bump_axes.set_title('cantilever calibration %s' % timestamp)
331 bump_axes.plot(bumps, 'g.-')
332 bump_axes.set_ylabel('photodiode sensitivity (V/m)')
333 T_axes.plot(temperatures, 'r.-')
334 T_axes.set_ylabel('temperature (K)')
335 vib_axes.plot(vibrations, 'b.-')
336 vib_axes.set_ylabel('thermal deflection variance (V^2)')
341 def calib_load_all(filename, group='/'):
342 "Load all data from a `calib()` run."
343 assert group.endswith('/'), group
344 bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load(
345 filename, group+'calibration/')
347 for i in range(calibration_config['num-bumps']):
348 (raw_bump,bump_config,z_channel_config,z_axis_config,
349 deflection_channel_config,processed_bump) = _bump_load(
350 filename=filename, group='%sbump/%d/' % (group, i))
351 bump_details.append({
352 'raw_bump': raw_bump,
353 'bump_config': bump_config,
354 'z_channel_config': z_channel_config,
355 'z_axis_config': z_axis_config,
356 'deflection_channel_config': deflection_channel_config,
357 'processed_bump': processed_bump,
359 temperature_details = []
360 for i in range(calibration_config['num-temperatures']):
361 (raw_temperature,temperature_config,processed_temperature
362 ) = _temperature_load(
363 filename=filename, group='%stemperature/%d/' % (group, i))
364 temperature_details.append({
365 'raw_temperature': raw_temperature,
366 'temperature_config': temperature_config,
367 'processed_temperature': processed_temperature,
369 vibration_details = []
370 for i in range(calibration_config['num-vibrations']):
371 (raw_vibration,vibration_config,deflection_channel_config,
372 processed_vibration) = _vibration_load(
373 filename=filename, group='%svibration/%d/' % (group, i))
374 vibration_details.append({
375 'raw_vibration': raw_vibration,
376 'vibration_config': vibration_config,
377 'deflection_channel_config': deflection_channel_config,
378 'processed_vibration': processed_vibration,
382 'bump_details': bump_details,
383 'temperatures': temperatures,
384 'temperature_details': temperature_details,
385 'vibrations': vibrations,
386 'vibration_details': vibration_details,
387 'calibration_config': calibration_config,
392 def calib_analyze_all(filename, group='/', maximum_relative_error=1e-5,
394 "(Re)analyze all data from a `calib()` run."
395 assert group.endswith('/'), group
396 bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load(
397 filename, group+'calibration/')
398 changed_bump = changed_temperature = changed_vibration = False
399 for i in range(calibration_config['num-bumps']):
400 bump_group = '%sbump/%d/' % (group, i)
401 (raw_bump,bump_config,z_channel_config,z_axis_config,
402 deflection_channel_config,processed_bump) = _bump_load(
403 filename=filename, group=bump_group)
404 sensitivity = _bump_analyze(
405 data=raw_bump, bump_config=bump_config,
406 z_channel_config=z_channel_config,
407 z_axis_config=z_axis_config,
408 deflection_channel_config=deflection_channel_config)
409 bumps[i] = sensitivity
410 rel_error = abs(sensitivity - processed_bump)/processed_bump
411 if rel_error > maximum_relative_error:
412 _LOG.warn(("new analysis doesn't match for bump %d: %g -> %g "
413 "(difference: %g, relative error: %g)")
414 % (i, processed_bump, sensitivity,
415 sensitivity-processed_bump, rel_error))
418 _bump_save(filename, bump_group, processed_bump=sensitivity)
419 for i in range(calibration_config['num-temperatures']):
420 temperature_group = '%stemperature/%d/' % (group, i)
421 (raw_temperature,temperature_config,processed_temperature
422 ) = _temperature_load(
423 filename=filename, group=temperature_group)
424 temperature = _temperature_analyze(
425 raw_temperature, temperature_config)
426 temperatures[i] = temperature
427 rel_error = abs(temperature - processed_temperature
428 )/processed_temperature
429 if rel_error > maximum_relative_error:
430 _LOG.warn(("new analysis doesn't match for temperature %d: "
431 "%g -> %g (difference: %g, relative error: %g)")
432 % (i, processed_temperature, temperature,
433 temperature-processed_temperature, rel_error))
434 changed_temperature = True
437 filename, temperature_group,
438 processed_T=temperature)
439 for i in range(calibration_config['num-vibrations']):
440 vibration_group = '%svibration/%d/' % (group, i)
441 (raw_vibration,vibration_config,deflection_channel_config,
442 processed_vibration) = _vibration_load(
443 filename=filename, group=vibration_group)
444 variance = _vibration_analyze(
445 deflection=raw_vibration, vibration_config=vibration_config,
446 deflection_channel_config=deflection_channel_config)
447 vibrations[i] = variance
448 rel_error = abs(variance - processed_vibration)/processed_vibration
449 if rel_error > maximum_relative_error:
450 _LOG.warn(("new analysis doesn't match for vibration %d: %g -> %g "
451 "(difference: %g, relative error: %g)")
452 % (i, processed_vibration, variance,
453 variance-processed_vibration, rel_error))
454 changed_vibration = True
457 filename, vibration_group, processed_vibration=variance)
459 calib_group = '%scalibration/' % group
461 if changed_bump and not dry_run:
462 calib_save(filename, calib_group, bumps=bumps)
463 if changed_temperature and not dry_run:
464 calib_save(filename, calib_group, temperatures=temperatures)
465 if changed_vibration and not dry_run:
466 calib_save(filename, calib_group, vibrations=vibrations)
468 new_k,new_k_s = calib_analyze(
469 bumps=bumps, temperatures=temperatures, vibrations=vibrations)
470 rel_error = abs(new_k-k)/k
471 if rel_error > maximum_relative_error:
472 _LOG.warn(("new analysis doesn't match for k: %g -> %g "
473 "(difference: %g, relative error: %g)")
474 % (k, new_k, new_k-k, rel_error))
476 calib_save(filename, calib_group, k=new_k)
477 rel_error = abs(new_k_s-k_s)/k_s
478 if rel_error > maximum_relative_error:
479 _LOG.warn(("new analysis doesn't match for k_s: %g -> %g "
480 "(difference: %g, relative error: %g)")
481 % (k_s, new_k_s, new_k_s-k_s, rel_error))
483 calib_save(filename, calib_group, k_s=new_k_s)
484 return (new_k, new_k_s)
486 def calib_plot_all(bumps, bump_details, temperatures, temperature_details,
487 vibrations, vibration_details, calibration_config, k, k_s,
488 maximum_relative_error=1e-5):
489 calib_plot(bumps, temperatures, vibrations)
490 for i,bump in enumerate(bump_details):
491 sensitivity = _bump_analyze(
492 data=bump['raw_bump'], bump_config=bump['bump_config'],
493 z_channel_config=bump['z_channel_config'],
494 z_axis_config=bump['z_axis_config'],
495 deflection_channel_config=bump['deflection_channel_config'],
497 rel_error = abs(sensitivity - bump['processed_bump']
498 )/bump['processed_bump']
499 if rel_error > maximum_relative_error:
500 _LOG.warn(("new analysis doesn't match for bump %d: %g != %g "
501 "(difference: %g, relative error: %g)")
502 % (i, sensitivity, bump['processed_bump'],
503 sensitivity-bump['processed_bump'], rel_error))
504 # nothing interesting to plot for temperatures...
505 for i,vibration in enumerate(vibration_details):
506 variance = _vibration_analyze(
507 deflection=vibration['raw_vibration'],
508 vibration_config=vibration['vibration_config'],
509 deflection_channel_config=vibration['deflection_channel_config'],
511 rel_error = abs(variance - vibration['processed_vibration']
512 )/vibration['processed_vibration']
513 if rel_error > maximum_relative_error:
514 _LOG.warn(("new analysis doesn't match for vibration %d: %g != %g "
515 "(difference: %g, relative error: %g)")
516 % (i, variance, vibration['processed_vibration'],
517 variance-vibration['processed_vibration'], rel_error))