calibcant/calibrate.py

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