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'):
186 _plot = plot # alternative name for use inside analyze_all()
189 def analyze_all(config, data, raw_data, maximum_relative_error=1e-5,
190 filename=None, group=None, plot=False, dry_run=False):
191 "(Re)analyze (and possibly plot) all data from a `calib()` run."
192 if not data.get('bump', None):
193 data['bump'] = _numpy.zeros((config['num-bumps'],), dtype=float)
194 if not data.get('temperature', None):
195 data['temperature'] = _numpy.zeros(
196 (config['num-temperatures'],), dtype=float)
197 if not data.get('vibrations', None):
198 data['vibration'] = _numpy.zeros(
199 (config['num-vibrations'],), dtype=float)
200 axis_config = config['afm']['piezo'].select_config(
202 attribute_value=config['afm']['main-axis'],
203 get_attribute=_get_axis_name)
204 input_config = config['afm']['piezo'].select_config(
205 setting_name='inputs', attribute_value='deflection')
206 bumps_changed = temperatures_changed = vibrations_changed = False
207 if not isinstance(group, _h5py.Group) and not dry_run:
208 f = _h5py.File(filename, mode)
209 group = _h5_create_group(f, group)
213 for i,bump in enumerate(raw_data['bump']):
214 data['bump'][i],changed = check_bump(
215 index=i, bump=bump, z_axis_config=axis_config,
216 deflection_channel_config=input_config, plot=plot,
217 maximum_relative_error=maximum_relative_error)
218 if changed and not dry_run:
220 bump_group = _h5_create_group(group, 'bump/{}'.format(i))
221 _bump_save(group=bump_group, processed=data['bump'][i])
222 for i,temperature in enumerate(raw_data['temperature']):
223 data['temperature'][i],changed = check_temperature(
224 index=i, temperature=temperature,
225 maximum_relative_error=maximum_relative_error)
226 if changed and not dry_run:
227 temperatures_changed = True
228 temperature_group = _h5_create_group(
229 group, 'temperature/{}'.format(i))
231 group=temerature_group, processed=data['temperature'][i])
232 for i,vibration in enumerate(raw_data['vibration']):
233 data['vibration'][i],changed = check_vibration(
234 index=i, vibration=vibration,
235 deflection_channel_config=input_config, plot=plot,
236 maximum_relative_error=maximum_relative_error)
237 if changed and not dry_run:
238 vibrations_changed = True
239 vibration_group = _h5_create_group(
240 group, 'vibration/{}'.format(i))
242 group=vibration_group, processed=data['vibration'])
243 k,k_s,changed = check_calibration(
244 k=data['processed']['spring_constant'],
245 k_s=data['processed']['spring_constant_deviation'],
247 temperatures=data['temperature'], vibrations=data['vibration'],
248 maximum_relative_error=maximum_relative_error)
249 if (changed or bumps_changed or temperatures_changed or
250 vibrations_changed) and not dry_run:
251 calibration_group = _h5_create_group(group, 'calibration')
253 calib_save(group=calibration_group, bump=data['bump'])
254 if temperatures_changed:
256 group=calibration_group, temperature=data['temperature'])
257 if vibrations_changed:
259 group=calibration_group, vibration=data['vibration'])
261 calib_save(group=calibration_group, k=k, k_s=k_s)
266 _plot(bumps=data['raw']['bump'],
267 temperatures=data['raw']['temperature'],
268 vibrations=data['raw']['vibration'])
271 def check_bump(index, bump, maximum_relative_error, **kwargs):
273 sensitivity = _bump_analyze(
274 config=bump['config']['bump'], data=bump['raw'], **kwargs)
275 if bump.get('processed', None) is None:
277 _LOG.warn('new analysis for bump {}: {}'.format(index, sensitivity))
279 rel_error = abs(sensitivity - bump['processed'])/bump['processed']
280 if rel_error > maximum_relative_error:
282 _LOG.warn(("new analysis doesn't match for bump {}: {} -> {} "
283 "(difference: {}, relative error: {})").format(
284 index, bump['processed'], sensitivity,
285 sensitivity-bump['processed'], rel_error))
286 return (sensitivity, changed)
288 def check_temperature(index, temperature, maximum_relative_error, **kwargs):
290 temp = _temperature_analyze(
291 config=temperature['config']['temperature'],
292 temperature=temperature['raw'], **kwargs)
293 if temperature.get('processed', None) is None:
295 _LOG.warn('new analysis for temperature {}: {}'.format(index, temp))
297 rel_error = abs(temp - temperature['processed']
298 )/temperature['processed']
299 if rel_error > maximum_relative_error:
301 _LOG.warn(("new analysis doesn't match for temperature "
302 "{} -> {} (difference: {}, relative error: {})"
304 index, temperature['processed'], temp,
305 temp-temperature['processed'], rel_error))
306 return (temp, changed)
308 def check_vibration(index, vibration, maximum_relative_error, **kwargs):
310 variance = _vibration_analyze(
311 config=vibration['config']['vibration'],
312 deflection=vibration['raw'], **kwargs)
313 if vibration.get('processed', None) is None:
315 _LOG.warn('new analysis for temperature {}: {}'.format(
318 rel_error = abs(variance-vibration['processed'])/vibration['processed']
319 if rel_error > maximum_relative_error:
320 _LOG.warn(("new analysis doesn't match for vibration {}: {} != {} "
321 "(difference: {}, relative error: {})").format(
322 index, variance, vibration['processed'],
323 variance-vibration['processed'], rel_error))
324 return (variance, changed)
326 def check_calibration(k, k_s, maximum_relative_error, **kwargs):
328 new_k,new_k_s = analyze(**kwargs)
331 _LOG.warn('new analysis for the spring constant: {}'.format(new_k))
333 rel_error = abs(new_k-k)/k
334 if rel_error > maximum_relative_error:
335 _LOG.warn(("new analysis doesn't match for the spring constant: "
336 "{} != {} (difference: {}, relative error: {})").format(
337 new_k, k, new_k-k, rel_error))
340 _LOG.warn('new analysis for the spring constant deviation: {}'.format(
343 rel_error = abs(new_k-k)/k
344 if rel_error > maximum_relative_error:
346 ("new analysis doesn't match for the spring constant deviation"
347 ": {} != {} (difference: {}, relative error: {})").format(
348 new_k_s, k_s, new_k_s-k_s, rel_error))
349 return (new_k, new_k_s, changed)