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.package_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, filename=None, group='/'):
 120     """Acquire data for calibrating a cantilever in one function.
 121
 122     Inputs:
 123       afm                 a pyafm.AFM instance
 124       calibration_config  a .config._CalibrationConfig instance
 125
 126     Outputs (all are arrays of recorded data):
 127       bumps measured (V_photodiode / nm_tip) proportionality constant
 128       Ts    measured temperature (K)
 129       vibs  measured V_photodiode variance (Volts**2) in free solution
 130
 131     The temperatures are collected after moving far from the surface
 132     but before and vibrations are measured to give everything time to
 133     settle after the big move.
 134     """
 135     assert group.endswith('/'), group
 136
 137     bumps = _zeros((calibration_config['num-bumps'],), dtype=_float)
 138     for i in range(calibration_config['num-bumps']):
 139         _LOG.info('acquire bump %d of %d' % (i, calibration_config['num-bumps']))
 140         bumps[i] = _bump(afm=afm, bump_config=calibration_config['bump'],
 141                          filename=filename, group='%sbump/%d/' % (group, i))
 142     _LOG.debug('bumps: %s' % bumps)
 143
 144     move_far_from_surface(
 145         afm.stepper, distance=calibration_config['vibration-spacing'])
 146
 147     Ts = _zeros((calibration_config['num-temperatures'],), dtype=_float)
 148     for i in range(calibration_config['num-temperatures']):
 149         _LOG.info('acquire T %d of %d'
 150                  % (i, calibration_config['num-temperatures']))
 151         Ts[i] = _T(
 152             get_T=afm.get_temperature,
 153             temperature_config=calibration_config['temperature'],
 154             filename=filename, group='%stemperature/%d/' % (group, i))
 155         _sleep(calibration_config['temperature-sleep'])
 156     _LOG.debug('temperatures: %s' % Ts)
 157
 158     # get vibs
 159     vibs = _zeros((calibration_config['num-vibrations'],), dtype=_float)
 160     for i in range(calibration_config['num-vibrations']):
 161         vibs[i] = _vib(
 162             piezo=afm.piezo, vibration_config=calibration_config['vibration'],
 163             filename=filename, group='%svibration/%d/' % (group, i))
 164     _LOG.debug('vibrations: %s' % vibs)
 165
 166     return (bumps, Ts, vibs)
 167
 168 def calib(afm, calibration_config, filename=None, group='/'):
 169     """Calibrate a cantilever in one function.
 170
 171     Inputs:
 172       (see `calib_acquire()`)
 173
 174     Outputs:
 175       k    cantilever spring constant (in N/m, or equivalently nN/nm)
 176       k_s  standard deviation in our estimate of k
 177
 178     >>> import os
 179     >>> from pprint import pprint
 180     >>> import tempfile
 181     >>> from h5config.storage.hdf5 import pprint_HDF5
 182     >>> from pycomedi.device import Device
 183     >>> from pycomedi.subdevice import StreamingSubdevice
 184     >>> from pycomedi.channel import AnalogChannel, DigitalChannel
 185     >>> from pycomedi.constant import AREF, IO_DIRECTION, SUBDEVICE_TYPE, UNIT
 186     >>> from pypiezo.afm import AFMPiezo
 187     >>> from pypiezo.base import PiezoAxis, InputChannel
 188     >>> from pypiezo.config import ChannelConfig, AxisConfig
 189     >>> from stepper import Stepper
 190     >>> from pyafm import AFM
 191     >>> from .config import (CalibrationConfig, BumpConfig,
 192     ...     TemperatureConfig, VibrationConfig)
 193     >>> from .analyze import calib_load_all
 194
 195     >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
 196     >>> os.close(fd)
 197
 198     >>> d = Device('/dev/comedi0')
 199     >>> d.open()
 200
 201     Setup an `AFMPiezo` instance.
 202
 203     >>> s_in = d.find_subdevice_by_type(SUBDEVICE_TYPE.ai,
 204     ...     factory=StreamingSubdevice)
 205     >>> s_out = d.find_subdevice_by_type(SUBDEVICE_TYPE.ao,
 206     ...     factory=StreamingSubdevice)
 207
 208     >>> axis_channel = s_out.channel(
 209     ...     0, factory=AnalogChannel, aref=AREF.ground)
 210     >>> input_channel = s_in.channel(0, factory=AnalogChannel, aref=AREF.diff)
 211     >>> for chan in [axis_channel, input_channel]:
 212     ...     chan.range = chan.find_range(unit=UNIT.volt, min=-10, max=10)
 213
 214     We set the minimum voltage for the `z` axis to -9 (a volt above
 215     the minimum possible voltage) to help with testing
 216     `.get_surface_position`.  Without this minimum voltage, small
 217     calibration errors could lead to a railed -10 V input for the
 218     first few surface approaching steps, which could lead to an
 219     `EdgeKink` error instead of a `FlatFit` error.
 220
 221     >>> axis_config = AxisConfig()
 222     >>> axis_config.update(
 223     ...     {'gain':20, 'sensitivity':8e-9, 'minimum':-9})
 224     >>> axis_channel_config = ChannelConfig()
 225     >>> input_channel_config = ChannelConfig()
 226
 227     >>> a = PiezoAxis(axis_config=axis_config,
 228     ...     axis_channel_config=axis_channel_config,
 229     ...     axis_channel=axis_channel, name='z')
 230     >>> a.setup_config()
 231
 232     >>> c = InputChannel(
 233     ...     channel_config=input_channel_config, channel=input_channel,
 234     ...     name='deflection')
 235     >>> c.setup_config()
 236
 237     >>> piezo = AFMPiezo(axes=[a], input_channels=[c])
 238
 239     Setup a `stepper` instance.
 240
 241     >>> s_d = d.find_subdevice_by_type(SUBDEVICE_TYPE.dio)
 242     >>> d_channels = [s_d.channel(i, factory=DigitalChannel)
 243     ...             for i in (0, 1, 2, 3)]
 244     >>> for chan in d_channels:
 245     ...     chan.dio_config(IO_DIRECTION.output)
 246
 247     >>> def write(value):
 248     ...     s_d.dio_bitfield(bits=value, write_mask=2**4-1)
 249
 250     >>> stepper = Stepper(write=write)
 251
 252     Setup an `AFM` instance.
 253
 254     >>> afm = AFM(piezo, stepper)
 255
 256     Test calibration:
 257
 258     >>> calibration_config = CalibrationConfig()
 259     >>> bump_config = BumpConfig()
 260     >>> temperature_config = TemperatureConfig()
 261     >>> vibration_config = VibrationConfig()
 262     >>> calib(afm, calibration_config, bump_config, temperature_config,
 263     ...     vibration_config, filename=filename, group='/')
 264     TODO: replace skipped example data with real-world values
 265     >>> pprint_HDF5(filename)  # doctest: +ELLIPSIS, +REPORT_UDIFF
 266     >>> everything = calib_load_all(filename, '/')
 267     >>> pprint(everything)
 268
 269     Close the Comedi device.
 270
 271     >>> d.close()
 272
 273     Cleanup our temporary config file.
 274
 275     os.remove(filename)
 276     """
 277     bumps, Ts, vibs = calib_acquire(
 278         afm, calibration_config, bump_config, temperature_config,
 279         vibration_config, filename=filename, group=group)
 280     # TODO: convert vib units?
 281     k,k_s = _calib_analyze(bumps, Ts, vibs)
 282     _calib_save(filename, group=group+'calibration/', bumps=bumps,
 283                 temperatures=Ts, vibrations=vibs,
 284                 calibration_config=calibration_config, k=k, k_s=k_s)
 285     return (k, k_s)