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, RuntimeError), e:
  34     _matplotlib = None
  35     _matplotlib_import_error = e
  36
  37 from . import LOG as _LOG
  38 from . import package_config as _package_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
  64 class EdgeKink (PoorFit):
  65     def __init__(self, kink, edge, window):
  66         self.kink = kink
  67         self.edge = edge
  68         self.window = window
  69         msg = 'no kink (kink %d not within %d of %d)' % (kink, window, edge)
  70         super(EdgeKink, self).__init__(self, msg)
  71
  72
  73 def _linspace(*args, **kwargs):
  74     dtype = kwargs.pop('dtype')
  75     out = _numpy.linspace(*args, **kwargs)
  76     return out.reshape((len(out), 1)).astype(dtype)
  77
  78 def _get_min_max_positions(piezo, axis_name, min_position=None,
  79                            max_position=None):
  80     output_axis = piezo.axis_by_name(axis_name)
  81     if min_position is None:
  82         min_position = _base.convert_volts_to_bits(
  83             output_axis.config['channel'],
  84             output_axis.config['minimum'])
  85     if max_position is None:
  86         max_position = _base.convert_volts_to_bits(
  87             output_axis.config['channel'],
  88             output_axis.config['maximum'])
  89     return (min_position, max_position)
  90
  91 def get_surface_position_data(piezo, axis_name, max_deflection, steps=2000,
  92                               frequency=10e3, min_position=None,
  93                               max_position=None):
  94     "Measure the distance to the surface"
  95     _LOG.debug('get surface position')
  96     orig_position = piezo.last_output[axis_name]
  97     # fully retract the piezo
  98     min_position,max_position = _get_min_max_positions(
  99         piezo, axis_name, min_position=min_position, max_position=max_position)
 100     _LOG.debug('retract the piezo to %d' % min_position)
 101     dtype = piezo.channel_dtype(axis_name, direction='output')
 102     out = _linspace(orig_position, min_position, steps, dtype=dtype)
 103     out = out.reshape((len(out), 1)).astype(
 104         piezo.channel_dtype(axis_name, direction='output'))
 105     ramp_kwargs = {
 106         'frequency': frequency,
 107         'output_names': [axis_name],
 108         'input_names': ['deflection'],
 109         }
 110     ret1 = piezo.named_ramp(data=out, **ramp_kwargs)
 111     # locate high deflection position
 112     _LOG.debug('approach until there is dangerous deflection (> %d)'
 113                % max_deflection)
 114     if SLEEP_DURING_SURF_POS_AQUISITION == True:
 115         _sleep(.2) # sleeping briefly seems to reduce bounce?
 116     mtpod = piezo.move_to_pos_or_def(
 117         axis_name=axis_name, position=max_position, deflection=max_deflection,
 118         step=(max_position-min_position)/steps, return_data=True)
 119     high_contact_pos = piezo.last_output[axis_name]
 120     # fully retract the piezo again
 121     _LOG.debug('retract the piezo to %d again' % min_position)
 122     if SLEEP_DURING_SURF_POS_AQUISITION == True:
 123         _sleep(.2)
 124     out = _linspace(high_contact_pos, min_position, steps, dtype=dtype)
 125     ret2 = piezo.named_ramp(data=out, **ramp_kwargs)
 126     # scan to the high contact position
 127     _LOG.debug('ramp in to the deflected position %d' % high_contact_pos)
 128     if SLEEP_DURING_SURF_POS_AQUISITION == True:
 129         _sleep(.2)
 130     out = _linspace(min_position, high_contact_pos, steps, dtype=dtype)
 131     data = piezo.named_ramp(data=out, **ramp_kwargs)
 132     if SLEEP_DURING_SURF_POS_AQUISITION == True:
 133         _sleep(.2)
 134     # return to the original position
 135     _LOG.debug('return to the original position %d' % orig_position)
 136     out = _linspace(high_contact_pos, orig_position, steps, dtype=dtype)
 137     ret3 = piezo.named_ramp(data=out, **ramp_kwargs)
 138     return {'ret1':ret1, 'mtpod':mtpod, 'ret2':ret2,
 139             'approach':data, 'ret3':ret3}
 140
 141 def bilinear(x, params):
 142     """bilinear fit for surface bumps.  Model has two linear regimes
 143     which meet at x=kink_position and have independent slopes.
 144
 145     `x` should be a `numpy.ndarray`.
 146     """
 147     left_offset,left_slope,kink_position,right_slope = params
 148     left_mask = x < kink_position
 149     right_mask = x >= kink_position # = not left_mask
 150     left_y = left_offset + left_slope*x
 151     right_y = (left_offset + left_slope*kink_position
 152                + right_slope*(x-kink_position))
 153     return left_mask * left_y + right_mask * right_y
 154
 155 def analyze_surface_position_data(
 156     ddict, min_slope_ratio=10.0, kink_window=None,
 157     return_all_parameters=False):
 158     """
 159
 160     min_slope_ratio : float
 161         Minimum ratio between the non-contact "left" slope and the
 162         contact "right" slope.
 163     kink_window : int (in bits) or None
 164         Raise `EdgeKink` if the kink is within `kink_window` of the
 165         minimum or maximum `z` position during the approach.  If
 166         `None`, a default value of 2% of the approach range is used.
 167     """
 168     # ususes ddict["approach"] for analysis
 169     # the others are just along to be plotted
 170     _LOG.debug('snalyze surface position data')
 171
 172     data = ddict['approach']
 173     # analyze data, using bilinear model
 174     # y = p0 + p1 x                for x <= p2
 175     #   = p0 + p1 p2 + p3 (x-p2)   for x >= p2
 176     dump_before_index = 0 # 25 # HACK!!
 177     # Generate a reasonable guess...
 178     start_pos = int(data['z'].min())
 179     final_pos = int(data['z'].max())
 180     start_def = int(data['deflection'].min())
 181     final_def = int(data['deflection'].max())
 182     # start_def and start_pos are probably for 2 different points
 183     _LOG.debug('min deflection %d, max deflection %d'
 184                % (start_def, final_def))
 185     _LOG.debug('min position %d, max position %d'
 186                % (start_pos, final_pos))
 187
 188     left_offset   = start_def
 189     left_slope    = 0
 190     kink_position = (final_pos+start_pos)/2.0
 191     right_slope   = 2.0*(final_def-start_def)/(final_pos-start_pos)
 192     pstart = [left_offset, left_slope, kink_position, right_slope]
 193     _LOG.debug('guessed params: %s' % pstart)
 194
 195     offset_scale = (final_pos - start_pos)/100
 196     left_slope_scale = right_slope/10
 197     kink_scale = (final_pos-start_pos)/100
 198     right_slope_scale = right_slope
 199     scale = [offset_scale, left_slope_scale, kink_scale, right_slope_scale]
 200     _LOG.debug('guessed scale: %s' % scale)
 201
 202     def residual(p, y, x):
 203         Y = bilinear(x, p)
 204         return Y-y
 205     params,cov,info,mesg,ier = _leastsq(
 206         residual, pstart,
 207         args=(data["deflection"][dump_before_index:],
 208               data["z"][dump_before_index:]),
 209         full_output=True, maxfev=10000)
 210     _LOG.debug('best fit parameters: %s' % (params,))
 211
 212     if _package_config['matplotlib']:
 213         if not _matplotlib:
 214             raise _matplotlib_import_error
 215         figure = _matplotlib_pyplot.figure()
 216         axes = figure.add_subplot(1, 1, 1)
 217         axes.hold(True)
 218         timestamp = _time.strftime('%H%M%S')
 219         axes.set_title('surf_pos %5g %5g %5g %5g' % tuple(params))
 220         def plot_dict(d, label):
 221             _pylab.plot(d["z"], d["deflection"], '.',label=label)
 222         for n,name in [('ret1', 'first retract'),
 223                        ('mtpod', 'single step in'),
 224                        ('ret2', 'second retract'),
 225                        ('approach', 'main approach'),
 226                        ('ret3', 'return to start')]:
 227             axes.plot(ddict[n]['z'], ddict[n]['deflection'], label=name)
 228         def fit_fn(x, params):
 229             if x <= params[2]:
 230                 return params[0] + params[1]*x
 231             else:
 232                 return (params[0] + params[1]*params[2]
 233                         + params[3]*(x-params[2]))
 234         axes.plot([start_pos, params[2], final_pos],
 235                   [fit_fn(start_pos, params), fit_fn(params[2], params),
 236                    fit_fn(final_pos, params)], '-',label='fit')
 237         #_pylab.legend(loc='best')
 238         figure.show()
 239
 240     # check that the fit is reasonable
 241     # params[1] is slope in non-contact region
 242     # params[2] is kink position
 243     # params[3] is slope in contact region
 244     if kink_window is None:
 245         kink_window = int(0.02*(final_pos-start_pos))
 246
 247     if abs(params[1]*min_slope_ratio) > abs(params[3]):
 248         raise FlatFit(left_slope=params[1], right_slope=params[3])
 249     if params[2] < start_pos+kink_window:
 250         raise EdgeKink(kink=params[2], edge=start_pos, window=kink_window)
 251     if params[2] > final_pos-kink_window:
 252         raise EdgeKink(kink=params[2], edge=final_pos, window=kink_window)
 253     _LOG.debug('surface position %s' % params[2])
 254     if return_all_parameters:
 255         return params
 256     return params[2]
 257
 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)