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.afm 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     >>> axis_channel_config['name'] = 'z'
 226     >>> axis_config['channel'] = axis_channel_config
 227     >>> input_channel_config = ChannelConfig()
 228     >>> input_channel_config['name'] = 'deflection'
 229
 230     >>> a = PiezoAxis(config=axis_config, axis_channel=axis_channel)
 231     >>> a.setup_config()
 232
 233     >>> c = InputChannel(config=input_channel_config, channel=input_channel)
 234     >>> c.setup_config()
 235
 236     >>> piezo = AFMPiezo(axes=[a], inputs=[c])
 237
 238     Setup a `stepper` instance.
 239
 240     >>> s_d = d.find_subdevice_by_type(SUBDEVICE_TYPE.dio)
 241     >>> d_channels = [s_d.channel(i, factory=DigitalChannel)
 242     ...             for i in (0, 1, 2, 3)]
 243     >>> for chan in d_channels:
 244     ...     chan.dio_config(IO_DIRECTION.output)
 245
 246     >>> def write(value):
 247     ...     s_d.dio_bitfield(bits=value, write_mask=2**4-1)
 248
 249     >>> stepper = Stepper(write=write)
 250
 251     Setup an `AFM` instance.
 252
 253     >>> afm = AFM(piezo, stepper)
 254
 255     Test calibration:
 256
 257     >>> calibration_config = CalibrationConfig()
 258     >>> calibration_config['bump'] = BumpConfig()
 259     >>> calibration_config['temperature'] = TemperatureConfig()
 260     >>> calibration_config['vibration'] = VibrationConfig()
 261     >>> calib(afm, calibration_config, filename=filename, group='/')
 262     TODO: replace skipped example data with real-world values
 263     >>> pprint_HDF5(filename)  # doctest: +ELLIPSIS, +REPORT_UDIFF
 264     /
 265       /bump
 266         /bump/0
 267           /bump/0/config
 268             /bump/0/config/bump
 269               <HDF5 dataset "far-steps": shape (), type "<i4">
 270                 200
 271               ...
 272             /bump/0/config/deflection
 273               /bump/0/config/deflection/channel
 274                 <HDF5 dataset "channel": shape (), type "<i4">
 275                   0
 276                 ...
 277             /bump/0/config/z
 278               /bump/0/config/z/axis
 279                 /bump/0/config/z/axis/channel
 280                   <HDF5 dataset "channel": shape (), type "<i4">
 281                     0
 282                   ...
 283                 <HDF5 dataset "gain": shape (), type "<i4">
 284                   20
 285                 ...
 286           <HDF5 dataset "processed": shape (), type "<f8">
 287             ...
 288           /bump/0/raw
 289             <HDF5 dataset "deflection": shape (2048,), type "<u2">
 290               [...]
 291             <HDF5 dataset "z": shape (2048,), type "<u2">
 292               [...]
 293         /bump/1
 294         ...
 295       /calibration
 296         /calibration/config
 297           /calibration/config/bump
 298             <HDF5 dataset "far-steps": shape (), type "<i4">
 299               200
 300             ...
 301           <HDF5 dataset "num-bumps": shape (), type "<i4">
 302             10
 303           ...
 304         /calibration/processed
 305           /calibration/processed/spring-constant
 306             <HDF5 dataset "data": shape (), type "<f8">
 307               ...
 308             <HDF5 dataset "standard-deviation": shape (), type "<f8">
 309               ...
 310             <HDF5 dataset "units": shape (), type "|S3">
 311               N/m
 312         /calibration/raw
 313           /calibration/raw/photodiode-sensitivity
 314             <HDF5 dataset "data": shape (10,), type "<f8">
 315               [...]
 316             <HDF5 dataset "units": shape (), type "|S3">
 317               V/m
 318           /calibration/raw/temperature
 319             <HDF5 dataset "data": shape (10,), type "<f8">
 320               [...]
 321             <HDF5 dataset "units": shape (), type "|S1">
 322               K
 323           /calibration/raw/thermal-vibration-variance
 324             <HDF5 dataset "data": shape (20,), type "<f8">
 325               [...]
 326             <HDF5 dataset "units": shape (), type "|S3">
 327               V^2
 328       /temperature
 329         /temperature/0
 330           /temperature/0/config
 331             <HDF5 dataset "default": shape (), type "|b1">
 332               False
 333             <HDF5 dataset "units": shape (), type "|S7">
 334               Celsius
 335           <HDF5 dataset "processed": shape (), type "<f8">
 336             295.15
 337           <HDF5 dataset "raw": shape (), type "<i4">
 338             22
 339         /temperature/1
 340         ...
 341       /vibration
 342         /vibration/0
 343           /vibration/0/config
 344             /vibration/0/config/deflection
 345               <HDF5 dataset "channel": shape (), type "<i4">
 346                 0
 347               ...
 348             /vibration/0/config/vibration
 349               <HDF5 dataset "chunk-size": shape (), type "<i4">
 350                 2048
 351               ...
 352           <HDF5 dataset "processed": shape (), type "<f8">
 353             ...
 354           /vibration/0/raw
 355             <HDF5 dataset "deflection": shape (65536,), type "<u2">
 356               [...]
 357         /vibration/1
 358         ...
 359         /vibration/19
 360           ...
 361           /vibration/19/raw
 362             <HDF5 dataset "deflection": shape (65536,), type "<u2">
 363               [...]
 364     >>> everything = calib_load_all(filename, '/')
 365     >>> pprint(everything)  # doctest: +ELLIPSIS, +REPORT_UDIFF
 366     {'bump_details': [{'bump_config': <BumpConfig ...>,
 367                        'deflection_channel_config': <ChannelConfig ...>,
 368                        'processed_bump': ...,
 369                        'raw_bump': {'deflection': array([...], dtype=uint16),
 370                                     'z': array([...], dtype=uint16)},
 371                        'z_axis_config': <AxisConfig ...>},
 372                       ...],
 373      'bumps': array([...]),
 374      'calibration_config': <CalibrationConfig ...>,
 375      'k': ...,
 376      'k_s': ...,
 377      'temperature_details': [{'processed_temperature': ...,
 378                               'raw_temperature': array(22),
 379                               'temperature_config': <TemperatureConfig ...>},
 380                              ...],
 381      'temperatures': array([...]),
 382      'vibration_details': [{'deflection_channel_config': <ChannelConfig ...>,
 383                             'processed_vibration': ...,
 384                             'raw_vibration': array([...], dtype=uint16),
 385                             'vibration_config': <VibrationConfig ...>},
 386                            ...],
 387      'vibrations': array([...])}
 388
 389     Close the Comedi device.
 390
 391     >>> d.close()
 392
 393     Cleanup our temporary config file.
 394
 395     os.remove(filename)
 396     """
 397     bumps, Ts, vibs = calib_acquire(
 398         afm, calibration_config, filename=filename, group=group)
 399     # TODO: convert vib units?
 400     k,k_s = _calib_analyze(bumps, Ts, vibs)
 401     _calib_save(filename, group=group+'calibration/', bumps=bumps,
 402                 temperatures=Ts, vibrations=vibs,
 403                 calibration_config=calibration_config, k=k, k_s=k_s)
 404     return (k, k_s)