pypiezo/surface.py

   1 # Copyright (C) 2008-2011 W. Trevor King <wking@drexel.edu>
   2 #
   3 # This file is part of pypiezo.
   4 #
   5 # pypiezo is free software; you can redistribute it and/or modify it
   6 # under the terms of the GNU General Public License as published by the
   7 # Free Software Foundation, either version 3 of the License, or (at your
   8 # option) any later version.
   9 #
  10 # pypiezo is distributed in the hope that it will be useful, but
  11 # WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 # General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU General Public License
  16 # along with pypiezo.  If not, see <http://www.gnu.org/licenses/>.
  17
  18 """Utilities detecting the position of the sample surface.
  19
  20 TODO: doctests
  21 """
  22
  23 SLEEP_DURING_SURF_POS_AQUISITION = False # doesn't help
  24
  25 #from time import sleep as _sleep
  26 import numpy as _numpy
  27 from scipy.optimize import leastsq as _leastsq
  28
  29 try:
  30     import matplotlib as _matplotlib
  31     import matplotlib.pyplot as _matplotlib_pyplot
  32     import time as _time  # for timestamping lines on plots
  33 except ImportError, e:
  34     _matplotlib = None
  35     _matplotlib_import_error = e
  36
  37 from . import LOG as _LOG
  38 from . import base_config as _base_config
  39 from . import base as _base
  40
  41
  42 class SurfaceError (Exception):
  43     pass
  44
  45
  46 class PoorFit (SurfaceError):
  47     pass
  48
  49
  50 class PoorGuess (PoorFit):
  51     pass
  52
  53
  54 class FlatFit (PoorFit):
  55     "Raised for slope in the non-contact region, or no slope in contact region"
  56     def __init__(self, left_slope, right_slope):
  57         self.left_slope = left_slope
  58         self.right_slope = right_slope
  59         msg = 'slopes not sufficiently different: %g and %g' % (
  60             left_slope, right_slope)
  61         super(FlatFit, self).__init__(msg)
  62
  63 class EdgeKink (PoorFit):
  64     def __init__(self, kink, edge, window):
  65         self.kink = kink
  66         self.edge = edge
  67         self.window = window
  68         msg = 'no kink (kink %d not within %d of %d)' % (kink, window, edge)
  69         super(EdgeKink, self).__init__(self, msg)
  70
  71
  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)
  76
  77 def _get_min_max_positions(piezo, axis_name, min_position=None,
  78                            max_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.axis_channel_config,
  83             output_axis.axis_config['minimum'])
  84     if max_position is None:
  85         max_position = _base.convert_volts_to_bits(
  86             output_axis.axis_channel_config,
  87             output_axis.axis_config['maximum'])
  88     return (min_position, max_position)
  89
  90 def get_surface_position_data(piezo, axis_name, max_deflection, steps=2000,
  91                               frequency=10e3, min_position=None,
  92                               max_position=None):
  93     "Measure the distance to the surface"
  94     _LOG.debug('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.debug('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'))
 104     ramp_kwargs = {
 105         'frequency': frequency,
 106         'output_names': [axis_name],
 107         'input_names': ['deflection'],
 108         }
 109     ret1 = piezo.named_ramp(data=out, **ramp_kwargs)
 110     # locate high deflection position
 111     _LOG.debug('approach until there is dangerous deflection (> %d)'
 112                % max_deflection)
 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.debug('retract the piezo to %d again' % min_position)
 121     if SLEEP_DURING_SURF_POS_AQUISITION == True:
 122         _sleep(.2)
 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.debug('ramp in to the deflected position %d' % high_contact_pos)
 127     if SLEEP_DURING_SURF_POS_AQUISITION == True:
 128         _sleep(.2)
 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:
 132         _sleep(.2)
 133     # return to the original position
 134     _LOG.debug('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}
 139
 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.
 143
 144     `x` should be a `numpy.ndarray`.
 145     """
 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
 153
 154 def analyze_surface_position_data(
 155     ddict, min_slope_ratio=10.0, kink_window=None, return_all_params=False):
 156     """
 157
 158     min_slope_ratio : float
 159         Minimum ratio between the non-contact "left" slope and the
 160         contact "right" slope.
 161     kink_window : int (in bits) or None
 162         Raise `EdgeKink` if the kink is within `kink_window` of the
 163         minimum or maximum `z` position during the approach.  If
 164         `None`, a default value of 2% of the approach range is used.
 165     """
 166     # ususes ddict["approach"] for analysis
 167     # the others are just along to be plotted
 168     _LOG.debug('snalyze surface position data')
 169
 170     data = ddict['approach']
 171     # analyze data, using bilinear model
 172     # y = p0 + p1 x                for x <= p2
 173     #   = p0 + p1 p2 + p3 (x-p2)   for x >= p2
 174     dump_before_index = 0 # 25 # HACK!!
 175     # Generate a reasonable guess...
 176     start_pos = int(data['z'].min())
 177     final_pos = int(data['z'].max())
 178     start_def = int(data['deflection'].min())
 179     final_def = int(data['deflection'].max())
 180     # start_def and start_pos are probably for 2 different points
 181     _LOG.debug('min deflection %d, max deflection %d'
 182                % (start_def, final_def))
 183     _LOG.debug('min position %d, max position %d'
 184                % (start_pos, final_pos))
 185
 186     left_offset   = start_def
 187     left_slope    = 0
 188     kink_position = (final_pos+start_pos)/2.0
 189     right_slope   = 2.0*(final_def-start_def)/(final_pos-start_pos)
 190     pstart = [left_offset, left_slope, kink_position, right_slope]
 191     _LOG.debug('guessed params: %s' % pstart)
 192
 193     offset_scale = (final_pos - start_pos)/100
 194     left_slope_scale = right_slope/10
 195     kink_scale = (final_pos-start_pos)/100
 196     right_slope_scale = right_slope
 197     scale = [offset_scale, left_slope_scale, kink_scale, right_slope_scale]
 198     _LOG.debug('guessed scale: %s' % scale)
 199
 200     def residual(p, y, x):
 201         Y = bilinear(x, p)
 202         return Y-y
 203     params,cov,info,mesg,ier = _leastsq(
 204         residual, pstart,
 205         args=(data["deflection"][dump_before_index:],
 206               data["z"][dump_before_index:]),
 207         full_output=True, maxfev=10000)
 208     _LOG.debug('best fit parameters: %s' % (params,))
 209
 210     if _base_config['matplotlib']:
 211         if not _matplotlib:
 212             raise _matplotlib_import_error
 213         figure = _matplotlib_pyplot.figure()
 214         axes = figure.add_subplot(1, 1, 1)
 215         axes.hold(True)
 216         timestamp = _time.strftime('%H%M%S')
 217         axes.set_title('surf_pos %5g %5g %5g %5g' % tuple(params))
 218         def plot_dict(d, label):
 219             _pylab.plot(d["z"], d["deflection"], '.',label=label)
 220         for n,name in [('ret1', 'first retract'),
 221                        ('mtpod', 'single step in'),
 222                        ('ret2', 'second retract'),
 223                        ('approach', 'main approach'),
 224                        ('ret3', 'return to start')]:
 225             axes.plot(ddict[n]['z'], ddict[n]['deflection'], label=name)
 226         def fit_fn(x, params):
 227             if x <= params[2]:
 228                 return params[0] + params[1]*x
 229             else:
 230                 return (params[0] + params[1]*params[2]
 231                         + params[3]*(x-params[2]))
 232         axes.plot([start_pos, params[2], final_pos],
 233                   [fit_fn(start_pos, params), fit_fn(params[2], params),
 234                    fit_fn(final_pos, params)], '-',label='fit')
 235         #_pylab.legend(loc='best')
 236         figure.show()
 237
 238     # check that the fit is reasonable
 239     # params[1] is slope in non-contact region
 240     # params[2] is kink position
 241     # params[3] is slope in contact region
 242     if kink_window is None:
 243         kink_window = int(0.02*(final_pos-start_pos))
 244
 245     if abs(params[1]*min_slope_ratio) > abs(params[3]):
 246         raise FlatFit(left_slope=params[1], right_slope=params[3])
 247     if params[2] < start_pos+kink_window:
 248         raise EdgeKink(kink=params[2], edge=start_pos, window=kink_window)
 249     if params[2] > final_pos-kink_window:
 250         raise EdgeKink(kink=params[2], edge=final_pos, window=kink_window)
 251     _LOG.debug('surface position %s' % params[2])
 252     if return_all_parameters:
 253         return params
 254     return params[2]
 255
 256 def get_surface_position(piezo, axis_name, max_deflection, **kwargs):
 257     ddict = get_surface_position_data(piezo, axis_name, max_deflection)
 258     return analyze_surface_position_data(ddict, **kwargs)