1 # Copyright (C) 2011-2012 W. Trevor King <wking@drexel.edu>
3 # This file is part of pypiezo.
5 # pypiezo is free software: you can redistribute it and/or modify it under the
6 # terms of the GNU General Public License as published by the Free Software
7 # Foundation, either version 3 of the License, or (at your option) any later
10 # pypiezo is distributed in the hope that it will be useful, but WITHOUT ANY
11 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12 # A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14 # You should have received a copy of the GNU General Public License along with
15 # pypiezo. If not, see <http://www.gnu.org/licenses/>.
17 """Utilities detecting the position of the sample surface.
22 SLEEP_DURING_SURF_POS_AQUISITION = False # doesn't help
24 #from time import sleep as _sleep
25 import numpy as _numpy
26 from scipy.optimize import leastsq as _leastsq
29 import matplotlib as _matplotlib
30 import matplotlib.pyplot as _matplotlib_pyplot
31 import time as _time # for timestamping lines on plots
32 except (ImportError, RuntimeError), e:
34 _matplotlib_import_error = e
36 from . import LOG as _LOG
37 from . import package_config as _package_config
38 from . import base as _base
41 class SurfaceError (Exception):
45 class PoorFit (SurfaceError):
49 class PoorGuess (PoorFit):
53 class FlatFit (PoorFit):
54 "Raised for slope in the non-contact region, or no slope in contact region"
55 def __init__(self, left_slope, right_slope):
56 self.left_slope = left_slope
57 self.right_slope = right_slope
58 msg = 'slopes not sufficiently different: %g and %g' % (
59 left_slope, right_slope)
60 super(FlatFit, self).__init__(msg)
63 class EdgeKink (PoorFit):
64 def __init__(self, kink, edge, window):
68 msg = 'no kink (kink %d not within %d of %d)' % (kink, window, edge)
69 super(EdgeKink, self).__init__(self, msg)
72 def _linspace(*args, **kwargs):
73 dtype = kwargs.pop('dtype')
74 out = _numpy.linspace(*args, **kwargs)
75 return out.reshape((len(out), 1)).astype(dtype)
77 def _get_min_max_positions(piezo, axis_name, min_position=None,
79 output_axis = piezo.axis_by_name(axis_name)
80 if min_position is None:
81 min_position = _base.convert_volts_to_bits(
82 output_axis.config['channel'],
83 output_axis.config['minimum'])
84 if max_position is None:
85 max_position = _base.convert_volts_to_bits(
86 output_axis.config['channel'],
87 output_axis.config['maximum'])
88 return (min_position, max_position)
90 def get_surface_position_data(piezo, axis_name, max_deflection, steps=2000,
91 frequency=10e3, min_position=None,
93 "Measure the distance to the surface"
94 _LOG.info('get surface position')
95 orig_position = piezo.last_output[axis_name]
96 # fully retract the piezo
97 min_position,max_position = _get_min_max_positions(
98 piezo, axis_name, min_position=min_position, max_position=max_position)
99 _LOG.info('retract the piezo to %d' % min_position)
100 dtype = piezo.channel_dtype(axis_name, direction='output')
101 out = _linspace(orig_position, min_position, steps, dtype=dtype)
102 out = out.reshape((len(out), 1)).astype(
103 piezo.channel_dtype(axis_name, direction='output'))
105 'frequency': frequency,
106 'output_names': [axis_name],
107 'input_names': ['deflection'],
109 ret1 = piezo.named_ramp(data=out, **ramp_kwargs)
110 # locate high deflection position
111 _LOG.info('approach until there is dangerous deflection (> %d)'
113 if SLEEP_DURING_SURF_POS_AQUISITION == True:
114 _sleep(.2) # sleeping briefly seems to reduce bounce?
115 mtpod = piezo.move_to_pos_or_def(
116 axis_name=axis_name, position=max_position, deflection=max_deflection,
117 step=(max_position-min_position)/steps, return_data=True)
118 high_contact_pos = piezo.last_output[axis_name]
119 # fully retract the piezo again
120 _LOG.info('retract the piezo to %d again' % min_position)
121 if SLEEP_DURING_SURF_POS_AQUISITION == True:
123 out = _linspace(high_contact_pos, min_position, steps, dtype=dtype)
124 ret2 = piezo.named_ramp(data=out, **ramp_kwargs)
125 # scan to the high contact position
126 _LOG.info('ramp in to the deflected position %d' % high_contact_pos)
127 if SLEEP_DURING_SURF_POS_AQUISITION == True:
129 out = _linspace(min_position, high_contact_pos, steps, dtype=dtype)
130 data = piezo.named_ramp(data=out, **ramp_kwargs)
131 if SLEEP_DURING_SURF_POS_AQUISITION == True:
133 # return to the original position
134 _LOG.info('return to the original position %d' % orig_position)
135 out = _linspace(high_contact_pos, orig_position, steps, dtype=dtype)
136 ret3 = piezo.named_ramp(data=out, **ramp_kwargs)
137 return {'ret1':ret1, 'mtpod':mtpod, 'ret2':ret2,
138 'approach':data, 'ret3':ret3}
140 def bilinear(x, params):
141 """bilinear fit for surface bumps. Model has two linear regimes
142 which meet at x=kink_position and have independent slopes.
144 `x` should be a `numpy.ndarray`.
146 left_offset,left_slope,kink_position,right_slope = params
147 left_mask = x < kink_position
148 right_mask = x >= kink_position # = not left_mask
149 left_y = left_offset + left_slope*x
150 right_y = (left_offset + left_slope*kink_position
151 + right_slope*(x-kink_position))
152 return left_mask * left_y + right_mask * right_y
154 def analyze_surface_position_data(
155 ddict, min_slope_ratio=10.0, kink_window=None,
156 return_all_parameters=False):
159 min_slope_ratio : float
160 Minimum ratio between the non-contact "left" slope and the
161 contact "right" slope.
162 kink_window : int (in bits) or None
163 Raise `EdgeKink` if the kink is within `kink_window` of the
164 minimum or maximum `z` position during the approach. If
165 `None`, a default value of 2% of the approach range is used.
167 # ususes ddict["approach"] for analysis
168 # the others are just along to be plotted
169 _LOG.info('analyze surface position data')
171 data = ddict['approach']
172 # analyze data, using bilinear model
173 # y = p0 + p1 x for x <= p2
174 # = p0 + p1 p2 + p3 (x-p2) for x >= p2
175 dump_before_index = 0 # 25 # HACK!!
176 # Generate a reasonable guess...
177 start_pos = int(data['z'].min())
178 final_pos = int(data['z'].max())
179 start_def = int(data['deflection'].min())
180 final_def = int(data['deflection'].max())
181 # start_def and start_pos are probably for 2 different points
182 _LOG.info('min deflection %d, max deflection %d'
183 % (start_def, final_def))
184 _LOG.info('min position %d, max position %d'
185 % (start_pos, final_pos))
187 left_offset = start_def
189 kink_position = (final_pos+start_pos)/2.0
190 right_slope = 2.0*(final_def-start_def)/(final_pos-start_pos)
191 pstart = [left_offset, left_slope, kink_position, right_slope]
192 _LOG.info('guessed params: %s' % pstart)
194 offset_scale = (final_pos - start_pos)/100
195 left_slope_scale = right_slope/10
196 kink_scale = (final_pos-start_pos)/100
197 right_slope_scale = right_slope
198 scale = [offset_scale, left_slope_scale, kink_scale, right_slope_scale]
199 _LOG.info('guessed scale: %s' % scale)
201 def residual(p, y, x):
204 params,cov,info,mesg,ier = _leastsq(
206 args=(data["deflection"][dump_before_index:],
207 data["z"][dump_before_index:]),
208 full_output=True, maxfev=10000)
209 _LOG.info('best fit parameters: %s' % (params,))
211 if _package_config['matplotlib']:
213 raise _matplotlib_import_error
214 figure = _matplotlib_pyplot.figure()
215 axes = figure.add_subplot(1, 1, 1)
217 timestamp = _time.strftime('%H%M%S')
218 axes.set_title('surf_pos %5g %5g %5g %5g' % tuple(params))
219 def plot_dict(d, label):
220 _pylab.plot(d["z"], d["deflection"], '.',label=label)
221 for n,name in [('ret1', 'first retract'),
222 ('mtpod', 'single step in'),
223 ('ret2', 'second retract'),
224 ('approach', 'main approach'),
225 ('ret3', 'return to start')]:
226 axes.plot(ddict[n]['z'], ddict[n]['deflection'], label=name)
227 def fit_fn(x, params):
229 return params[0] + params[1]*x
231 return (params[0] + params[1]*params[2]
232 + params[3]*(x-params[2]))
233 axes.plot([start_pos, params[2], final_pos],
234 [fit_fn(start_pos, params), fit_fn(params[2], params),
235 fit_fn(final_pos, params)], '-',label='fit')
236 #_pylab.legend(loc='best')
239 # check that the fit is reasonable
240 # params[1] is slope in non-contact region
241 # params[2] is kink position
242 # params[3] is slope in contact region
243 if kink_window is None:
244 kink_window = int(0.02*(final_pos-start_pos))
246 if abs(params[1]*min_slope_ratio) > abs(params[3]):
247 raise FlatFit(left_slope=params[1], right_slope=params[3])
248 if params[2] < start_pos+kink_window:
249 raise EdgeKink(kink=params[2], edge=start_pos, window=kink_window)
250 if params[2] > final_pos-kink_window:
251 raise EdgeKink(kink=params[2], edge=final_pos, window=kink_window)
252 _LOG.info('surface position %s' % params[2])
253 _matplotlib_pyplot.show()
254 if return_all_parameters:
258 def get_surface_position(piezo, axis_name, max_deflection, **kwargs):
259 ddict = get_surface_position_data(piezo, axis_name, max_deflection)
260 return analyze_surface_position_data(ddict, **kwargs)