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 """Surface bump analysis (measures photodiode sensitivity).
23 Separate the more general `bump_analyze()` from the other `bump_*()`
24 functions in calibcant.
26 The relevant physical quantities are:
28 * `Vzp_out` Output z-piezo voltage (what we generate)
29 * `Vzp` Applied z-piezo voltage (after external amplification)
30 * `Zp` The z-piezo position
31 * `Zcant` The cantilever vertical deflection
32 * `Vphoto` The photodiode vertical deflection voltage (what we measure)
34 Which are related by the parameters:
36 * `zp_gain` Vzp_out / Vzp
37 * `zp_sensitivity` Zp / Vzp
38 * `photo_sensitivity` Vphoto / Zcant
40 `photo_sensitivity` is measured by bumping the cantilever against the
41 surface, where `Zp = Zcant` (see `calibrate.bump_acquire()`). The
42 measured slope `Vphoto/Vout` is converted to `photo_sensitivity` with
46 >>> from pprint import pprint
49 >>> from .config import HDF5_BumpConfig
50 >>> from pypiezo.config import pprint_HDF5, HDF5_ChannelConfig, HDF5_AxisConfig
52 >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
55 >>> bump_config = HDF5_BumpConfig(filename=filename, group='/bump/config/')
56 >>> bump_config.save()
57 >>> z_channel_config = HDF5_ChannelConfig(filename=None)
58 >>> z_channel_config['maxdata'] = 200
59 >>> z_channel_config['conversion-coefficients'] = (0,1)
60 >>> z_channel_config['conversion-origin'] = 0
61 >>> z_axis_config = HDF5_AxisConfig(filename=None)
62 >>> deflection_channel_config = HDF5_ChannelConfig(filename=None)
63 >>> deflection_channel_config['maxdata'] = 200
64 >>> deflection_channel_config['conversion-coefficients'] = (0,1)
65 >>> deflection_channel_config['conversion-origin'] = 0
68 ... 'z': numpy.arange(100, dtype=numpy.uint16),
69 ... 'deflection': numpy.arange(100, dtype=numpy.uint16),
71 >>> raw_bump['deflection'][:50] = 50
72 >>> processed_bump = bump_analyze(
73 ... raw_bump, bump_config, z_channel_config, z_axis_config,
74 ... deflection_channel_config)
75 >>> bump_plot(data=raw_bump) # TODO: convert to V and m
76 >>> bump_save(filename=filename, group='/bump/', raw_bump=raw_bump,
77 ... z_channel_config=z_channel_config, z_axis_config=z_axis_config,
78 ... deflection_channel_config=deflection_channel_config,
79 ... processed_bump=processed_bump)
81 >>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
85 /bump/config/deflection
86 /bump/config/deflection/channel
87 <HDF5 dataset "channel": shape (), type "|S1">
89 <HDF5 dataset "conversion-coefficients": shape (), type "|S4">
91 <HDF5 dataset "conversion-origin": shape (), type "|S1">
93 <HDF5 dataset "device": shape (), type "|S12">
95 <HDF5 dataset "inverse-conversion-coefficients": shape (), type "|S1">
97 <HDF5 dataset "inverse-conversion-origin": shape (), type "|S3">
99 <HDF5 dataset "maxdata": shape (), type "|S3">
101 <HDF5 dataset "range": shape (), type "|S1">
103 <HDF5 dataset "subdevice": shape (), type "|S2">
105 <HDF5 dataset "far-steps": shape (), type "|S3">
107 <HDF5 dataset "initial-position": shape (), type "|S6">
109 <HDF5 dataset "model": shape (), type "|S9">
111 <HDF5 dataset "push-depth": shape (), type "|S5">
113 <HDF5 dataset "push-speed": shape (), type "|S5">
115 <HDF5 dataset "samples": shape (), type "|S4">
117 <HDF5 dataset "setpoint": shape (), type "|S3">
121 <HDF5 dataset "gain": shape (), type "|S3">
123 <HDF5 dataset "maximum": shape (), type "|S4">
125 <HDF5 dataset "minimum": shape (), type "|S4">
127 <HDF5 dataset "sensitivity": shape (), type "|S3">
129 /bump/config/z/channel
130 <HDF5 dataset "channel": shape (), type "|S1">
132 <HDF5 dataset "conversion-coefficients": shape (), type "|S4">
134 <HDF5 dataset "conversion-origin": shape (), type "|S1">
136 <HDF5 dataset "device": shape (), type "|S12">
138 <HDF5 dataset "inverse-conversion-coefficients": shape (), type "|S1">
140 <HDF5 dataset "inverse-conversion-origin": shape (), type "|S3">
142 <HDF5 dataset "maxdata": shape (), type "|S3">
144 <HDF5 dataset "range": shape (), type "|S1">
146 <HDF5 dataset "subdevice": shape (), type "|S2">
148 <HDF5 dataset "processed": shape (), type "<f8">
151 <HDF5 dataset "deflection": shape (100,), type "<u2">
152 [50 50 ... 50 51 52 ... 97 98 99]
153 <HDF5 dataset "z": shape (100,), type "<u2">
154 [ 0 1 2 3 ... 97 98 99]
156 >>> (raw_bump,bump_config,z_channel_config,z_axis_config,
157 ... deflection_channel_config,processed_bump) = bump_load(
158 ... filename=filename, group='/bump/')
160 >>> pprint(raw_bump) # doctest: +ELLIPSIS
161 {'deflection': array([50, 50, ... 51, 52, 53, ..., 97, 98, 99], dtype=uint16),
162 'z': array([ 0, 1, 2, ..., 97, 98, 99], dtype=uint16)}
163 >>> processed_bump # doctest: +ELLIPSIS
166 >>> os.remove(filename)
170 import numpy as _numpy
171 from scipy.optimize import leastsq as _leastsq
174 import matplotlib as _matplotlib
175 import matplotlib.pyplot as _matplotlib_pyplot
176 import time as _time # for timestamping lines on plots
177 except (ImportError, RuntimeError), e:
179 _matplotlib_import_error = e
181 from pypiezo.base import convert_bits_to_volts as _convert_bits_to_volts
182 from pypiezo.base import convert_bits_to_meters as _convert_bits_to_meters
183 from pypiezo.config import HDF5_ChannelConfig as _HDF5_ChannelConfig
184 from pypiezo.config import HDF5_AxisConfig as _HDF5_AxisConfig
185 from pypiezo.config import h5_create_group as _h5_create_group
187 from . import LOG as _LOG
188 from . import base_config as _base_config
189 from .config import Linear as _Linear
190 from .config import Quadratic as _Quadratic
191 from .config import HDF5_BumpConfig as _HDF5_BumpConfig
194 def bump_analyze(data, bump_config, z_channel_config, z_axis_config,
195 deflection_channel_config, plot=False):
196 """Return the slope of the bump.
199 data dictionary of data in DAC/ADC bits
200 bump_config `.config._BumpConfig` instance
201 z_channel_config z `pypiezo.config._ChannelConfig` instance
202 z_axis_config z `pypiezo.config._AxisConfig` instance
203 deflection_channel_config
204 deflection `pypiezo.config._ChannelConfig` instance
205 plot boolean overriding matplotlib config setting.
207 photo_sensitivity (Vphoto/Zcant) in Volts/m
209 Checks for strong correlation (r-value) and low randomness chance
212 z = _convert_bits_to_meters(z_channel_config, z_axis_config, data['z'])
213 deflection = _convert_bits_to_volts(
214 deflection_channel_config, data['deflection'])
215 if bump_config['model'] == _Linear:
217 'param_guesser': limited_linear_param_guess,
218 'model': limited_linear,
219 'sensitivity_from_fit_params': limited_linear_sensitivity,
223 'param_guesser': limited_quadratic_param_guess,
224 'model': limited_quadratic,
225 'sensitivity_from_fit_params': limited_quadratic_sensitivity,
227 photo_sensitivity = bump_fit(z, deflection, plot=plot, **kwargs)
228 return photo_sensitivity
230 def limited_linear(x, params):
233 flat region (off-surface)
234 linear region (in-contact)
235 flat region (high-voltage-rail)
237 x_contact (x value for the surface-contact kink)
238 y_contact (y value for the surface-contact kink)
239 slope (dy/dx at the surface-contact kink)
241 high_voltage_rail = 2**16 - 1 # bits
242 x_contact,y_contact,slope = params
243 y = slope*(x-x_contact) + y_contact
244 y = _numpy.clip(y, y_contact, high_voltage_rail)
247 def limited_linear_param_guess(x, y):
249 Guess rough parameters for a limited_linear model. Assumes the
250 bump approaches (raising the deflection as it does so) first.
251 Retracting after the approach is optional. Approximates the contact
252 position and an on-surface (high) position by finding first crossings
253 of thresholds 0.3 and 0.7 of the y value's total range. Not the
254 most efficient algorithm, but it seems fairly robust.
256 y_contact = float(y.min())
257 y_max = float(y.max())
259 y_low = y_contact + 0.3 * (y_max-y_contact)
260 y_high = y_contact + 0.7 * (y_max-y_contact)
267 x_contact = float(x[i_low])
268 x_high = float(x[i_high])
269 slope = (y_high - y_contact) / (x_high - x_contact)
270 return (x_contact, y_contact, slope)
272 def limited_linear_sensitivity(params):
274 Return the estimated sensitivity to small deflections according to
275 limited_linear fit parameters.
280 def limited_quadratic(x, params):
283 flat region (off-surface)
284 quadratic region (in-contact)
285 flat region (high-voltage-rail)
287 x_contact (x value for the surface-contact kink)
288 y_contact (y value for the surface-contact kink)
289 slope (dy/dx at the surface-contact kink)
290 quad (d**2 y / dx**2, allow decreasing sensitivity with increased x)
292 high_voltage_rail = 2**16 - 1 # bits
293 x_contact,y_contact,slope,quad = params
294 y = slope*(x-x_contact) + quad*(x-x_contact)**2+ y_contact
295 y = _numpy.clip(y, y_contact, high_voltage_rail)
298 def limited_quadratic_param_guess(x, y):
300 Guess rough parameters for a limited_quadratic model. Assumes the
301 bump approaches (raising the deflection as it does so) first.
302 Retracting after the approach is optional. Approximates the contact
303 position and an on-surface (high) position by finding first crossings
304 of thresholds 0.3 and 0.7 of the y value's total range. Not the
305 most efficient algorithm, but it seems fairly robust.
307 x_contact,y_contact,linear_slope = limited_linear_param_guess(x,y)
308 contact_depth = x.max() - x_contact
309 slope = linear_slope / 2
310 quad = slope / contact_depth
311 # The above slope and quad were chosen so
313 # x*slope + x**2 * slope == x * linear_slope
314 return (x_contact, y_contact, slope, quad)
316 def limited_quadratic_sensitivity(params):
318 Return the estimated sensitivity to small deflections according to
319 limited_quadratic fit parameters.
324 def bump_fit(z, deflection,
325 param_guesser=limited_quadratic_param_guess,
326 model=limited_quadratic,
327 sensitivity_from_fit_params=limited_quadratic_sensitivity,
329 """Fit a aurface bump and return the photodiode sensitivitiy.
332 z piezo position in meters
333 deflection photodiode deflection in volts
334 param_guesser function that guesses initial model parameters
335 model parametric model of deflection vs. z
336 sensitivity_from_fit_params given fit params, return the sensitivity
337 plot boolean overriding matplotlib config setting.
339 photodiode_sensitivity photodiode volts per piezo meter
341 _LOG.debug('fit bump data with model %s' % model)
342 def residual(p, deflection, z):
343 return model(z, p) - deflection
344 param_guess = param_guesser(z, deflection)
345 p,cov,info,mesg,ier = _leastsq(
346 residual, param_guess, args=(deflection, z), full_output=True,
348 _LOG.debug('fitted params: %s' % p)
349 _LOG.debug('covariance matrix: %s' % cov)
350 #_LOG.debug('info: %s' % info)
351 _LOG.debug('message: %s' % mesg)
353 _LOG.debug('solution converged')
355 _LOG.debug('solution did not converge')
356 if plot or _base_config['matplotlib']:
357 yguess = model(z, param_guess)
359 bump_plot({'z': z, 'deflection': deflection}, yguess=yguess, yfit=yfit)
360 return sensitivity_from_fit_params(p)
362 def bump_save(filename, group='/', raw_bump=None, bump_config=None,
363 z_channel_config=None, z_axis_config=None,
364 deflection_channel_config=None, processed_bump=None):
365 with _h5py.File(filename, 'a') as f:
366 cwg = _h5_create_group(f, group)
367 if raw_bump is not None:
373 del cwg['raw/deflection']
376 cwg['raw/z'] = raw_bump['z']
377 cwg['raw/deflection'] = raw_bump['deflection']
378 for config,key in [(bump_config, 'config/bump'),
379 (z_channel_config, 'config/z/channel'),
380 (z_axis_config, 'config/z/axis'),
381 (deflection_channel_config,
382 'config/deflection/channel')]:
385 config_cwg = _h5_create_group(cwg, key)
386 config.save(group=config_cwg)
387 if processed_bump is not None:
392 cwg['processed'] = processed_bump
394 def bump_load(filename, group='/'):
395 assert group.endswith('/')
396 raw_bump = processed_bump = None
398 with _h5py.File(filename, 'a') as f:
401 'z': f[group+'raw/z'][...],
402 'deflection': f[group+'raw/deflection'][...],
406 for Config,key in [(_HDF5_BumpConfig, 'config/bump'),
407 (_HDF5_ChannelConfig, 'config/z/channel'),
408 (_HDF5_AxisConfig, 'config/z/axis'),
409 (_HDF5_ChannelConfig, 'config/deflection/channel')]:
410 config = Config(filename=filename, group=group+key)
411 configs.append(config)
413 processed_bump = f[group+'processed'][...]
418 ret.append(processed_bump)
419 for config in configs:
423 def bump_plot(data, yguess=None, yfit=None):
424 "Plot the bump (Vphoto vs Vzp)"
426 raise _matplotlib_import_error
427 figure = _matplotlib_pyplot.figure()
429 axes = figure.add_subplot(1, 1, 1)
431 axes = figure.add_subplot(2, 1, 1)
432 residual_axes = figure.add_subplot(2, 1, 2)
433 timestamp = _time.strftime('%H%M%S')
436 axes.plot(data['z'], data['deflection'], '.', label='bump')
438 axes.plot(data['z'], yguess, 'g-', label='guess')
440 axes.plot(data['z'], yfit, 'r-', label='fit')
442 axes.set_title('bump surface %s' % timestamp)
443 #axes.legend(loc='upper left')
444 axes.set_xlabel('Z piezo (meters)')
445 axes.set_ylabel('Photodiode (Volts)')
447 # second subplot for residual
448 residual_axes.plot(data['z'], data['deflection'] - yfit,
449 'r-', label='residual')
450 #residual_axes.legend(loc='upper right')
451 residual_axes.set_xlabel('Z piezo (meters)')
452 residual_axes.set_ylabel('Photodiode (Volts)')