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 `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 >>> 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['deflection'][:50] = 50
73 >>> processed = analyze(
74 ... config=config, data=raw, z_axis_config=z_axis_config,
75 ... deflection_channel_config=deflection_channel_config)
76 >>> plot(data=raw) # TODO: convert to V and m
77 >>> save(filename=filename, group='/bump/',
78 ... config=config, raw=raw, processed=processed)
80 >>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
85 <HDF5 dataset "far-steps": shape (), type "<i4">
87 <HDF5 dataset "initial-position": shape (), type "<f8">
89 <HDF5 dataset "min-slope-ratio": shape (), type "<f8">
91 <HDF5 dataset "model": shape (), type "|S9">
93 <HDF5 dataset "push-depth": shape (), type "<f8">
95 <HDF5 dataset "push-speed": shape (), type "<f8">
97 <HDF5 dataset "samples": shape (), type "<i4">
99 <HDF5 dataset "setpoint": shape (), type "<f8">
102 <HDF5 dataset "data": shape (), type "<f8">
104 <HDF5 dataset "units": shape (), type "|S3">
108 <HDF5 dataset "data": shape (100,), type "<u2">
109 [50 50 ... 50 51 52 ... 97 98 99]
110 <HDF5 dataset "units": shape (), type "|S4">
113 <HDF5 dataset "data": shape (100,), type "<u2">
114 [ 0 1 2 3 ... 97 98 99]
115 <HDF5 dataset "units": shape (), type "|S4">
118 >>> data = load(filename=filename, group='/bump/')
120 >>> pprint(data) # doctest: +ELLIPSIS, +REPORT_UDIFF
121 {'config': {'bump': <BumpConfig ...>},
122 'processed': 1.00...,
123 'raw': {'deflection': array([50, 50, ..., 52, 53, ..., 98, 99], dtype=uint16),
124 'z': array([ 0, 1, 2, ..., 97, 98, 99], dtype=uint16)}}
126 >>> os.remove(filename)
129 import numpy as _numpy
130 from scipy.optimize import leastsq as _leastsq
133 import matplotlib as _matplotlib
134 import matplotlib.pyplot as _matplotlib_pyplot
135 import time as _time # for timestamping lines on plots
136 except (ImportError, RuntimeError), e:
138 _matplotlib_import_error = e
140 from pypiezo.base import convert_bits_to_volts as _convert_bits_to_volts
141 from pypiezo.base import convert_bits_to_meters as _convert_bits_to_meters
142 from pypiezo.config import AxisConfig as _AxisConfig
143 from pypiezo.config import InputChannelConfig as _InputChannelConfig
145 from . import LOG as _LOG
146 from . import package_config as _package_config
147 from .config import Linear as _Linear
148 from .config import Quadratic as _Quadratic
149 from .config import BumpConfig as _BumpConfig
150 from .util import SaveSpec as _SaveSpec
151 from .util import save as _save
152 from .util import load as _load
155 def analyze(config, data, z_axis_config,
156 deflection_channel_config, plot=False):
157 """Return the slope of the bump.
160 data dictionary of data in DAC/ADC bits
161 config `.config._BumpConfig` instance
162 z_axis_config z `pypiezo.config.AxisConfig` instance
163 deflection_channel_config
164 deflection `pypiezo.config.InputChannelConfig` instance
165 plot boolean overriding matplotlib config setting.
167 photo_sensitivity (Vphoto/Zcant) in Volts/m
169 Checks for strong correlation (r-value) and low randomness chance
172 z = _convert_bits_to_meters(z_axis_config, data['z'])
173 deflection = _convert_bits_to_volts(
174 deflection_channel_config, data['deflection'])
175 high_voltage_rail = _convert_bits_to_volts(
176 deflection_channel_config, deflection_channel_config['maxdata'])
177 if config['model'] == _Linear:
179 'param_guesser': limited_linear_param_guess,
180 'model': limited_linear,
181 'sensitivity_from_fit_params': limited_linear_sensitivity,
185 'param_guesser': limited_quadratic_param_guess,
186 'model': limited_quadratic,
187 'sensitivity_from_fit_params': limited_quadratic_sensitivity,
189 photo_sensitivity = fit(
190 z, deflection, high_voltage_rail=high_voltage_rail, plot=plot,
192 return photo_sensitivity
194 def limited_linear(x, params, high_voltage_rail):
197 flat region (off-surface)
198 linear region (in-contact)
199 flat region (high-voltage-rail)
201 x_contact (x value for the surface-contact kink)
202 y_contact (y value for the surface-contact kink)
203 slope (dy/dx at the surface-contact kink)
205 x_contact,y_contact,slope = params
206 off_surface_mask = x <= x_contact
207 on_surface_mask = x > x_contact
208 y = (off_surface_mask * y_contact +
209 on_surface_mask * (y_contact + slope*(x-x_contact)))
210 y = _numpy.clip(y, y_contact, high_voltage_rail)
213 def limited_linear_param_guess(x, y):
215 Guess rough parameters for a limited_linear model. Assumes the
216 bump approaches (raising the deflection as it does so) first.
217 Retracting after the approach is optional. Approximates the contact
218 position and an on-surface (high) position by finding first crossings
219 of thresholds 0.3 and 0.7 of the y value's total range. Not the
220 most efficient algorithm, but it seems fairly robust.
222 y_contact = float(y.min())
223 y_max = float(y.max())
225 y_low = y_contact + 0.3 * (y_max-y_contact)
226 y_high = y_contact + 0.7 * (y_max-y_contact)
233 x_contact = float(x[i_low])
234 x_high = float(x[i_high])
235 if x_high == x_contact: # things must be pretty flat
236 x_contact = (x_contact + x[0]) / 2
237 slope = (y_high - y_contact) / (x_high - x_contact)
238 return (x_contact, y_contact, slope)
240 def limited_linear_sensitivity(params):
242 Return the estimated sensitivity to small deflections according to
243 limited_linear fit parameters.
248 def limited_quadratic(x, params, high_voltage_rail):
251 flat region (off-surface)
252 quadratic region (in-contact)
253 flat region (high-voltage-rail)
255 x_contact (x value for the surface-contact kink)
256 y_contact (y value for the surface-contact kink)
257 slope (dy/dx at the surface-contact kink)
258 quad (d**2 y / dx**2, allow decreasing sensitivity with increased x)
260 x_contact,y_contact,slope,quad = params
261 off_surface_mask = x <= x_contact
262 on_surface_mask = x > x_contact
263 y = (off_surface_mask * y_contact +
265 y_contact + slope*(x-x_contact) + quad*(x-x_contact)**2))
266 y = _numpy.clip(y, y_contact, high_voltage_rail)
269 def limited_quadratic_param_guess(x, y):
271 Guess rough parameters for a limited_quadratic model. Assumes the
272 bump approaches (raising the deflection as it does so) first.
273 Retracting after the approach is optional. Approximates the contact
274 position and an on-surface (high) position by finding first crossings
275 of thresholds 0.3 and 0.7 of the y value's total range. Not the
276 most efficient algorithm, but it seems fairly robust.
278 x_contact,y_contact,linear_slope = limited_linear_param_guess(x,y)
279 contact_depth = x.max() - x_contact
280 slope = linear_slope / 2
281 quad = slope / contact_depth
282 # The above slope and quad were chosen so
284 # x*slope + x**2 * slope == x * linear_slope
285 return (x_contact, y_contact, slope, quad)
287 def limited_quadratic_sensitivity(params):
289 Return the estimated sensitivity to small deflections according to
290 limited_quadratic fit parameters.
295 def fit(z, deflection, high_voltage_rail,
296 param_guesser=limited_quadratic_param_guess,
297 model=limited_quadratic,
298 sensitivity_from_fit_params=limited_quadratic_sensitivity,
300 """Fit a aurface bump and return the photodiode sensitivitiy.
303 z piezo position in meters
304 deflection photodiode deflection in volts
305 param_guesser function that guesses initial model parameters
306 model parametric model of deflection vs. z
307 sensitivity_from_fit_params given fit params, return the sensitivity
308 plot boolean overriding matplotlib config setting.
310 photodiode_sensitivity photodiode volts per piezo meter
312 _LOG.debug('fit bump data with model %s' % model)
313 def residual(p, deflection, z):
314 return model(z, p, high_voltage_rail=high_voltage_rail) - deflection
315 param_guess = param_guesser(z, deflection)
317 p,cov,info,mesg,ier = _leastsq(
318 residual, param_guess, args=(deflection, z), full_output=True,
321 zd = _numpy.ndarray(list(z.shape) + [2], dtype=z.dtype)
324 _numpy.savetxt('/tmp/z-deflection.dat', zd, delimiter='\t')
326 _LOG.debug('fitted params: %s' % p)
327 _LOG.debug('covariance matrix: %s' % cov)
328 #_LOG.debug('info: %s' % info)
329 _LOG.debug('message: %s' % mesg)
331 _LOG.debug('solution converged')
333 _LOG.debug('solution did not converge')
334 if plot or _package_config['matplotlib']:
335 yguess = model(z, param_guess, high_voltage_rail=high_voltage_rail)
336 yfit = model(z, p, high_voltage_rail=high_voltage_rail)
337 _plot({'z': z, 'deflection': deflection}, yguess=yguess, yfit=yfit)
338 return sensitivity_from_fit_params(p)
340 def save(filename=None, group='/', config=None, z_axis_config=None,
341 deflection_channel_config=None, raw=None, processed=None):
343 _SaveSpec(item=config, relpath='config/bump', config=_BumpConfig),
344 _SaveSpec(item=z_axis_config, relpath='config/z', config=_AxisConfig),
345 _SaveSpec(item=deflection_channel_config, relpath='config/deflection',
346 config=_InputChannelConfig),
347 _SaveSpec(item=processed, relpath='processed', units='V/m'),
350 for key in raw.keys():
351 specs.append(_SaveSpec(
352 item=raw[key], relpath='raw/{}'.format(key), array=True,
354 _save(filename=filename, group=group, specs=specs)
356 def load(filename=None, group='/'):
358 _SaveSpec(key=('config', 'bump'), relpath='config/bump',
360 _SaveSpec(key=('config', 'z_axis_config'), relpath='config/z',
362 _SaveSpec(key=('config', 'deflection_channel_config'),
363 relpath='config/deflection', config=_InputChannelConfig),
364 _SaveSpec(key=('raw', 'z'), relpath='raw/z', array=True, units='bits'),
365 _SaveSpec(key=('raw', 'deflection'), relpath='raw/deflection',
366 array=True, units='bits'),
367 _SaveSpec(key=('processed',), relpath='processed', units='V/m'),
369 return _load(filename=filename, group=group, specs=specs)
371 def plot(data, yguess=None, yfit=None):
372 "Plot the bump (Vphoto vs Vzp)"
374 raise _matplotlib_import_error
375 figure = _matplotlib_pyplot.figure()
377 axes = figure.add_subplot(1, 1, 1)
379 axes = figure.add_subplot(2, 1, 1)
380 residual_axes = figure.add_subplot(2, 1, 2)
381 timestamp = _time.strftime('%H%M%S')
384 axes.plot(data['z'], data['deflection'], '.', label='bump')
386 axes.plot(data['z'], yguess, 'g-', label='guess')
388 axes.plot(data['z'], yfit, 'r-', label='fit')
390 axes.autoscale(tight=True)
391 axes.set_title('bump surface %s' % timestamp)
392 #axes.legend(loc='upper left')
393 axes.set_xlabel('Z piezo (meters)')
394 axes.set_ylabel('Photodiode (Volts)')
396 # second subplot for residual
397 residual_axes.plot(data['z'], data['deflection'] - yfit,
398 'r-', label='residual')
399 residual_axes.autoscale(tight=True)
400 #residual_axes.legend(loc='upper right')
401 residual_axes.set_xlabel('Z piezo (meters)')
402 residual_axes.set_ylabel('Photodiode (Volts)')
403 if hasattr(figure, 'show'):
406 _plot = plot # alternative name for use inside fit()