calibcant/calibrate.py

   1 # calibcant - tools for thermally calibrating AFM cantilevers
   2 #
   3 # Copyright (C) 2008-2011 W. Trevor King <wking@drexel.edu>
   4 #
   5 # This file is part of calibcant.
   6 #
   7 # calibcant is free software: you can redistribute it and/or
   8 # modify it under the terms of the GNU Lesser General Public
   9 # License as published by the Free Software Foundation, either
  10 # version 3 of the License, or (at your option) any later version.
  11 #
  12 # calibcant is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU Lesser General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU Lesser General Public
  18 # License along with calibcant.  If not, see
  19 # <http://www.gnu.org/licenses/>.
  20
  21 """Acquire and analyze cantilever calibration data.
  22
  23 The relevent physical quantities are:
  24
  25 * Vzp_out  Output z-piezo voltage (what we generate)
  26 * Vzp      Applied z-piezo voltage (after external ZPGAIN)
  27 * Zp       The z-piezo position
  28 * Zcant    The cantilever vertical deflection
  29 * Vphoto   The photodiode vertical deflection voltage (what we measure)
  30 * Fcant    The force on the cantilever
  31 * T        The temperature of the cantilever and surrounding solution
  32 *          (another thing we measure or guess)
  33 * k_b      Boltzmann's constant
  34
  35 Which are related by the parameters:
  36
  37 * zpGain           Vzp_out / Vzp
  38 * zpSensitivity    Zp / Vzp
  39 * photoSensitivity Vphoto / Zcant
  40 * k_cant           Fcant / Zcant
  41
  42 Cantilever calibration assumes a pre-calibrated z-piezo
  43 (zpSensitivity) and a amplifier (zpGain).  In our lab, the z-piezo is
  44 calibrated by imaging a calibration sample, which has features with
  45 well defined sizes, and the gain is set with a knob on the Nanoscope.
  46
  47 photoSensitivity is measured by bumping the cantilever against the
  48 surface, where Zp = Zcant (see bump_acquire() and the bump_analyze
  49 submodule).
  50
  51 k_cant is measured by watching the cantilever vibrate in free solution
  52 (see the vib_acquire() and the vib_analyze submodule).  The average
  53 energy of the cantilever in the vertical direction is given by the
  54 equipartition theorem.
  55
  56 .. math::  \frac{1}{2} k_b T = \frac{1}{2} k_cant <Zcant**2>
  57
  58 so
  59
  60 .. math::   k_cant = \frac{k_b T}{Zcant**2}
  61
  62 but
  63
  64 .. math::   Zcant = \frac{Vphoto}{photoSensitivity}
  65
  66 so
  67
  68 .. math:: k_cant = \frac{k_b T * photoSensitivty^2}{<Vphoto**2>}
  69
  70 We measured photoSensitivity with the surface bumps.  We can either
  71 measure T using an external function (see temperature.py), or just
  72 estimate it (see T_acquire() and the T_analyze submodule).  Guessing
  73 room temp ~22 deg C is actually fairly reasonable.  Assuming the
  74 actual fluid temperature is within +/- 5 deg, the error in the spring
  75 constant k_cant is within 5/273.15 ~= 2%.  A time series of Vphoto
  76 while we're far from the surface and not changing Vzp_out will give us
  77 the average variance <Vphoto**2>.
  78
  79 We do all these measurements a few times to estimate statistical
  80 errors.
  81
  82 The functions are layed out in the families::
  83
  84   bump_*(), vib_*(), T_*(), and calib_*()
  85
  86 For each family, * can be any of:
  87
  88 * acquire       get real-world data
  89 * save         store real-world data to disk
  90 * load         get real-world data from disk
  91 * analyze      interperate the real-world data.
  92 * plot         show a nice graphic to convince people we're working :p
  93
  94 A family name without any `_*` extension (e.g. `bump()`), runs
  95 `*_acquire()`, `*_analyze()`, and `*_save()`.  `*_analyze()` will run
  96 `*_plot()` if `matplotlib` is set in `calibcant.base_config`.
  97 """
  98
  99 from numpy import zeros as _zeros
 100 from numpy import float as _float
 101 from time import sleep as _sleep
 102
 103 from . import LOG as _LOG
 104
 105 from .bump import bump as _bump
 106 from .T import T as _T
 107 from .vib import vib as _vib
 108 from .analyze import calib_analyze as _calib_analyze
 109 from .analyze import calib_save as _calib_save
 110
 111
 112 def move_far_from_surface(stepper, distance):
 113     """Step back approximately `distance` meters.
 114     """
 115     steps = int(distance/stepper.step_size)
 116     _LOG.info('step back %d steps (~%g m)' % (steps, distance))
 117     stepper.step_relative(-steps)
 118
 119 def calib_acquire(afm, calibration_config, bump_config, temperature_config,
 120                   vibration_config, filename=None, group='/'):
 121     """Acquire data for calibrating a cantilever in one function.
 122
 123     Inputs:
 124       afm                 a pyafm.AFM instance
 125       calibration_config  a .config._CalibrationConfig instance
 126       bump_config         a .config._BumpConfig instance
 127       temperature_config            a .config._TConfig instance
 128       vibration_config    a .config._VibrationConfig instance
 129
 130     Outputs (all are arrays of recorded data):
 131       bumps measured (V_photodiode / nm_tip) proportionality constant
 132       Ts    measured temperature (K)
 133       vibs  measured V_photodiode variance (Volts**2) in free solution
 134
 135     The temperatures are collected after moving far from the surface
 136     but before and vibrations are measured to give everything time to
 137     settle after the big move.
 138     """
 139     assert group.endswith('/'), group
 140
 141     bumps = _zeros((calibration_config['num-bumps'],), dtype=_float)
 142     for i in range(calibration_config['num-bumps']):
 143         _LOG.info('acquire bump %d of %d' % (i, calibration_config['num-bumps']))
 144         bumps[i] = _bump(afm=afm, bump_config=bump_config,
 145                          filename=filename, group='%sbump/%d/' % (group, i))
 146     _LOG.debug('bumps: %s' % bumps)
 147
 148     move_far_from_surface(
 149         afm.stepper, distance=calibration_config['vibration-spacing'])
 150
 151     Ts = _zeros((calibration_config['num-temperatures'],), dtype=_float)
 152     for i in range(calibration_config['num-temperatures']):
 153         _LOG.info('acquire T %d of %d'
 154                  % (i, calibration_config['num-temperatures']))
 155         Ts[i] = _T(
 156             get_T=afm.get_temperature, temperature_config=temperature_config,
 157             filename=filename, group='%stemperature/%d/' % (group, i))
 158         _sleep(calibration_config['temperature-sleep'])
 159     _LOG.debug('temperatures: %s' % Ts)
 160
 161     # get vibs
 162     vibs = _zeros((calibration_config['num-vibrations'],), dtype=_float)
 163     for i in range(calibration_config['num-vibrations']):
 164         vibs[i] = _vib(
 165             piezo=afm.piezo, vibration_config=vibration_config,
 166             filename=filename, group='%svibration/%d/' % (group, i))
 167     _LOG.debug('vibrations: %s' % vibs)
 168
 169     return (bumps, Ts, vibs)
 170
 171 def calib(afm, calibration_config, bump_config, temperature_config,
 172           vibration_config, filename=None, group='/'):
 173     """Calibrate a cantilever in one function.
 174
 175     Inputs:
 176       (see `calib_acquire()`)
 177
 178     Outputs:
 179       k    cantilever spring constant (in N/m, or equivalently nN/nm)
 180       k_s  standard deviation in our estimate of k
 181
 182     >>> import os
 183     >>> from pprint import pprint
 184     >>> import tempfile
 185     >>> from pycomedi.device import Device
 186     >>> from pycomedi.subdevice import StreamingSubdevice
 187     >>> from pycomedi.channel import AnalogChannel, DigitalChannel
 188     >>> from pycomedi.constant import AREF, IO_DIRECTION, SUBDEVICE_TYPE, UNIT
 189     >>> from pypiezo.afm import AFMPiezo
 190     >>> from pypiezo.base import PiezoAxis, InputChannel
 191     >>> from pypiezo.config import (HDF5_ChannelConfig, HDF5_AxisConfig,
 192     ...     pprint_HDF5)
 193     >>> from stepper import Stepper
 194     >>> from pyafm import AFM
 195     >>> from .config import (HDF5_CalibrationConfig, HDF5_BumpConfig,
 196     ...     HDF5_TemperatureConfig, HDF5_VibrationConfig)
 197     >>> from .analyze import calib_load_all
 198
 199     >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
 200     >>> os.close(fd)
 201
 202     >>> d = Device('/dev/comedi0')
 203     >>> d.open()
 204
 205     Setup an `AFMPiezo` instance.
 206
 207     >>> s_in = d.find_subdevice_by_type(SUBDEVICE_TYPE.ai,
 208     ...     factory=StreamingSubdevice)
 209     >>> s_out = d.find_subdevice_by_type(SUBDEVICE_TYPE.ao,
 210     ...     factory=StreamingSubdevice)
 211
 212     >>> axis_channel = s_out.channel(
 213     ...     0, factory=AnalogChannel, aref=AREF.ground)
 214     >>> input_channel = s_in.channel(0, factory=AnalogChannel, aref=AREF.diff)
 215     >>> for chan in [axis_channel, input_channel]:
 216     ...     chan.range = chan.find_range(unit=UNIT.volt, min=-10, max=10)
 217
 218     We set the minimum voltage for the `z` axis to -9 (a volt above
 219     the minimum possible voltage) to help with testing
 220     `.get_surface_position`.  Without this minimum voltage, small
 221     calibration errors could lead to a railed -10 V input for the
 222     first few surface approaching steps, which could lead to an
 223     `EdgeKink` error instead of a `FlatFit` error.
 224
 225     >>> axis_config = HDF5_AxisConfig(filename, '/bump/config/z/axis')
 226     >>> axis_config.update(
 227     ...     {'gain':20, 'sensitivity':8e-9, 'minimum':-9})
 228     >>> axis_channel_config = HDF5_ChannelConfig(
 229     ...     filename, '/bump/config/z/channel')
 230     >>> input_channel_config = HDF5_ChannelConfig(
 231     ...     filename, '/bump/config/deflection/channel')
 232
 233     >>> a = PiezoAxis(axis_config=axis_config,
 234     ...     axis_channel_config=axis_channel_config,
 235     ...     axis_channel=axis_channel, name='z')
 236     >>> a.setup_config()
 237
 238     >>> c = InputChannel(
 239     ...     channel_config=input_channel_config, channel=input_channel,
 240     ...     name='deflection')
 241     >>> c.setup_config()
 242
 243     >>> piezo = AFMPiezo(axes=[a], input_channels=[c])
 244
 245     Setup a `stepper` instance.
 246
 247     >>> s_d = d.find_subdevice_by_type(SUBDEVICE_TYPE.dio)
 248     >>> d_channels = [s_d.channel(i, factory=DigitalChannel)
 249     ...             for i in (0, 1, 2, 3)]
 250     >>> for chan in d_channels:
 251     ...     chan.dio_config(IO_DIRECTION.output)
 252
 253     >>> def write(value):
 254     ...     s_d.dio_bitfield(bits=value, write_mask=2**4-1)
 255
 256     >>> stepper = Stepper(write=write)
 257
 258     Setup an `AFM` instance.
 259
 260     >>> afm = AFM(piezo, stepper)
 261
 262     Test calibration:
 263
 264     >>> calibration_config = HDF5_CalibrationConfig(
 265     ...     filename=filename, group='/bump/config/calibration/')
 266     >>> bump_config = HDF5_BumpConfig(
 267     ...     filename=filename, group='/bump/config/bump/')
 268     >>> temperature_config = HDF5_TemperatureConfig(
 269     ...     filename=filename, group='/bump/config/temperature/')
 270     >>> vibration_config = HDF5_VibrationConfig(
 271     ...     filename=filename, group='/bump/config/vibration')
 272     >>> calib(afm, calibration_config, bump_config, temperature_config,
 273     ...     vibration_config, filename=filename, group='/')
 274     TODO: replace skipped example data with real-world values
 275     >>> pprint_HDF5(filename)  # doctest: +ELLIPSIS, +REPORT_UDIFF
 276     >>> everything = calib_load_all(filename, '/')
 277     >>> pprint(everything)
 278
 279     Close the Comedi device.
 280
 281     >>> d.close()
 282
 283     Cleanup our temporary config file.
 284
 285     os.remove(filename)
 286     """
 287     bumps, Ts, vibs = calib_acquire(
 288         afm, calibration_config, bump_config, temperature_config,
 289         vibration_config, filename=filename, group=group)
 290     # TODO: convert vib units?
 291     k,k_s = _calib_analyze(bumps, Ts, vibs)
 292     _calib_save(filename, group=group+'calibration/', bumps=bumps, Ts=Ts,
 293                 vibs=vibs, calibration_config=calibration_config, k=k, k_s=k_s)
 294     return (k, k_s)