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     >>> abs(e.right_slope-1) < 0.1
 161     True
 162     >>> abs(e.left_slope-1) < 0.1
 163     True
 164
 165     >>> d.close()
 166     """
 167     def _deflection_channel(self):
 168         return self.channel_by_name(name='deflection', direction='input')
 169
 170     def read_deflection(self):
 171         """Return sensor deflection in bits.
 172
 173         TODO: explain how bit <-> volt conversion will work for this
 174         "virtual" channel.
 175         """
 176         return self._deflection_channel().data_read()
 177
 178     def deflection_dtype(self):
 179         "Return a Numpy dtype suitable for deflection bit values."
 180         return self._deflection_channel().subdevice.get_dtype()
 181
 182     def move_to_pos_or_def(self, axis_name, position=None, deflection=None,
 183                            step=1, return_data=False, pre_move_steps=0,
 184                            frequency=None):
 185         """TODO
 186
 187         pre_move_steps : int
 188             number of 'null' steps to take before moving (confirming a
 189             stable input deflection).
 190         frequency : float
 191             The target step frequency in hertz.  If `Null`, go as fast
 192             as possible.  Note that this is software timing, so it
 193             should not be relied upon for precise results.
 194         """
 195         if position is None and deflection is None:
 196             raise ValueError('must specify position, deflection, or both')
 197
 198         if return_data or _package_config['matplotlib']:
 199             aquire_data = True
 200         else:
 201             aquire_data = False
 202
 203         if position is None:
 204             # default to the extreme value in the step direction
 205             if step > 0:
 206                 axis = self.axis_by_name(axis_name)
 207                 position = axis.axis_channel.get_maxdata()
 208             else:
 209                 position = 0
 210         elif deflection is None:
 211             # default to the extreme value
 212             channel = self._deflection_channel(self)
 213             deflection = channel.get_maxdata()
 214
 215         if step == 0:
 216             raise ValueError('must have non-zero step size')
 217         elif step < 0 and position > self.last_output[axis_name]:
 218             step = -step
 219         elif step > 0 and position < self.last_output[axis_name]:
 220             step = -step
 221
 222         log_string = (
 223             'move to position %d or deflection %g on axis %s in steps of %d'
 224             % (position, deflection, axis_name, step))
 225         _LOG.debug(log_string)
 226         current_deflection = self.read_deflection()
 227         log_string = 'current position %d and deflection %g' % (
 228             self.last_output[axis_name], current_deflection)
 229         _LOG.debug(log_string)
 230
 231         if aquire_data:
 232             def_array=[current_deflection]
 233             pos_array=[self.last_output[axis_name]]
 234         for i in range(pre_move_steps):
 235             self.jump(axis_name, piezo.last_output[axis_name])
 236             delection = self.read_deflection()
 237             if aquire_data:
 238                 def_array.append(current_deflection)
 239                 pos_array.append(self.last_output[axis_name])
 240         if frequency is not None:
 241             time_step = 1./frequency
 242             next_time = _time.time() + time_step
 243         # step in until we hit our target position or exceed our target deflection
 244         while (self.last_output[axis_name] != position and
 245                current_deflection < deflection):
 246             dist_to = position - self.last_output[axis_name]
 247             if abs(dist_to) < abs(step):
 248                 jump_to = position
 249             else:
 250                 jump_to = self.last_output[axis_name] + step
 251             self.jump(axis_name, jump_to)
 252             current_deflection = self.read_deflection()
 253             log_string = (
 254                 'current z piezo position %6d, current deflection %6d'
 255                 % (current_deflection, self.last_output[axis_name]))
 256             _LOG.debug(log_string)
 257             if aquire_data:
 258                 def_array.append(current_deflection)
 259                 pos_array.append(self.last_output[axis_name])
 260             if frequency is not None:
 261                 now = _time.time()
 262                 if now < next_time:
 263                     _time.sleep(next_time - now)
 264                 next_time += time_step
 265
 266         log_string = (
 267             'move to position %d or deflection %g on axis %s complete'
 268             % (position, deflection, axis_name))
 269         _LOG.debug(log_string)
 270         log_string = 'current position %d and deflection %g' % (
 271             self.last_output[axis_name], current_deflection)
 272         _LOG.debug(log_string)
 273         if _package_config['matplotlib']:
 274             if not _matplotlib:
 275                 raise _matplotlib_import_error
 276             figure = _matplotlib_pyplot.figure()
 277             axes = figure.add_subplot(1, 1, 1)
 278             axes.hold(True)
 279             timestamp = _time.strftime('%H%M%S')
 280             axes.set_title('step approach %s' % timestamp)
 281             axes.plot(pos_array, def_array, '.', label=timestamp)
 282             #_pylab.legend(loc='best')
 283             figure.show()
 284
 285         if return_data:
 286             data = {
 287                 axis_name:_numpy.array(
 288                     pos_array, dtype=self.channel_dtype(
 289                         axis_name, direction='output')),
 290                 'deflection':_numpy.array(
 291                     def_array, dtype=self.deflection_dtype()),
 292                 }
 293             return data
 294
 295     def wiggle_for_interference(
 296         self, config, plot=True, filename=None, group='/',
 297         keypress_test_mode=False):
 298         """Output a sine wave and measure interference.
 299
 300         With a poorly focused or aligned laser, leaked laser light
 301         reflecting off the surface may interfere with the light
 302         reflected off the cantilever, causing distance-dependent
 303         interference with a period roughly half the laser's
 304         wavelength.  This method wiggles the cantilever near the
 305         surface and monitors the magnitude of deflection oscillation,
 306         allowing the operator to adjust the laser alignment in order
 307         to minimize the interference.
 308
 309         Modern commercial AFMs with computer-aligned lasers must do
 310         something like this automatically.
 311         """
 312         if _package_config['matplotlib']:
 313             plot = True
 314         if config['wavelength'] and config['amplitude']:
 315             log_string = \
 316                 'use either laser_wavelength or amplitude, but not both'
 317             _LOG.warn(log_string)
 318
 319         if None in (config['amplitude'], config['offset']):
 320             output_axis = self.axis_by_name(config['axis'])
 321             maxdata = output_axis.axis_channel.get_maxdata()
 322             midpoint = int(maxdata/2)
 323             if config['offset'] is None:
 324                 offset = midpoint
 325                 _LOG.debug(('generated offset for interference wiggle: {}'
 326                             ).format(config['offset']))
 327             if config['amplitude'] is None:
 328                 if config['offset'] <= midpoint:
 329                     max_amplitude = int(config['offset'])
 330                 else:
 331                     max_amplitude = int(maxdata-config['offset'])
 332                 offset_meters = _base.convert_bits_to_meters(
 333                     output_axis.config, config['offset'])
 334                 if config['wavelength'] is None:
 335                     config['amplitude'] = 0.5*max_amplitude
 336                 else:
 337                     bit_wavelength = _base.convert_meters_to_bits(
 338                         output_axis.config,
 339                         offset_meters + config['wavelength']
 340                         ) - config['offset']
 341                     config['amplitude'] = 2*bit_wavelength
 342                 _LOG.debug(('generated amplitude for interference wiggle: {}'
 343                             ).format(config['amplitude']))
 344                 if config['amplitude'] > max_amplitude:
 345                     raise ValueError(
 346                         'no room for a two wavelength wiggle ({} > {})'.format(
 347                             config['amplitude'], max_amplitude))
 348
 349         n = config['samples']  # samples in a single oscillation
 350         scan_frequency = config['frequency'] * n
 351         out = (config['amplitude']
 352                * _numpy.sin(_numpy.arange(2*n)*2*_numpy.pi/n)
 353                + config['offset'])
 354         # 2*n for 2 periods, so you can judge precision
 355         out = out.reshape((len(out), 1)).astype(
 356             self.channel_dtype(config['axis'], direction='output'))
 357
 358         _LOG.debug('oscillate for interference wiggle ({})'.format(config))
 359
 360         if filename:
 361             if not _h5py:
 362                 raise _h5py_import_error
 363             if not output_axis:  # from amplitude/offset setup
 364                 output_axis = afm.piezo.axis_by_name(config['axis'])
 365             input_channel = afm.piezo.input_channel_by_name(config['input'])
 366             with _h5py.File(filename, 'w') as f:
 367                 cwg = _h5_create_group(f, group)
 368                 storage = _HDF5_Storage()
 369                 for config_,key in [
 370                     (config, 'config/wiggle'),
 371                     (output.axis.config,
 372                      'config/{}/axis'.format(config['axis'])),
 373                     (input_channel.config,
 374                      'config/{}/channel'.format(config['input']))]:
 375                     if config_ is None:
 376                         continue
 377                     config_cwg = _h5_create_group(cwg, key)
 378                     storage.save(config=config, group=config_cwg)
 379         if plot:
 380             if not _matplotlib:
 381                 raise _matplotlib_import_error
 382             interactive = _matplotlib.is_interactive()
 383             _matplotlib.interactive(True)
 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             _matplotlib_pyplot.draw()
 392             _matplotlib_pyplot.show()
 393         cycle = 0
 394         c = _CheckForKeypress(test_mode=keypress_test_mode)
 395         while c.input() == None:
 396             # input will need processing for multi-segment AFMs...
 397             data = self.ramp(
 398                 out, scan_frequency, output_names=[config['axis']],
 399                 input_names=[config['input']])
 400             _LOG.debug('completed a wiggle round')
 401             if filename:
 402                 timestamp = ('{0}-{1:02d}-{2:02d}T{3:02d}-{4:02d}-{5:02d}'
 403                              ).format(*_time.localtime())
 404                 with _h5py.File(filename, 'a') as f:
 405                     wiggle_group = _h5_create_group(f, group)
 406                     cwg = _h5_create_group(
 407                         wiggle_group, 'wiggle/{}'.format(cycle))
 408                     cwg['time'] = timestamp
 409                     cwg['raw/{}'.format(config['axis'])] = out
 410                     cwg['raw/{}'.format(config['input'])] = data
 411             if plot:
 412                 plot_p[0].set_ydata(data[:,0])
 413                 axes.set_ylim([data.min(), data.max()])
 414                 _matplotlib_pyplot.draw()
 415             cycle += 1
 416         if plot:
 417             _matplotlib.interactive(interactive)
 418         self.last_output[config['axis']] = out[-1,0]
 419         _LOG.debug('interference wiggle complete')
 420
 421     get_surface_position = _surface.get_surface_position
 422
 423
 424 #def ramp
 425 #        if USE_ABCD_DEFLECTION :
 426 #            for i in range(4) : # i is the photodiode element (0->A, 1->B, ...)
 427 #                self.curIn[i] = out["Deflection segment"][i][-1]
 428 #        else :
 429 #            self.curIn[self.chan_info.def_ind] = out["deflection"][-1]
 430
 431
 432 #class FourSegmentAFM (AFM):
 433 #    def read_deflection(self):
 434 #        "Return sensor deflection in bits."
 435 #        A = int(self.curIn[self.chan_info.def_ind[0]])
 436 #        B = int(self.curIn[self.chan_info.def_ind[1]])
 437 #        C = int(self.curIn[self.chan_info.def_ind[2]])
 438 #        D = int(self.curIn[self.chan_info.def_ind[3]])
 439 #        df = float((A+B)-(C+D))/(A+B+C+D)
 440 #        dfout = int(df * 2**15) + 2**15
 441 #        if TEXT_VERBOSE :
 442 #            print "Current deflection %d (%d, %d, %d, %d)" \
 443 #                % (dfout, A, B, C, D)
 444 #        return dfout
 445
 446
 447 #def test_smoothness(zp, plotVerbose=True):
 448 #    posA = 20000
 449 #    posB = 50000
 450 #    setpoint = zp.def_V2in(3)
 451 #    steps = 200
 452 #    outfreq = 1e5
 453 #    outarray = linspace(posB, posA, 1000)
 454 #    indata=[]
 455 #    outdata=[]
 456 #    curVals = zp.jumpToPos(posA)
 457 #    zp.pCurVals(curVals)
 458 #    _sleep(1) # let jitters die down
 459 #    for i in range(10):
 460 #        print "ramp %d to %d" % (zp.curPos(), posB)
 461 #        curVals, data = moveToPosOrDef(zp, posB, setpoint, step=steps,
 462 #                                       return_data = True)
 463 #        indata.append(data)
 464 #        out = zp.ramp(outarray, outfreq)
 465 #        outdata.append(out)
 466 #    if plotVerbose:
 467 #        from pylab import figure, plot, title, legend, hold, subplot
 468 #    if PYLAB_VERBOSE or plotVerbose:
 469 #        _import_pylab()
 470 #        _pylab.figure(BASE_FIG_NUM+4)
 471 #        for i in range(10):
 472 #            _pylab.plot(indata[i]['z'],
 473 #                        indata[i]['deflection'], '+--', label='in')
 474 #            _pylab.plot(outdata[i]['z'],
 475 #                        outdata[i]['deflection'], '.-', label='out')
 476 #        _pylab.title('test smoothness (step in, ramp out)')
 477 #        #_pylab.legend(loc='best')
 478 #
 479 #def test():
 480 #    import z_piezo
 481 #    zp = z_piezo.z_piezo()
 482 #    curVals = zp.moveToPosOrDef(zp.pos_nm2out(600), defl=zp.curDef()+6000, step=(zp.pos_nm2out(10)-zp.pos_nm2out(0)))
 483 #    if TEXT_VERBOSE:
 484 #        zp.pCurVals(curVals)
 485 #    pos = zp.getSurfPos(maxDefl=zp.curDef()+6000)
 486 #    if TEXT_VERBOSE:
 487 #        print "Surface at %g nm", pos
 488 #    print "success"
 489 #    if PYLAB_VERBOSE and _final_flush_plot != None:
 490 #        _final_flush_plot()
 491