pypiezo/afm.py

   1 # Copyright (C) 2011-2012 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 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
   8 # version.
   9 #
  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.
  13 #
  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/>.
  16
  17 "Control of a piezo-based atomic force microscope."
  18
  19 import time as _time
  20
  21 import numpy as _numpy
  22
  23 try:
  24     import matplotlib as _matplotlib
  25     import matplotlib.pyplot as _matplotlib_pyplot
  26 except (ImportError, RuntimeError), e:
  27     _matplotlib = None
  28     _matplotlib_import_error = e
  29
  30 try:
  31     import h5py as _h5py
  32     import h5config as _h5config
  33     from h5config.storage.hdf5 import HDF5_Storage as _HDF5_Storage
  34     from h5config.storage.hdf5 import h5_create_group as _h5_create_group
  35 except ImportError, e:
  36     _h5py = None
  37     _h5py_import_error = e
  38
  39 from curses_check_for_keypress import CheckForKeypress as _CheckForKeypress
  40
  41 from . import LOG as _LOG
  42 from . import base as _base
  43 from . import package_config as _package_config
  44 from . import surface as _surface
  45
  46
  47 class AFMPiezo (_base.Piezo):
  48     """A piezo-controlled atomic force microscope.
  49
  50     This particular class expects a single input channel for measuring
  51     deflection.  Other subclasses provide support for multi-segment
  52     deflection measurements.
  53
  54     >>> from pprint import pprint
  55     >>> from pycomedi.device import Device
  56     >>> from pycomedi.subdevice import StreamingSubdevice
  57     >>> from pycomedi.channel import AnalogChannel
  58     >>> from pycomedi.constant import AREF, SUBDEVICE_TYPE, UNIT
  59     >>> from . import config
  60     >>> from . import surface
  61
  62     >>> d = Device('/dev/comedi0')
  63     >>> d.open()
  64
  65     >>> s_in = d.find_subdevice_by_type(SUBDEVICE_TYPE.ai,
  66     ...     factory=StreamingSubdevice)
  67     >>> s_out = d.find_subdevice_by_type(SUBDEVICE_TYPE.ao,
  68     ...     factory=StreamingSubdevice)
  69
  70     >>> axis_channel = s_out.channel(
  71     ...     0, factory=AnalogChannel, aref=AREF.ground)
  72     >>> input_channel = s_in.channel(0, factory=AnalogChannel, aref=AREF.diff)
  73     >>> for chan in [axis_channel, input_channel]:
  74     ...     chan.range = chan.find_range(unit=UNIT.volt, min=-10, max=10)
  75
  76     We set the minimum voltage for the `z` axis to -9 (a volt above
  77     the minimum possible voltage) to help with testing
  78     `.get_surface_position`.  Without this minimum voltage, small
  79     calibration errors could lead to a railed -10 V input for the
  80     first few surface approaching steps, which could lead to an
  81     `EdgeKink` error instead of a `FlatFit` error.
  82
  83     >>> axis_config = config.AxisConfig()
  84     >>> axis_config.update(
  85     ...     {'gain':20, 'sensitivity':8e-9, 'minimum':-9})
  86     >>> axis_config['channel'] = config.OutputChannelConfig()
  87     >>> axis_config['channel']['name'] = 'z'
  88     >>> input_config = config.InputChannelConfig()
  89     >>> input_config['name'] = 'deflection'
  90
  91     >>> a = _base.PiezoAxis(config=axis_config, axis_channel=axis_channel)
  92     >>> a.setup_config()
  93
  94     >>> c = _base.InputChannel(config=input_config, channel=input_channel)
  95     >>> c.setup_config()
  96
  97     >>> p = AFMPiezo(axes=[a], inputs=[c], name='Molly')
  98
  99     >>> deflection = p.read_deflection()
 100     >>> deflection  # doctest: +SKIP
 101     34494L
 102     >>> p.deflection_dtype()
 103     <type 'numpy.uint16'>
 104
 105     We need to know where we are before we can move somewhere
 106     smoothly.
 107
 108     >>> pos = _base.convert_volts_to_bits(p.config.select_config(
 109     ...     'axes', 'z', get_attribute=_base.get_axis_name)['channel'], 0)
 110     >>> p.jump('z', pos)
 111
 112     Usually `.move_to_pos_or_def` is used to approach the surface, but
 113     for testing we assume the z output channel is connected directly
 114     into the deflection input channel.
 115
 116     >>> target_pos = _base.convert_volts_to_bits(
 117     ...     p.config.select_config('axes', 'z',
 118     ...     get_attribute=_base.get_axis_name)['channel'], 2)
 119     >>> step = int((target_pos - pos)/5)
 120     >>> target_def = _base.convert_volts_to_bits(
 121     ...     p.config.select_config('inputs', 'deflection'), 3)
 122     >>> data = p.move_to_pos_or_def('z', target_pos, target_def, step=step,
 123     ...     return_data=True)
 124     >>> p.last_output == {'z': int(target_pos)}
 125     True
 126     >>> pprint(data)  # doctest: +SKIP
 127     {'deflection':
 128        array([32655, 33967, 35280, 36593, 37905, 39218, 39222], dtype=uint16),
 129      'z':
 130        array([32767, 34077, 35387, 36697, 38007, 39317, 39321], dtype=uint16)}
 131
 132     That was a working position-limited approach.  Now move back to
 133     the center and try a deflection-limited approach.
 134
 135     >>> p.jump('z', pos)
 136     >>> target_def = _base.convert_volts_to_bits(
 137     ...     p.config.select_config('inputs', 'deflection'), 1)
 138     >>> data = p.move_to_pos_or_def('z', target_pos, target_def, step=step,
 139     ...     return_data=True)
 140     >>> print (p.last_output['z'] < int(target_pos))
 141     True
 142     >>> pprint(data)  # doctest: +SKIP
 143     {'deflection': array([32655, 33968, 35281, 36593], dtype=uint16),
 144      'z': array([32767, 34077, 35387, 36697], dtype=uint16)}
 145
 146     >>> wiggle_config = config.WiggleConfig()
 147     >>> wiggle_config['offset'] = p.last_output['z']
 148     >>> wiggle_config['wavelength'] = 650e-9
 149     >>> p.wiggle_for_interference(config=wiggle_config,
 150     ...     keypress_test_mode=True)
 151     Press any key to continue
 152
 153     >>> try:
 154     ...     p.get_surface_position('z', max_deflection=target_def)
 155     ... except surface.FlatFit, e:
 156     ...     print 'got FlatFit'
 157     got FlatFit
 158     >>> print e  # doctest: +SKIP
 159     slopes not sufficiently different: 1.0021 and 1.0021
 160     >>> e.right_slope
 161     >>> abs(e.right_slope-1) < 0.1
 162     True
 163     >>> e.left_slope
 164     >>> abs(e.left_slope-1) < 0.1
 165     True
 166
 167     >>> d.close()
 168     """
 169     def _deflection_channel(self):
 170         return self.channel_by_name(name='deflection', direction='input')
 171
 172     def read_deflection(self):
 173         """Return sensor deflection in bits.
 174
 175         TODO: explain how bit <-> volt conversion will work for this
 176         "virtual" channel.
 177         """
 178         return self._deflection_channel().data_read()
 179
 180     def deflection_dtype(self):
 181         "Return a Numpy dtype suitable for deflection bit values."
 182         return self._deflection_channel().subdevice.get_dtype()
 183
 184     def move_to_pos_or_def(self, axis_name, position=None, deflection=None,
 185                            step=1, return_data=False, pre_move_steps=0,
 186                            frequency=None):
 187         """TODO
 188
 189         pre_move_steps : int
 190             number of 'null' steps to take before moving (confirming a
 191             stable input deflection).
 192         frequency : float
 193             The target step frequency in hertz.  If `Null`, go as fast
 194             as possible.  Note that this is software timing, so it
 195             should not be relied upon for precise results.
 196         """
 197         if position is None and deflection is None:
 198             raise ValueError('must specify position, deflection, or both')
 199
 200         if return_data or _package_config['matplotlib']:
 201             aquire_data = True
 202         else:
 203             aquire_data = False
 204
 205         if position is None:
 206             # default to the extreme value in the step direction
 207             if step > 0:
 208                 axis = self.axis_by_name(axis_name)
 209                 position = axis.axis_channel.get_maxdata()
 210             else:
 211                 position = 0
 212         elif deflection is None:
 213             # default to the extreme value
 214             channel = self._deflection_channel(self)
 215             deflection = channel.get_maxdata()
 216
 217         if step == 0:
 218             raise ValueError('must have non-zero step size')
 219         elif step < 0 and position > self.last_output[axis_name]:
 220             step = -step
 221         elif step > 0 and position < self.last_output[axis_name]:
 222             step = -step
 223
 224         log_string = (
 225             'move to position %d or deflection %g on axis %s in steps of %d'
 226             % (position, deflection, axis_name, step))
 227         _LOG.debug(log_string)
 228         current_deflection = self.read_deflection()
 229         log_string = 'current position %d and deflection %g' % (
 230             self.last_output[axis_name], current_deflection)
 231         _LOG.debug(log_string)
 232
 233         if aquire_data:
 234             def_array=[current_deflection]
 235             pos_array=[self.last_output[axis_name]]
 236         for i in range(pre_move_steps):
 237             self.jump(axis_name, piezo.last_output[axis_name])
 238             delection = self.read_deflection()
 239             if aquire_data:
 240                 def_array.append(current_deflection)
 241                 pos_array.append(self.last_output[axis_name])
 242         if frequency is not None:
 243             time_step = 1./frequency
 244             next_time = _time.time() + time_step
 245         # step in until we hit our target position or exceed our target deflection
 246         while (self.last_output[axis_name] != position and
 247                current_deflection < deflection):
 248             dist_to = position - self.last_output[axis_name]
 249             if abs(dist_to) < abs(step):
 250                 jump_to = position
 251             else:
 252                 jump_to = self.last_output[axis_name] + step
 253             self.jump(axis_name, jump_to)
 254             current_deflection = self.read_deflection()
 255             log_string = (
 256                 'current z piezo position %6d, current deflection %6d'
 257                 % (current_deflection, self.last_output[axis_name]))
 258             _LOG.debug(log_string)
 259             if aquire_data:
 260                 def_array.append(current_deflection)
 261                 pos_array.append(self.last_output[axis_name])
 262             if frequency is not None:
 263                 now = _time.time()
 264                 if now < next_time:
 265                     _time.sleep(next_time - now)
 266                 next_time += time_step
 267
 268         log_string = (
 269             'move to position %d or deflection %g on axis %s complete'
 270             % (position, deflection, axis_name))
 271         _LOG.debug(log_string)
 272         log_string = 'current position %d and deflection %g' % (
 273             self.last_output[axis_name], current_deflection)
 274         _LOG.debug(log_string)
 275         if _package_config['matplotlib']:
 276             if not _matplotlib:
 277                 raise _matplotlib_import_error
 278             figure = _matplotlib_pyplot.figure()
 279             axes = figure.add_subplot(1, 1, 1)
 280             axes.hold(True)
 281             timestamp = _time.strftime('%H%M%S')
 282             axes.set_title('step approach %s' % timestamp)
 283             axes.plot(pos_array, def_array, '.', label=timestamp)
 284             #_pylab.legend(loc='best')
 285             figure.show()
 286
 287         if return_data:
 288             data = {
 289                 axis_name:_numpy.array(
 290                     pos_array, dtype=self.channel_dtype(
 291                         axis_name, direction='output')),
 292                 'deflection':_numpy.array(
 293                     def_array, dtype=self.deflection_dtype()),
 294                 }
 295             return data
 296
 297     def wiggle_for_interference(
 298         self, config, plot=True, filename=None, group='/',
 299         keypress_test_mode=False):
 300         """Output a sine wave and measure interference.
 301
 302         With a poorly focused or aligned laser, leaked laser light
 303         reflecting off the surface may interfere with the light
 304         reflected off the cantilever, causing distance-dependent
 305         interference with a period roughly half the laser's
 306         wavelength.  This method wiggles the cantilever near the
 307         surface and monitors the magnitude of deflection oscillation,
 308         allowing the operator to adjust the laser alignment in order
 309         to minimize the interference.
 310
 311         Modern commercial AFMs with computer-aligned lasers must do
 312         something like this automatically.
 313         """
 314         if _package_config['matplotlib']:
 315             plot = True
 316         if config['wavelength'] and config['amplitude']:
 317             log_string = \
 318                 'use either laser_wavelength or amplitude, but not both'
 319             _LOG.warn(log_string)
 320
 321         if None in (config['amplitude'], config['offset']):
 322             output_axis = self.axis_by_name(config['axis'])
 323             maxdata = output_axis.axis_channel.get_maxdata()
 324             midpoint = int(maxdata/2)
 325             if config['offset'] is None:
 326                 offset = midpoint
 327                 _LOG.debug(('generated offset for interference wiggle: {}'
 328                             ).format(config['offset']))
 329             if config['amplitude'] is None:
 330                 if config['offset'] <= midpoint:
 331                     max_amplitude = int(config['offset'])
 332                 else:
 333                     max_amplitude = int(maxdata-config['offset'])
 334                 offset_meters = _base.convert_bits_to_meters(
 335                     output_axis.config, config['offset'])
 336                 if config['wavelength'] is None:
 337                     config['amplitude'] = 0.5*max_amplitude
 338                 else:
 339                     bit_wavelength = _base.convert_meters_to_bits(
 340                         output_axis.config,
 341                         offset_meters + config['wavelength']
 342                         ) - config['offset']
 343                     config['amplitude'] = 2*bit_wavelength
 344                 _LOG.debug(('generated amplitude for interference wiggle: {}'
 345                             ).format(config['amplitude']))
 346                 if config['amplitude'] > max_amplitude:
 347                     raise ValueError(
 348                         'no room for a two wavelength wiggle ({} > {})'.format(
 349                             config['amplitude'], max_amplitude))
 350
 351         n = config['samples']  # samples in a single oscillation
 352         scan_frequency = config['frequency'] * n
 353         out = (config['amplitude']
 354                * _numpy.sin(_numpy.arange(2*n)*2*_numpy.pi/n)
 355                + config['offset'])
 356         # 2*n for 2 periods, so you can judge precision
 357         out = out.reshape((len(out), 1)).astype(
 358             self.channel_dtype(config['axis'], direction='output'))
 359
 360         _LOG.debug('oscillate for interference wiggle ({})'.format(config))
 361
 362         if filename:
 363             if not _h5py:
 364                 raise _h5py_import_error
 365             if not output_axis:  # from amplitude/offset setup
 366                 output_axis = afm.piezo.axis_by_name(config['axis'])
 367             input_channel = afm.piezo.input_channel_by_name(config['input'])
 368             with _h5py.File(filename, 'w') as f:
 369                 cwg = _h5_create_group(f, group)
 370                 storage = _HDF5_Storage()
 371                 for config_,key in [
 372                     (config, 'config/wiggle'),
 373                     (output.axis.config,
 374                      'config/{}/axis'.format(config['axis'])),
 375                     (input_channel.config,
 376                      'config/{}/channel'.format(config['input']))]:
 377                     if config_ is None:
 378                         continue
 379                     config_cwg = _h5_create_group(cwg, key)
 380                     storage.save(config=config, group=config_cwg)
 381         if plot:
 382             if not _matplotlib:
 383                 raise _matplotlib_import_error
 384             figure = _matplotlib_pyplot.figure()
 385             axes = figure.add_subplot(1, 1, 1)
 386             axes.hold(False)
 387             timestamp = _time.strftime('%H%M%S')
 388             axes.set_title('wiggle for interference %s' % timestamp)
 389             plot_p = axes.plot(out, out, 'b.-')
 390             figure.show()
 391         cycle = 0
 392         c = _CheckForKeypress(test_mode=keypress_test_mode)
 393         while c.input() == None:
 394             # input will need processing for multi-segment AFMs...
 395             data = self.ramp(
 396                 out, scan_frequency, output_names=[config['axis']],
 397                 input_names=[config['input']])
 398             _LOG.debug('completed a wiggle round')
 399             if filename:
 400                 timestamp = ('{0}-{1:02d}-{2:02d}T{3:02d}-{4:02d}-{5:02d}'
 401                              ).format(*_time.localtime())
 402                 with _h5py.File(filename, 'a') as f:
 403                     wiggle_group = _h5_create_group(f, group)
 404                     cwg = _h5_create_group(
 405                         wiggle_group, 'wiggle/{}'.format(cycle))
 406                     cwg['time'] = timestamp
 407                     cwg['raw/{}'.format(config['axis'])] = out
 408                     cwg['raw/{}'.format(config['input'])] = data
 409             if plot:
 410                 plot_p[0].set_ydata(data[:,0])
 411                 axes.set_ylim([data.min(), data.max()])
 412                 #_flush_plot()
 413             cycle += 1
 414         self.last_output[config['axis']] = out[-1,0]
 415         _LOG.debug('interference wiggle complete')
 416
 417     get_surface_position = _surface.get_surface_position
 418
 419
 420 #def ramp
 421 #        if USE_ABCD_DEFLECTION :
 422 #            for i in range(4) : # i is the photodiode element (0->A, 1->B, ...)
 423 #                self.curIn[i] = out["Deflection segment"][i][-1]
 424 #        else :
 425 #            self.curIn[self.chan_info.def_ind] = out["deflection"][-1]
 426
 427
 428 #class FourSegmentAFM (AFM):
 429 #    def read_deflection(self):
 430 #        "Return sensor deflection in bits."
 431 #        A = int(self.curIn[self.chan_info.def_ind[0]])
 432 #        B = int(self.curIn[self.chan_info.def_ind[1]])
 433 #        C = int(self.curIn[self.chan_info.def_ind[2]])
 434 #        D = int(self.curIn[self.chan_info.def_ind[3]])
 435 #        df = float((A+B)-(C+D))/(A+B+C+D)
 436 #        dfout = int(df * 2**15) + 2**15
 437 #        if TEXT_VERBOSE :
 438 #            print "Current deflection %d (%d, %d, %d, %d)" \
 439 #                % (dfout, A, B, C, D)
 440 #        return dfout
 441
 442
 443 #def test_smoothness(zp, plotVerbose=True):
 444 #    posA = 20000
 445 #    posB = 50000
 446 #    setpoint = zp.def_V2in(3)
 447 #    steps = 200
 448 #    outfreq = 1e5
 449 #    outarray = linspace(posB, posA, 1000)
 450 #    indata=[]
 451 #    outdata=[]
 452 #    curVals = zp.jumpToPos(posA)
 453 #    zp.pCurVals(curVals)
 454 #    _sleep(1) # let jitters die down
 455 #    for i in range(10):
 456 #        print "ramp %d to %d" % (zp.curPos(), posB)
 457 #        curVals, data = moveToPosOrDef(zp, posB, setpoint, step=steps,
 458 #                                       return_data = True)
 459 #        indata.append(data)
 460 #        out = zp.ramp(outarray, outfreq)
 461 #        outdata.append(out)
 462 #    if plotVerbose:
 463 #        from pylab import figure, plot, title, legend, hold, subplot
 464 #    if PYLAB_VERBOSE or plotVerbose:
 465 #        _import_pylab()
 466 #        _pylab.figure(BASE_FIG_NUM+4)
 467 #        for i in range(10):
 468 #            _pylab.plot(indata[i]['z'],
 469 #                        indata[i]['deflection'], '+--', label='in')
 470 #            _pylab.plot(outdata[i]['z'],
 471 #                        outdata[i]['deflection'], '.-', label='out')
 472 #        _pylab.title('test smoothness (step in, ramp out)')
 473 #        #_pylab.legend(loc='best')
 474 #
 475 #def test():
 476 #    import z_piezo
 477 #    zp = z_piezo.z_piezo()
 478 #    curVals = zp.moveToPosOrDef(zp.pos_nm2out(600), defl=zp.curDef()+6000, step=(zp.pos_nm2out(10)-zp.pos_nm2out(0)))
 479 #    if TEXT_VERBOSE:
 480 #        zp.pCurVals(curVals)
 481 #    pos = zp.getSurfPos(maxDefl=zp.curDef()+6000)
 482 #    if TEXT_VERBOSE:
 483 #        print "Surface at %g nm", pos
 484 #    print "success"
 485 #    if PYLAB_VERBOSE and _final_flush_plot != None:
 486 #        _final_flush_plot()
 487