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 """Surface bump analysis (measures photodiode sensitivity).
21 Separate the more general `bump_analyze()` from the other `bump_*()`
22 functions in calibcant.
24 The relevant physical quantities are:
26 * `Vzp_out` Output z-piezo voltage (what we generate)
27 * `Vzp` Applied z-piezo voltage (after external amplification)
28 * `Zp` The z-piezo position
29 * `Zcant` The cantilever vertical deflection
30 * `Vphoto` The photodiode vertical deflection voltage (what we measure)
32 Which are related by the parameters:
34 * `zp_gain` Vzp_out / Vzp
35 * `zp_sensitivity` Zp / Vzp
36 * `photo_sensitivity` Vphoto / Zcant
38 `photo_sensitivity` is measured by bumping the cantilever against the
39 surface, where `Zp = Zcant` (see `calibrate.bump_acquire()`). The
40 measured slope `Vphoto/Vout` is converted to `photo_sensitivity` with
44 >>> from pprint import pprint
47 >>> from h5config.storage.hdf5 import pprint_HDF5
48 >>> from pypiezo.config import ChannelConfig, AxisConfig
49 >>> from .config import BumpConfig
51 >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
54 >>> bump_config = BumpConfig()
55 >>> z_channel_config = ChannelConfig()
56 >>> z_channel_config['name'] = 'z'
57 >>> z_channel_config['maxdata'] = 200
58 >>> z_channel_config['conversion-coefficients'] = (0,1)
59 >>> z_channel_config['conversion-origin'] = 0
60 >>> z_axis_config = AxisConfig()
61 >>> z_axis_config['channel'] = z_channel_config
62 >>> deflection_channel_config = ChannelConfig()
63 >>> deflection_channel_config['name'] = 'deflection'
64 >>> deflection_channel_config['maxdata'] = 200
65 >>> deflection_channel_config['conversion-coefficients'] = (0,1)
66 >>> deflection_channel_config['conversion-origin'] = 0
69 ... 'z': numpy.arange(100, dtype=numpy.uint16),
70 ... 'deflection': numpy.arange(100, dtype=numpy.uint16),
72 >>> raw_bump['deflection'][:50] = 50
73 >>> processed_bump = bump_analyze(
74 ... raw_bump, bump_config, z_axis_config, 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 ... bump_config=bump_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
86 <HDF5 dataset "far-steps": shape (), type "<i4">
88 <HDF5 dataset "initial-position": shape (), type "<f8">
90 <HDF5 dataset "model": shape (), type "|S9">
92 <HDF5 dataset "push-depth": shape (), type "<f8">
94 <HDF5 dataset "push-speed": shape (), type "<f8">
96 <HDF5 dataset "samples": shape (), type "<i4">
98 <HDF5 dataset "setpoint": shape (), type "<f8">
100 /bump/config/deflection
101 /bump/config/deflection/channel
102 <HDF5 dataset "channel": shape (), type "<i4">
104 <HDF5 dataset "conversion-coefficients": shape (2,), type "<i4">
106 <HDF5 dataset "conversion-origin": shape (), type "<i4">
108 <HDF5 dataset "device": shape (), type "|S12">
110 <HDF5 dataset "inverse-conversion-coefficients": shape (), type "|S1">
112 <HDF5 dataset "inverse-conversion-origin": shape (), type "<f8">
114 <HDF5 dataset "maxdata": shape (), type "<i4">
116 <HDF5 dataset "name": shape (), type "|S10">
118 <HDF5 dataset "range": shape (), type "<i4">
120 <HDF5 dataset "subdevice": shape (), type "<i4">
124 /bump/config/z/axis/channel
125 <HDF5 dataset "channel": shape (), type "<i4">
127 <HDF5 dataset "conversion-coefficients": shape (2,), type "<i4">
129 <HDF5 dataset "conversion-origin": shape (), type "<i4">
131 <HDF5 dataset "device": shape (), type "|S12">
133 <HDF5 dataset "inverse-conversion-coefficients": shape (), type "|S1">
135 <HDF5 dataset "inverse-conversion-origin": shape (), type "<f8">
137 <HDF5 dataset "maxdata": shape (), type "<i4">
139 <HDF5 dataset "name": shape (), type "|S1">
141 <HDF5 dataset "range": shape (), type "<i4">
143 <HDF5 dataset "subdevice": shape (), type "<i4">
145 <HDF5 dataset "gain": shape (), type "<f8">
147 <HDF5 dataset "maximum": shape (), type "|S4">
149 <HDF5 dataset "minimum": shape (), type "|S4">
151 <HDF5 dataset "monitor": shape (), type "|S1">
153 <HDF5 dataset "sensitivity": shape (), type "<f8">
155 <HDF5 dataset "processed": shape (), type "<f8">
158 <HDF5 dataset "deflection": shape (100,), type "<u2">
159 [50 50 ... 50 51 52 ... 97 98 99]
160 <HDF5 dataset "z": shape (100,), type "<u2">
161 [ 0 1 2 3 ... 97 98 99]
163 >>> (raw_bump,bump_config,z_axis_config,deflection_channel_config,
164 ... processed_bump) = bump_load(filename=filename, group='/bump/')
166 >>> pprint(raw_bump) # doctest: +ELLIPSIS
167 {'deflection': array([50, 50, ... 51, 52, 53, ..., 97, 98, 99], dtype=uint16),
168 'z': array([ 0, 1, 2, ..., 97, 98, 99], dtype=uint16)}
169 >>> processed_bump # doctest: +ELLIPSIS
172 >>> os.remove(filename)
176 import numpy as _numpy
177 from scipy.optimize import leastsq as _leastsq
180 import matplotlib as _matplotlib
181 import matplotlib.pyplot as _matplotlib_pyplot
182 import time as _time # for timestamping lines on plots
183 except (ImportError, RuntimeError), e:
185 _matplotlib_import_error = e
187 from h5config.storage.hdf5 import HDF5_Storage as _HDF5_Storage
188 from h5config.storage.hdf5 import h5_create_group as _h5_create_group
189 from pypiezo.base import convert_bits_to_volts as _convert_bits_to_volts
190 from pypiezo.base import convert_bits_to_meters as _convert_bits_to_meters
191 from pypiezo.config import ChannelConfig as _ChannelConfig
192 from pypiezo.config import AxisConfig as _AxisConfig
194 from . import LOG as _LOG
195 from . import package_config as _package_config
196 from .config import Linear as _Linear
197 from .config import Quadratic as _Quadratic
198 from .config import BumpConfig as _BumpConfig
201 def bump_analyze(data, bump_config, z_axis_config,
202 deflection_channel_config, plot=False):
203 """Return the slope of the bump.
206 data dictionary of data in DAC/ADC bits
207 bump_config `.config._BumpConfig` instance
208 z_axis_config z `pypiezo.config.AxisConfig` instance
209 deflection_channel_config
210 deflection `pypiezo.config.ChannelConfig` instance
211 plot boolean overriding matplotlib config setting.
213 photo_sensitivity (Vphoto/Zcant) in Volts/m
215 Checks for strong correlation (r-value) and low randomness chance
218 z = _convert_bits_to_meters(z_axis_config, data['z'])
219 deflection = _convert_bits_to_volts(
220 deflection_channel_config, data['deflection'])
221 high_voltage_rail = _convert_bits_to_volts(
222 deflection_channel_config, deflection_channel_config['maxdata'])
223 if bump_config['model'] == _Linear:
225 'param_guesser': limited_linear_param_guess,
226 'model': limited_linear,
227 'sensitivity_from_fit_params': limited_linear_sensitivity,
231 'param_guesser': limited_quadratic_param_guess,
232 'model': limited_quadratic,
233 'sensitivity_from_fit_params': limited_quadratic_sensitivity,
235 photo_sensitivity = bump_fit(
236 z, deflection, high_voltage_rail=high_voltage_rail, plot=plot,
238 return photo_sensitivity
240 def limited_linear(x, params, high_voltage_rail):
243 flat region (off-surface)
244 linear region (in-contact)
245 flat region (high-voltage-rail)
247 x_contact (x value for the surface-contact kink)
248 y_contact (y value for the surface-contact kink)
249 slope (dy/dx at the surface-contact kink)
251 x_contact,y_contact,slope = params
252 off_surface_mask = x <= x_contact
253 on_surface_mask = x > x_contact
254 y = (off_surface_mask * y_contact +
255 on_surface_mask * (y_contact + slope*(x-x_contact)))
256 y = _numpy.clip(y, y_contact, high_voltage_rail)
259 def limited_linear_param_guess(x, y):
261 Guess rough parameters for a limited_linear model. Assumes the
262 bump approaches (raising the deflection as it does so) first.
263 Retracting after the approach is optional. Approximates the contact
264 position and an on-surface (high) position by finding first crossings
265 of thresholds 0.3 and 0.7 of the y value's total range. Not the
266 most efficient algorithm, but it seems fairly robust.
268 y_contact = float(y.min())
269 y_max = float(y.max())
271 y_low = y_contact + 0.3 * (y_max-y_contact)
272 y_high = y_contact + 0.7 * (y_max-y_contact)
279 x_contact = float(x[i_low])
280 x_high = float(x[i_high])
281 if x_high == x_contact: # things must be pretty flat
282 x_contact = (x_contact + x[0]) / 2
283 slope = (y_high - y_contact) / (x_high - x_contact)
284 return (x_contact, y_contact, slope)
286 def limited_linear_sensitivity(params):
288 Return the estimated sensitivity to small deflections according to
289 limited_linear fit parameters.
294 def limited_quadratic(x, params, high_voltage_rail):
297 flat region (off-surface)
298 quadratic region (in-contact)
299 flat region (high-voltage-rail)
301 x_contact (x value for the surface-contact kink)
302 y_contact (y value for the surface-contact kink)
303 slope (dy/dx at the surface-contact kink)
304 quad (d**2 y / dx**2, allow decreasing sensitivity with increased x)
306 x_contact,y_contact,slope,quad = params
307 off_surface_mask = x <= x_contact
308 on_surface_mask = x > x_contact
309 y = (off_surface_mask * y_contact +
311 y_contact + slope*(x-x_contact) + quad*(x-x_contact)**2))
312 y = _numpy.clip(y, y_contact, high_voltage_rail)
315 def limited_quadratic_param_guess(x, y):
317 Guess rough parameters for a limited_quadratic model. Assumes the
318 bump approaches (raising the deflection as it does so) first.
319 Retracting after the approach is optional. Approximates the contact
320 position and an on-surface (high) position by finding first crossings
321 of thresholds 0.3 and 0.7 of the y value's total range. Not the
322 most efficient algorithm, but it seems fairly robust.
324 x_contact,y_contact,linear_slope = limited_linear_param_guess(x,y)
325 contact_depth = x.max() - x_contact
326 slope = linear_slope / 2
327 quad = slope / contact_depth
328 # The above slope and quad were chosen so
330 # x*slope + x**2 * slope == x * linear_slope
331 return (x_contact, y_contact, slope, quad)
333 def limited_quadratic_sensitivity(params):
335 Return the estimated sensitivity to small deflections according to
336 limited_quadratic fit parameters.
341 def bump_fit(z, deflection, high_voltage_rail,
342 param_guesser=limited_quadratic_param_guess,
343 model=limited_quadratic,
344 sensitivity_from_fit_params=limited_quadratic_sensitivity,
346 """Fit a aurface bump and return the photodiode sensitivitiy.
349 z piezo position in meters
350 deflection photodiode deflection in volts
351 param_guesser function that guesses initial model parameters
352 model parametric model of deflection vs. z
353 sensitivity_from_fit_params given fit params, return the sensitivity
354 plot boolean overriding matplotlib config setting.
356 photodiode_sensitivity photodiode volts per piezo meter
358 _LOG.debug('fit bump data with model %s' % model)
359 def residual(p, deflection, z):
360 return model(z, p, high_voltage_rail=high_voltage_rail) - deflection
361 param_guess = param_guesser(z, deflection)
362 p,cov,info,mesg,ier = _leastsq(
363 residual, param_guess, args=(deflection, z), full_output=True,
365 _LOG.debug('fitted params: %s' % p)
366 _LOG.debug('covariance matrix: %s' % cov)
367 #_LOG.debug('info: %s' % info)
368 _LOG.debug('message: %s' % mesg)
370 _LOG.debug('solution converged')
372 _LOG.debug('solution did not converge')
373 if plot or _package_config['matplotlib']:
374 yguess = model(z, param_guess, high_voltage_rail=high_voltage_rail)
375 yfit = model(z, p, high_voltage_rail=high_voltage_rail)
376 bump_plot({'z': z, 'deflection': deflection}, yguess=yguess, yfit=yfit)
377 return sensitivity_from_fit_params(p)
379 def bump_save(filename, group='/', raw_bump=None, bump_config=None,
380 z_axis_config=None, deflection_channel_config=None,
381 processed_bump=None):
382 with _h5py.File(filename, 'a') as f:
383 cwg = _h5_create_group(f, group)
384 if raw_bump is not None:
385 for k in ['z', 'deflection']:
387 del cwg['raw/{}'.format(k)]
390 cwg['raw/z'] = raw_bump['z']
391 cwg['raw/deflection'] = raw_bump['deflection']
392 storage = _HDF5_Storage()
393 for config,key in [(bump_config, 'config/bump'),
394 (z_axis_config, 'config/z/axis'),
395 (deflection_channel_config,
396 'config/deflection/channel')]:
399 config_cwg = _h5_create_group(cwg, key)
400 storage.save(config=config, group=config_cwg)
401 if processed_bump is not None:
406 cwg['processed'] = processed_bump
408 def bump_load(filename, group='/'):
409 assert group.endswith('/')
410 raw_bump = processed_bump = None
412 with _h5py.File(filename, 'a') as f:
415 'z': f[group+'raw/z'][...],
416 'deflection': f[group+'raw/deflection'][...],
420 for Config,key in [(_BumpConfig, 'config/bump'),
421 (_AxisConfig, 'config/z/axis'),
422 (_ChannelConfig, 'config/deflection/channel')]:
423 config = Config(storage=_HDF5_Storage(
424 filename=filename, group=group+key))
425 configs.append(config)
427 processed_bump = float(f[group+'processed'][...])
432 ret.append(processed_bump)
433 for config in configs:
437 def bump_plot(data, yguess=None, yfit=None):
438 "Plot the bump (Vphoto vs Vzp)"
440 raise _matplotlib_import_error
441 figure = _matplotlib_pyplot.figure()
443 axes = figure.add_subplot(1, 1, 1)
445 axes = figure.add_subplot(2, 1, 1)
446 residual_axes = figure.add_subplot(2, 1, 2)
447 timestamp = _time.strftime('%H%M%S')
450 axes.plot(data['z'], data['deflection'], '.', label='bump')
452 axes.plot(data['z'], yguess, 'g-', label='guess')
454 axes.plot(data['z'], yfit, 'r-', label='fit')
456 axes.set_title('bump surface %s' % timestamp)
457 #axes.legend(loc='upper left')
458 axes.set_xlabel('Z piezo (meters)')
459 axes.set_ylabel('Photodiode (Volts)')
461 # second subplot for residual
462 residual_axes.plot(data['z'], data['deflection'] - yfit,
463 'r-', label='residual')
464 #residual_axes.legend(loc='upper right')
465 residual_axes.set_xlabel('Z piezo (meters)')
466 residual_axes.set_ylabel('Photodiode (Volts)')
467 if hasattr(figure, 'show'):