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.hdf5 import pprint_HDF5
48 >>> from .config import HDF5_CalibrationConfig
50 >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
53 >>> calibration_config = HDF5_CalibrationConfig(filename, '/calib/config/')
54 >>> bumps = numpy.array((15.9e6, 16.9e6, 16.3e6))
55 >>> temperatures = numpy.array((295, 295.2, 294.8))
56 >>> vibrations = numpy.array((2.20e-5, 2.22e-5, 2.21e-5))
58 >>> k,k_s = calib_analyze(bumps=bumps, temperatures=temperatures,
59 ... vibrations=vibrations)
60 >>> (k, k_s) # doctest: +ELLIPSIS
61 (0.0493..., 0.00248...)
63 Most of the error in this example comes from uncertainty in the
64 photodiode sensitivity (bumps).
66 >>> k_s/k # doctest: +ELLIPSIS
68 >>> bumps.std()/bumps.mean() # doctest: +ELLIPSIS
70 >>> temperatures.std()/temperatures.mean() # doctest: +ELLIPSIS
72 >>> vibrations.std()/vibrations.mean() # doctest: +ELLIPSIS
75 >>> calib_save(filename=filename, group='/calib/',
76 ... bumps=bumps, temperatures=temperatures, vibrations=vibrations,
77 ... calibration_config=calibration_config, k=k, k_s=k_s)
78 >>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
82 <HDF5 dataset "num-bumps": shape (), type "|S2">
84 <HDF5 dataset "num-temperatures": shape (), type "|S2">
86 <HDF5 dataset "num-vibrations": shape (), type "|S2">
88 <HDF5 dataset "temperature-sleep": shape (), type "|S1">
90 <HDF5 dataset "vibration-spacing": shape (), type "|S5">
93 /calib/processed/spring-constant
94 <HDF5 dataset "data": shape (), type "<f8">
96 <HDF5 dataset "standard-deviation": shape (), type "<f8">
98 <HDF5 dataset "units": shape (), type "|S3">
101 /calib/raw/photodiode-sensitivity
102 <HDF5 dataset "data": shape (3,), type "<f8">
103 [ 15900000. 16900000. 16300000.]
104 <HDF5 dataset "units": shape (), type "|S3">
106 /calib/raw/temperature
107 <HDF5 dataset "data": shape (3,), type "<f8">
109 <HDF5 dataset "units": shape (), type "|S1">
111 /calib/raw/thermal-vibration-variance
112 <HDF5 dataset "data": shape (3,), type "<f8">
113 [ 2.20...e-05 2.220...e-05 2.210...e-05]
114 <HDF5 dataset "units": shape (), type "|S3">
117 >>> bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load(
118 ... filename=filename, group='/calib/')
119 >>> (k, k_s) # doctest: +ELLIPSIS
120 (0.0493..., 0.00248...)
122 >>> os.remove(filename)
126 import numpy as _numpy
128 from scipy.constants import Boltzmann as _kB # in J/K
130 from scipy.constants import Bolzmann as _kB # in J/K
131 # Bolzmann -> Boltzmann patch submitted:
132 # http://projects.scipy.org/scipy/ticket/1417
133 # Fixed in scipy commit 4716d91, Apr 2, 2011, during work after v0.9.0rc5.
136 import matplotlib as _matplotlib
137 import matplotlib.pyplot as _matplotlib_pyplot
138 import time as _time # for timestamping lines on plots
139 except (ImportError, RuntimeError), e:
141 _matplotlib_import_error = e
143 from h5config.hdf5 import h5_create_group as _h5_create_group
145 from . import LOG as _LOG
146 from . import package_config as _package_config
147 from .bump_analyze import bump_analyze as _bump_analyze
148 from .bump_analyze import bump_load as _bump_load
149 from .bump_analyze import bump_save as _bump_save
150 from .config import HDF5_CalibrationConfig as _HDF5_CalibrationConfig
151 from .T_analyze import T_analyze as _temperature_analyze
152 from .T_analyze import T_load as _temperature_load
153 from .T_analyze import T_save as _temperature_save
154 from .vib_analyze import vib_analyze as _vibration_analyze
155 from .vib_analyze import vib_load as _vibration_load
156 from .vib_analyze import vib_save as _vibration_save
159 def calib_analyze(bumps, temperatures, vibrations):
160 """Analyze data from `get_calibration_data()`
162 Inputs (all are arrays of recorded data):
163 bumps measured (V_photodiode / nm_tip) proportionality constant
164 temperatures measured temperature (K)
165 vibrations measured V_photodiode variance in free solution (V**2)
167 k cantilever spring constant (in N/m, or equivalently nN/nm)
168 k_s standard deviation in our estimate of k
172 We're assuming vib is mostly from thermal cantilever vibrations
173 (and then only from vibrations in the single vertical degree of
174 freedom), and not from other noise sources.
176 If the error is large, check the relative errors
177 (`x.std()/x.mean()`)of your input arrays. If one of them is
178 small, don't bother repeating that measurment too often. If one
179 is large, try repeating that measurement more. Remember that you
180 need enough samples to have a valid error estimate in the first
181 place, and that none of this addresses any systematic errors.
183 ps_m = bumps.mean() # ps for photo-sensitivity
185 T_m = temperatures.mean()
186 T_s = temperatures.std()
187 v2_m = vibrations.mean() # average voltage variance
188 v2_s = vibrations.std()
190 # Vphoto / photo_sensitivity = x
191 # k = kB T / <x**2> = kB T photo_sensitivity**2 / Vphoto_var
193 # units, photo_sensitivity = Vphoto/(Zcant in m),
194 # so Vphoto/photo_sensitivity = Zcant in m
195 # so k = J/K * K / m^2 = J / m^2 = N/m
196 k = _kB * T_m * ps_m**2 / v2_m
198 # propogation of errors
202 dk_ps = 2*k/ps_m * ps_s
204 dk_v = -k/v2_m * v2_s
206 k_s = _numpy.sqrt(dk_T**2 + dk_ps**2 + dk_v**2)
208 _LOG.info('variable (units) : '
209 'mean +/- std. dev. (relative error)')
210 _LOG.info('cantilever k (N/m) : %g +/- %g (%g)' % (k, k_s, k_s/k))
211 _LOG.info('photo sensitivity (V/m) : %g +/- %g (%g)'
212 % (ps_m, ps_s, ps_s/ps_m))
213 _LOG.info('T (K) : %g +/- %g (%g)'
214 % (T_m, T_s, T_s/T_m))
215 _LOG.info('vibration variance (V^2) : %g +/- %g (%g)'
216 % (v2_m, v2_s, v2_s/v2_m))
218 if _package_config['matplotlib']:
219 calib_plot(bumps, temperatures, vibrations)
223 def calib_save(filename, group='/', bumps=None, temperatures=None,
224 vibrations=None, calibration_config=None, k=None, k_s=None):
225 with _h5py.File(filename, 'a') as f:
226 cwg = _h5_create_group(f, group)
227 if calibration_config is not None:
228 config_cwg = _h5_create_group(cwg, 'config')
229 calibration_config.save(group=config_cwg)
230 if bumps is not None:
232 del cwg['raw/photodiode-sensitivity/data']
236 del cwg['raw/photodiode-sensitivity/units']
239 cwg['raw/photodiode-sensitivity/data'] = bumps
240 cwg['raw/photodiode-sensitivity/units'] = 'V/m'
241 if temperatures is not None:
243 del cwg['raw/temperature/data']
247 del cwg['raw/temperature/units']
250 cwg['raw/temperature/data'] = temperatures
251 cwg['raw/temperature/units'] = 'K'
252 if vibrations is not None:
254 del cwg['raw/thermal-vibration-variance/data']
258 del cwg['raw/thermal-vibration-variance/units']
261 cwg['raw/thermal-vibration-variance/data'] = vibrations
262 cwg['raw/thermal-vibration-variance/units'] = 'V^2'
265 del cwg['processed/spring-constant/data']
269 del cwg['processed/spring-constant/units']
272 cwg['processed/spring-constant/data'] = k
273 cwg['processed/spring-constant/units'] = 'N/m'
276 del cwg['processed/spring-constant/standard-deviation']
279 cwg['processed/spring-constant/standard-deviation'] = k_s
281 def calib_load(filename, group='/'):
282 assert group.endswith('/')
283 bumps = temperatures = vibrations = k = k_s = None
285 with _h5py.File(filename, 'a') as f:
287 bumps = f[group+'raw/photodiode-sensitivity/data'][...]
291 temperatures = f[group+'raw/temperature/data'][...]
295 vibrations = f[group+'raw/thermal-vibration-variance/data'][...]
299 k = f[group+'processed/spring-constant/data'][...]
303 k_s = f[group+'processed/spring-constant/standard-deviation'][...]
306 calibration_config = _HDF5_CalibrationConfig(
307 filename=filename, group=group+'config/')
308 calibration_config.load()
309 return (bumps, temperatures, vibrations, calibration_config, k, k_s)
311 def calib_plot(bumps, temperatures, vibrations):
313 raise _matplotlib_import_error
314 figure = _matplotlib_pyplot.figure()
316 bump_axes = figure.add_subplot(3, 1, 1)
317 T_axes = figure.add_subplot(3, 1, 2)
318 vib_axes = figure.add_subplot(3, 1, 3)
320 timestamp = _time.strftime('%H%M%S')
321 bump_axes.set_title('cantilever calibration %s' % timestamp)
323 bump_axes.plot(bumps, 'g.-')
324 bump_axes.set_ylabel('photodiode sensitivity (V/m)')
325 T_axes.plot(temperatures, 'r.-')
326 T_axes.set_ylabel('temperature (K)')
327 vib_axes.plot(vibrations, 'b.-')
328 vib_axes.set_ylabel('thermal deflection variance (V^2)')
333 def calib_load_all(filename, group='/'):
334 "Load all data from a `calib()` run."
335 assert group.endswith('/'), group
336 bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load(
337 filename, group+'calibration/')
339 for i in range(calibration_config['num-bumps']):
340 (raw_bump,bump_config,z_channel_config,z_axis_config,
341 deflection_channel_config,processed_bump) = _bump_load(
342 filename=filename, group='%sbump/%d/' % (group, i))
343 bump_details.append({
344 'raw_bump': raw_bump,
345 'bump_config': bump_config,
346 'z_channel_config': z_channel_config,
347 'z_axis_config': z_axis_config,
348 'deflection_channel_config': deflection_channel_config,
349 'processed_bump': processed_bump,
351 temperature_details = []
352 for i in range(calibration_config['num-temperatures']):
353 (raw_temperature,temperature_config,processed_temperature
354 ) = _temperature_load(
355 filename=filename, group='%stemperature/%d/' % (group, i))
356 temperature_details.append({
357 'raw_temperature': raw_temperature,
358 'temperature_config': temperature_config,
359 'processed_temperature': processed_temperature,
361 vibration_details = []
362 for i in range(calibration_config['num-vibrations']):
363 (raw_vibration,vibration_config,deflection_channel_config,
364 processed_vibration) = _vibration_load(
365 filename=filename, group='%svibration/%d/' % (group, i))
366 vibration_details.append({
367 'raw_vibration': raw_vibration,
368 'vibration_config': vibration_config,
369 'deflection_channel_config': deflection_channel_config,
370 'processed_vibration': processed_vibration,
374 'bump_details': bump_details,
375 'temperatures': temperatures,
376 'temperature_details': temperature_details,
377 'vibrations': vibrations,
378 'vibration_details': vibration_details,
379 'calibration_config': calibration_config,
384 def calib_analyze_all(filename, group='/', maximum_relative_error=1e-5,
386 "(Re)analyze all data from a `calib()` run."
387 assert group.endswith('/'), group
388 bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load(
389 filename, group+'calibration/')
390 changed_bump = changed_temperature = changed_vibration = False
391 for i in range(calibration_config['num-bumps']):
392 bump_group = '%sbump/%d/' % (group, i)
393 (raw_bump,bump_config,z_channel_config,z_axis_config,
394 deflection_channel_config,processed_bump) = _bump_load(
395 filename=filename, group=bump_group)
396 sensitivity = _bump_analyze(
397 data=raw_bump, bump_config=bump_config,
398 z_channel_config=z_channel_config,
399 z_axis_config=z_axis_config,
400 deflection_channel_config=deflection_channel_config)
401 bumps[i] = sensitivity
402 rel_error = abs(sensitivity - processed_bump)/processed_bump
403 if rel_error > maximum_relative_error:
404 _LOG.warn(("new analysis doesn't match for bump %d: %g -> %g "
405 "(difference: %g, relative error: %g)")
406 % (i, processed_bump, sensitivity,
407 sensitivity-processed_bump, rel_error))
410 _bump_save(filename, bump_group, processed_bump=sensitivity)
411 for i in range(calibration_config['num-temperatures']):
412 temperature_group = '%stemperature/%d/' % (group, i)
413 (raw_temperature,temperature_config,processed_temperature
414 ) = _temperature_load(
415 filename=filename, group=temperature_group)
416 temperature = _temperature_analyze(
417 raw_temperature, temperature_config)
418 temperatures[i] = temperature
419 rel_error = abs(temperature - processed_temperature
420 )/processed_temperature
421 if rel_error > maximum_relative_error:
422 _LOG.warn(("new analysis doesn't match for temperature %d: "
423 "%g -> %g (difference: %g, relative error: %g)")
424 % (i, processed_temperature, temperature,
425 temperature-processed_temperature, rel_error))
426 changed_temperature = True
429 filename, temperature_group,
430 processed_T=temperature)
431 for i in range(calibration_config['num-vibrations']):
432 vibration_group = '%svibration/%d/' % (group, i)
433 (raw_vibration,vibration_config,deflection_channel_config,
434 processed_vibration) = _vibration_load(
435 filename=filename, group=vibration_group)
436 variance = _vibration_analyze(
437 deflection=raw_vibration, vibration_config=vibration_config,
438 deflection_channel_config=deflection_channel_config)
439 vibrations[i] = variance
440 rel_error = abs(variance - processed_vibration)/processed_vibration
441 if rel_error > maximum_relative_error:
442 _LOG.warn(("new analysis doesn't match for vibration %d: %g -> %g "
443 "(difference: %g, relative error: %g)")
444 % (i, processed_vibration, variance,
445 variance-processed_vibration, rel_error))
446 changed_vibration = True
449 filename, vibration_group, processed_vibration=variance)
451 calib_group = '%scalibration/' % group
453 if changed_bump and not dry_run:
454 calib_save(filename, calib_group, bumps=bumps)
455 if changed_temperature and not dry_run:
456 calib_save(filename, calib_group, temperatures=temperatures)
457 if changed_vibration and not dry_run:
458 calib_save(filename, calib_group, vibrations=vibrations)
460 new_k,new_k_s = calib_analyze(
461 bumps=bumps, temperatures=temperatures, vibrations=vibrations)
462 rel_error = abs(new_k-k)/k
463 if rel_error > maximum_relative_error:
464 _LOG.warn(("new analysis doesn't match for k: %g -> %g "
465 "(difference: %g, relative error: %g)")
466 % (k, new_k, new_k-k, rel_error))
468 calib_save(filename, calib_group, k=new_k)
469 rel_error = abs(new_k_s-k_s)/k_s
470 if rel_error > maximum_relative_error:
471 _LOG.warn(("new analysis doesn't match for k_s: %g -> %g "
472 "(difference: %g, relative error: %g)")
473 % (k_s, new_k_s, new_k_s-k_s, rel_error))
475 calib_save(filename, calib_group, k_s=new_k_s)
476 return (new_k, new_k_s)
478 def calib_plot_all(bumps, bump_details, temperatures, temperature_details,
479 vibrations, vibration_details, calibration_config, k, k_s,
480 maximum_relative_error=1e-5):
481 calib_plot(bumps, temperatures, vibrations)
482 for i,bump in enumerate(bump_details):
483 sensitivity = _bump_analyze(
484 data=bump['raw_bump'], bump_config=bump['bump_config'],
485 z_channel_config=bump['z_channel_config'],
486 z_axis_config=bump['z_axis_config'],
487 deflection_channel_config=bump['deflection_channel_config'],
489 rel_error = abs(sensitivity - bump['processed_bump']
490 )/bump['processed_bump']
491 if rel_error > maximum_relative_error:
492 _LOG.warn(("new analysis doesn't match for bump %d: %g != %g "
493 "(difference: %g, relative error: %g)")
494 % (i, sensitivity, bump['processed_bump'],
495 sensitivity-bump['processed_bump'], rel_error))
496 # nothing interesting to plot for temperatures...
497 for i,vibration in enumerate(vibration_details):
498 variance = _vibration_analyze(
499 deflection=vibration['raw_vibration'],
500 vibration_config=vibration['vibration_config'],
501 deflection_channel_config=vibration['deflection_channel_config'],
503 rel_error = abs(variance - vibration['processed_vibration']
504 )/vibration['processed_vibration']
505 if rel_error > maximum_relative_error:
506 _LOG.warn(("new analysis doesn't match for vibration %d: %g != %g "
507 "(difference: %g, relative error: %g)")
508 % (i, variance, vibration['processed_vibration'],
509 variance-vibration['processed_vibration'], rel_error))