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.plot(bumps, 'g.-')
178 bump_axes.set_ylabel('photodiode sensitivity (V/m)')
179 T_axes.plot(temperatures, 'r.-')
180 T_axes.set_ylabel('temperature (K)')
181 vib_axes.plot(vibrations, 'b.-')
182 vib_axes.set_ylabel('thermal deflection variance (V^2)')
184 if hasattr(figure, 'show'):
187 _plot = plot # alternative name for use inside analyze_all()
190 def analyze_all(config, data, raw_data, maximum_relative_error=1e-5,
191 filename=None, group=None, plot=False, dry_run=False):
192 "(Re)analyze (and possibly plot) all data from a `calib()` run."
193 if not data.get('bump', None):
194 data['bump'] = _numpy.zeros((config['num-bumps'],), dtype=float)
195 if not data.get('temperature', None):
196 data['temperature'] = _numpy.zeros(
197 (config['num-temperatures'],), dtype=float)
198 if not data.get('vibrations', None):
199 data['vibration'] = _numpy.zeros(
200 (config['num-vibrations'],), dtype=float)
201 axis_config = config['afm']['piezo'].select_config(
203 attribute_value=config['afm']['main-axis'],
204 get_attribute=_get_axis_name)
205 input_config = config['afm']['piezo'].select_config(
206 setting_name='inputs', attribute_value='deflection')
207 bumps_changed = temperatures_changed = vibrations_changed = False
208 if not isinstance(group, _h5py.Group) and not dry_run:
209 f = _h5py.File(filename, mode)
210 group = _h5_create_group(f, group)
214 for i,bump in enumerate(raw_data['bump']):
215 data['bump'][i],changed = check_bump(
216 index=i, bump=bump, z_axis_config=axis_config,
217 deflection_channel_config=input_config, plot=plot,
218 maximum_relative_error=maximum_relative_error)
219 if changed and not dry_run:
221 bump_group = _h5_create_group(group, 'bump/{}'.format(i))
222 _bump_save(group=bump_group, processed=data['bump'][i])
223 for i,temperature in enumerate(raw_data['temperature']):
224 data['temperature'][i],changed = check_temperature(
225 index=i, temperature=temperature,
226 maximum_relative_error=maximum_relative_error)
227 if changed and not dry_run:
228 temperatures_changed = True
229 temperature_group = _h5_create_group(
230 group, 'temperature/{}'.format(i))
232 group=temerature_group, processed=data['temperature'][i])
233 for i,vibration in enumerate(raw_data['vibration']):
234 data['vibration'][i],changed = check_vibration(
235 index=i, vibration=vibration,
236 deflection_channel_config=input_config, plot=plot,
237 maximum_relative_error=maximum_relative_error)
238 if changed and not dry_run:
239 vibrations_changed = True
240 vibration_group = _h5_create_group(
241 group, 'vibration/{}'.format(i))
243 group=vibration_group, processed=data['vibration'])
244 k,k_s,changed = check_calibration(
245 k=data['processed']['spring_constant'],
246 k_s=data['processed']['spring_constant_deviation'],
248 temperatures=data['temperature'], vibrations=data['vibration'],
249 maximum_relative_error=maximum_relative_error)
250 if (changed or bumps_changed or temperatures_changed or
251 vibrations_changed) and not dry_run:
252 calibration_group = _h5_create_group(group, 'calibration')
254 calib_save(group=calibration_group, bump=data['bump'])
255 if temperatures_changed:
257 group=calibration_group, temperature=data['temperature'])
258 if vibrations_changed:
260 group=calibration_group, vibration=data['vibration'])
262 calib_save(group=calibration_group, k=k, k_s=k_s)
267 _plot(bumps=data['raw']['bump'],
268 temperatures=data['raw']['temperature'],
269 vibrations=data['raw']['vibration'])
272 def check_bump(index, bump, maximum_relative_error, **kwargs):
274 sensitivity = _bump_analyze(
275 config=bump['config']['bump'], data=bump['raw'], **kwargs)
276 if bump.get('processed', None) is None:
278 _LOG.warn('new analysis for bump {}: {}'.format(index, sensitivity))
280 rel_error = abs(sensitivity - bump['processed'])/bump['processed']
281 if rel_error > maximum_relative_error:
283 _LOG.warn(("new analysis doesn't match for bump {}: {} -> {} "
284 "(difference: {}, relative error: {})").format(
285 index, bump['processed'], sensitivity,
286 sensitivity-bump['processed'], rel_error))
287 return (sensitivity, changed)
289 def check_temperature(index, temperature, maximum_relative_error, **kwargs):
291 temp = _temperature_analyze(
292 config=temperature['config']['temperature'],
293 temperature=temperature['raw'], **kwargs)
294 if temperature.get('processed', None) is None:
296 _LOG.warn('new analysis for temperature {}: {}'.format(index, temp))
298 rel_error = abs(temp - temperature['processed']
299 )/temperature['processed']
300 if rel_error > maximum_relative_error:
302 _LOG.warn(("new analysis doesn't match for temperature "
303 "{} -> {} (difference: {}, relative error: {})"
305 index, temperature['processed'], temp,
306 temp-temperature['processed'], rel_error))
307 return (temp, changed)
309 def check_vibration(index, vibration, maximum_relative_error, **kwargs):
311 variance = _vibration_analyze(
312 config=vibration['config']['vibration'],
313 deflection=vibration['raw'], **kwargs)
314 if vibration.get('processed', None) is None:
316 _LOG.warn('new analysis for temperature {}: {}'.format(
319 rel_error = abs(variance-vibration['processed'])/vibration['processed']
320 if rel_error > maximum_relative_error:
321 _LOG.warn(("new analysis doesn't match for vibration {}: {} != {} "
322 "(difference: {}, relative error: {})").format(
323 index, variance, vibration['processed'],
324 variance-vibration['processed'], rel_error))
325 return (variance, changed)
327 def check_calibration(k, k_s, maximum_relative_error, **kwargs):
329 new_k,new_k_s = analyze(**kwargs)
332 _LOG.warn('new analysis for the spring constant: {}'.format(new_k))
334 rel_error = abs(new_k-k)/k
335 if rel_error > maximum_relative_error:
336 _LOG.warn(("new analysis doesn't match for the spring constant: "
337 "{} != {} (difference: {}, relative error: {})").format(
338 new_k, k, new_k-k, rel_error))
341 _LOG.warn('new analysis for the spring constant deviation: {}'.format(
344 rel_error = abs(new_k-k)/k
345 if rel_error > maximum_relative_error:
347 ("new analysis doesn't match for the spring constant deviation"
348 ": {} != {} (difference: {}, relative error: {})").format(
349 new_k_s, k_s, new_k_s-k_s, rel_error))
350 return (new_k, new_k_s, changed)