1 # calibcant - tools for thermally calibrating AFM cantilevers
3 # Copyright (C) 2008-2012 W. Trevor King <wking@drexel.edu>
5 # This file is part of calibcant.
7 # calibcant is free software: you can redistribute it and/or modify it under
8 # the terms of the GNU General Public License as published by the Free Software
9 # Foundation, either version 3 of the License, or (at your option) any later
12 # calibcant is distributed in the hope that it will be useful, but WITHOUT ANY
13 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
14 # A PARTICULAR PURPOSE. See the GNU General Public License for more details.
16 # You should have received a copy of the GNU General Public License along with
17 # calibcant. If not, see <http://www.gnu.org/licenses/>.
19 """Calculate `k` from arrays of bumps, temperatures, and vibrations.
21 Separate the more general `analyze()` from the other calibration
22 functions in calibcant.
24 The relevent physical quantities are :
25 Vzp_out Output z-piezo voltage (what we generate)
26 Vzp Applied z-piezo voltage (after external ZPGAIN)
27 Zp The z-piezo position
28 Zcant The cantilever vertical deflection
29 Vphoto The photodiode vertical deflection voltage (what we measure)
30 Fcant The force on the cantilever
31 T The temperature of the cantilever and surrounding solution
32 (another thing we measure)
33 k_b Boltzmann's constant
35 Which are related by the parameters:
37 zp_sensitivity Zp / Vzp
38 photo_sensitivity Vphoto / Zcant
43 >>> from .config import CalibrateConfig
45 >>> config = CalibrateConfig()
46 >>> bumps = numpy.array((15.9e6, 16.9e6, 16.3e6))
47 >>> temperatures = numpy.array((295, 295.2, 294.8))
48 >>> vibrations = numpy.array((2.20e-5, 2.22e-5, 2.21e-5))
50 >>> k,k_s = analyze(bumps=bumps, temperatures=temperatures,
51 ... vibrations=vibrations)
52 >>> (k, k_s) # doctest: +ELLIPSIS
53 (0.0493..., 0.00248...)
55 Most of the error in this example comes from uncertainty in the
56 photodiode sensitivity (bumps).
58 >>> k_s/k # doctest: +ELLIPSIS
60 >>> bumps.std()/bumps.mean() # doctest: +ELLIPSIS
62 >>> temperatures.std()/temperatures.mean() # doctest: +ELLIPSIS
64 >>> vibrations.std()/vibrations.mean() # doctest: +ELLIPSIS
69 import numpy as _numpy
71 from scipy.constants import Boltzmann as _kB # in J/K
73 from scipy.constants import Bolzmann as _kB # in J/K
74 # Bolzmann -> Boltzmann patch submitted:
75 # http://projects.scipy.org/scipy/ticket/1417
76 # Fixed in scipy commit 4716d91, Apr 2, 2011, during work after v0.9.0rc5.
79 import matplotlib as _matplotlib
80 import matplotlib.pyplot as _matplotlib_pyplot
81 import time as _time # for timestamping lines on plots
82 except (ImportError, RuntimeError), e:
84 _matplotlib_import_error = e
86 from h5config.storage.hdf5 import h5_create_group as _h5_create_group
87 from pypiezo.base import get_axis_name as _get_axis_name
89 from . import LOG as _LOG
90 from . import package_config as _package_config
92 from .bump_analyze import analyze as _bump_analyze
93 from .bump_analyze import save as _bump_save
94 from .temperature_analyze import analyze as _temperature_analyze
95 from .temperature_analyze import save as _temperature_save
96 from .vibration_analyze import analyze as _vibration_analyze
97 from .vibration_analyze import save as _vibration_save
100 def analyze(bumps, temperatures, vibrations):
101 """Analyze data from `get_calibration_data()`
103 Inputs (all are arrays of recorded data):
104 bumps measured (V_photodiode / nm_tip) proportionality constant
105 temperatures measured temperature (K)
106 vibrations measured V_photodiode variance in free solution (V**2)
108 k cantilever spring constant (in N/m, or equivalently nN/nm)
109 k_s standard deviation in our estimate of k
113 We're assuming vib is mostly from thermal cantilever vibrations
114 (and then only from vibrations in the single vertical degree of
115 freedom), and not from other noise sources.
117 If the error is large, check the relative errors
118 (`x.std()/x.mean()`)of your input arrays. If one of them is
119 small, don't bother repeating that measurment too often. If one
120 is large, try repeating that measurement more. Remember that you
121 need enough samples to have a valid error estimate in the first
122 place, and that none of this addresses any systematic errors.
124 ps_m = bumps.mean() # ps for photo-sensitivity
126 T_m = temperatures.mean()
127 T_s = temperatures.std()
128 v2_m = vibrations.mean() # average voltage variance
129 v2_s = vibrations.std()
131 # Vphoto / photo_sensitivity = x
132 # k = kB T / <x**2> = kB T photo_sensitivity**2 / Vphoto_var
134 # units, photo_sensitivity = Vphoto/(Zcant in m),
135 # so Vphoto/photo_sensitivity = Zcant in m
136 # so k = J/K * K / m^2 = J / m^2 = N/m
137 k = _kB * T_m * ps_m**2 / v2_m
139 # propogation of errors
143 dk_ps = 2*k/ps_m * ps_s
145 dk_v = -k/v2_m * v2_s
147 k_s = _numpy.sqrt(dk_T**2 + dk_ps**2 + dk_v**2)
149 _LOG.info('variable (units) : '
150 'mean +/- std. dev. (relative error)')
151 _LOG.info('cantilever k (N/m) : %g +/- %g (%g)' % (k, k_s, k_s/k))
152 _LOG.info('photo sensitivity (V/m) : %g +/- %g (%g)'
153 % (ps_m, ps_s, ps_s/ps_m))
154 _LOG.info('T (K) : %g +/- %g (%g)'
155 % (T_m, T_s, T_s/T_m))
156 _LOG.info('vibration variance (V^2) : %g +/- %g (%g)'
157 % (v2_m, v2_s, v2_s/v2_m))
159 if _package_config['matplotlib']:
160 plot(bumps, temperatures, vibrations)
165 def plot(bumps, temperatures, vibrations):
167 raise _matplotlib_import_error
168 figure = _matplotlib_pyplot.figure()
170 bump_axes = figure.add_subplot(3, 1, 1)
171 T_axes = figure.add_subplot(3, 1, 2)
172 vib_axes = figure.add_subplot(3, 1, 3)
174 timestamp = _time.strftime('%H%M%S')
175 bump_axes.set_title('cantilever calibration %s' % timestamp)
177 bump_axes.autoscale(tight=True)
178 bump_axes.plot(bumps, 'g.-')
179 bump_axes.set_ylabel('photodiode sensitivity (V/m)')
180 T_axes.autoscale(tight=True)
181 T_axes.plot(temperatures, 'r.-')
182 T_axes.set_ylabel('temperature (K)')
183 vib_axes.autoscale(tight=True)
184 vib_axes.plot(vibrations, 'b.-')
185 vib_axes.set_ylabel('thermal deflection variance (V^2)')
187 if hasattr(figure, 'show'):
190 _plot = plot # alternative name for use inside analyze_all()
193 def analyze_all(config, data, raw_data, maximum_relative_error=1e-5,
194 filename=None, group=None, plot=False, dry_run=False):
195 "(Re)analyze (and possibly plot) all data from a `calib()` run."
196 if not data.get('bump', None):
197 data['bump'] = _numpy.zeros((config['num-bumps'],), dtype=float)
198 if not data.get('temperature', None):
199 data['temperature'] = _numpy.zeros(
200 (config['num-temperatures'],), dtype=float)
201 if not data.get('vibrations', None):
202 data['vibration'] = _numpy.zeros(
203 (config['num-vibrations'],), dtype=float)
204 axis_config = config['afm']['piezo'].select_config(
206 attribute_value=config['afm']['main-axis'],
207 get_attribute=_get_axis_name)
208 input_config = config['afm']['piezo'].select_config(
209 setting_name='inputs', attribute_value='deflection')
210 bumps_changed = temperatures_changed = vibrations_changed = False
211 if not isinstance(group, _h5py.Group) and not dry_run:
212 f = _h5py.File(filename, mode='a')
213 group = _h5_create_group(f, group)
217 for i,bump in enumerate(raw_data['bump']):
218 data['bump'][i],changed = check_bump(
219 index=i, bump=bump, z_axis_config=axis_config,
220 deflection_channel_config=input_config, plot=plot,
221 maximum_relative_error=maximum_relative_error)
222 if changed and not dry_run:
224 bump_group = _h5_create_group(group, 'bump/{}'.format(i))
225 _bump_save(group=bump_group, processed=data['bump'][i])
226 for i,temperature in enumerate(raw_data['temperature']):
227 data['temperature'][i],changed = check_temperature(
228 index=i, temperature=temperature,
229 maximum_relative_error=maximum_relative_error)
230 if changed and not dry_run:
231 temperatures_changed = True
232 temperature_group = _h5_create_group(
233 group, 'temperature/{}'.format(i))
235 group=temerature_group, processed=data['temperature'][i])
236 for i,vibration in enumerate(raw_data['vibration']):
237 data['vibration'][i],changed = check_vibration(
238 index=i, vibration=vibration,
239 deflection_channel_config=input_config, plot=plot,
240 maximum_relative_error=maximum_relative_error)
241 if changed and not dry_run:
242 vibrations_changed = True
243 vibration_group = _h5_create_group(
244 group, 'vibration/{}'.format(i))
246 group=vibration_group, processed=data['vibration'])
247 k,k_s,changed = check_calibration(
248 k=data['processed']['spring_constant'],
249 k_s=data['processed']['spring_constant_deviation'],
251 temperatures=data['temperature'], vibrations=data['vibration'],
252 maximum_relative_error=maximum_relative_error)
253 if (changed or bumps_changed or temperatures_changed or
254 vibrations_changed) and not dry_run:
255 calibration_group = _h5_create_group(group, 'calibration')
257 calib_save(group=calibration_group, bump=data['bump'])
258 if temperatures_changed:
260 group=calibration_group, temperature=data['temperature'])
261 if vibrations_changed:
263 group=calibration_group, vibration=data['vibration'])
265 calib_save(group=calibration_group, k=k, k_s=k_s)
270 _plot(bumps=data['raw']['bump'],
271 temperatures=data['raw']['temperature'],
272 vibrations=data['raw']['vibration'])
275 def check_bump(index, bump, maximum_relative_error, **kwargs):
277 sensitivity = _bump_analyze(
278 config=bump['config']['bump'], data=bump['raw'], **kwargs)
279 if bump.get('processed', None) is None:
281 _LOG.warn('new analysis for bump {}: {}'.format(index, sensitivity))
283 rel_error = abs(sensitivity - bump['processed'])/bump['processed']
284 if rel_error > maximum_relative_error:
286 _LOG.warn(("new analysis doesn't match for bump {}: {} -> {} "
287 "(difference: {}, relative error: {})").format(
288 index, bump['processed'], sensitivity,
289 sensitivity-bump['processed'], rel_error))
290 return (sensitivity, changed)
292 def check_temperature(index, temperature, maximum_relative_error, **kwargs):
294 temp = _temperature_analyze(
295 config=temperature['config']['temperature'],
296 temperature=temperature['raw'], **kwargs)
297 if temperature.get('processed', None) is None:
299 _LOG.warn('new analysis for temperature {}: {}'.format(index, temp))
301 rel_error = abs(temp - temperature['processed']
302 )/temperature['processed']
303 if rel_error > maximum_relative_error:
305 _LOG.warn(("new analysis doesn't match for temperature "
306 "{} -> {} (difference: {}, relative error: {})"
308 index, temperature['processed'], temp,
309 temp-temperature['processed'], rel_error))
310 return (temp, changed)
312 def check_vibration(index, vibration, maximum_relative_error, **kwargs):
314 variance = _vibration_analyze(
315 config=vibration['config']['vibration'],
316 deflection=vibration['raw'], **kwargs)
317 if vibration.get('processed', None) is None:
319 _LOG.warn('new analysis for temperature {}: {}'.format(
322 rel_error = abs(variance-vibration['processed'])/vibration['processed']
323 if rel_error > maximum_relative_error:
324 _LOG.warn(("new analysis doesn't match for vibration {}: {} != {} "
325 "(difference: {}, relative error: {})").format(
326 index, variance, vibration['processed'],
327 variance-vibration['processed'], rel_error))
328 return (variance, changed)
330 def check_calibration(k, k_s, maximum_relative_error, **kwargs):
332 new_k,new_k_s = analyze(**kwargs)
335 _LOG.warn('new analysis for the spring constant: {}'.format(new_k))
337 rel_error = abs(new_k-k)/k
338 if rel_error > maximum_relative_error:
339 _LOG.warn(("new analysis doesn't match for the spring constant: "
340 "{} != {} (difference: {}, relative error: {})").format(
341 new_k, k, new_k-k, rel_error))
344 _LOG.warn('new analysis for the spring constant deviation: {}'.format(
347 rel_error = abs(new_k-k)/k
348 if rel_error > maximum_relative_error:
350 ("new analysis doesn't match for the spring constant deviation"
351 ": {} != {} (difference: {}, relative error: {})").format(
352 new_k_s, k_s, new_k_s-k_s, rel_error))
353 return (new_k, new_k_s, changed)