-*.pyc
+MANIFEST
build/
dist/
calibcant.egg-info/
+*.pyc
Calculates the spring constant `k` of a cantilever using the
equipartition theorem:
- 1/2 k_B T = 1/2 K x^2
+
+.. math:: \frac{1}{2} k_B T = \frac{1}{2} K x^2
+
The analysis uses the expected power spectral density (PSD) of a
damped simple harmonic oscillator to filter out noise from the
measured PSD.
+Installation
+============
+
Packages
-========
+--------
Gentoo
-------
+~~~~~~
I've packaged calibcant for Gentoo. You need layman_ and my `wtk
overlay`_. Install with::
Dependencies
-============
+------------
-calibcant requires the following Python modules:
+If you're installing by hand or packaging calibcant for another
+distribution, you'll need the following dependencies:
-* linfit (depends in turn on scipy.stats.linregress)
-* TODO
+=========== ================= =====================
+Package Debian_ Gentoo_
+=========== ================= =====================
+Numpy_ python-numpy dev-python/numpy
+Scipy_ python-scipy sci-libs/scipy
+H5Py_ python-h5py dev-python/h5py
+Matplotlib_ python-matplotlib dev-python/matplotlib
+Nose_ python-nose dev-python/nose
+FFT_tools_ dev-python/FFT-tools
+Pypiezo_ sci-libs/pypiezo
+=========== ================= =====================
-* Numpy_
-* Scipy_
-* PyYAML_ (for saving and loading playlists)
-* Matplotlib_ (for generating plots)
-* wxPython_ (for the GUI)
+You'll also need my pyafm_ and stepper_ packages, or suitable
+replacements.
-.. _Python: http://www.python.org/
-.. _Numpy: http://numpy.scipy.org/
-.. _Scipy: http://www.scipy.org/
-.. _PyYAML: http://pyyaml.org/
-.. _Matplotlib: http://matplotlib.sourceforge.net/
-.. _wxPython: http://www.wxpython.org/
+Installing by hand
+------------------
+Calibcant is available as a Git_ repository::
-Getting the source
-==================
+ $ git clone http://www.physics.drexel.edu/~wking/code/git/calibcant.git
-calibcant is available as a Git_ repository::
+See the homepage_ for details. To install the checkout, run the
+standard::
- $ git clone http://www.physics.drexel.edu/~wking/code/git/calibcant.git
+ $ python setup.py install
-There are also periodic bundled releases. For example, get version
-0.5 as a gzipped tarball with::
- $ wget http://www.physics.drexel.edu/~wking/code/python/calibcant-0.5.tar.gz
- $ tar -xzvf calibcant-0.5.tar.gz
+Usage
+=====
-.. _Git: http://git-scm.com/
+See the module docstrings for simple examples.
-Installation
-============
+Testing
+=======
-After downloading, change to the source directory and run::
+Run internal unit tests with::
- $ python setup.py install
+ $ nosetests --with-doctest --doctest-tests calibcant
-to install calibcant. Run::
- $ python setup.py install --help
+Licence
+=======
-to see a list of installation options you may want to configure.
+This project is distributed under the `GNU General Public License
+Version 3`_ or greater.
-Usage
-=====
+Author
+======
-You can... While there are a number of command-line scripts
-TODO
+W. Trevor King
+wking@drexel.edu
+Copyright 2007-2011
+
+
+.. _layman: http://layman.sourceforge.net/
+.. _wtk overlay:
+ http://www.physics.drexel.edu/~wking/unfolding-disasters/posts/Gentoo_overlay
+.. _Debian: http://www.debian.org/
+.. _Gentoo: http://www.gentoo.org/
+.. _NumPy: http://numpy.scipy.org/
+.. _SciPy: http://www.scipy.org/
+.. _H5Py: http://code.google.com/p/h5py/
+.. _Matplotlib: http://matplotlib.sourceforge.net/
+.. _Nose: http://somethingaboutorange.com/mrl/projects/nose/
+.. _FFT_tools:
+ http://www.physics.drexel.edu/~wking/unfolding-disasters/posts/FFT-tools/
+.. _Pypiezo:
+ http://www.physics.drexel.edu/~wking/unfolding-disasters/posts/pypiezo/
+.. _Git: http://git-scm.com/
+.. _homepage:
+ http://www.physics.drexel.edu/~wking/unfolding-disasters/posts/calibcant/
+.. _GNU General Public License Version 3: http://www.gnu.org/licenses/gpl.txt
--- /dev/null
+#!/usr/bin/env python
+# calibcant - tools for thermally calibrating AFM cantilevers
+#
+# Copyright (C) 2008-2011 W. Trevor King <wking@drexel.edu>
+#
+# This file is part of calibcant.
+#
+# calibcant is free software: you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation, either
+# version 3 of the License, or (at your option) any later version.
+#
+# calibcant is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with calibcant. If not, see
+# <http://www.gnu.org/licenses/>.
+
+"""Load a calibration file and (re)analyze from the raw data.
+"""
+
+from optparse import OptionParser
+
+from matplotlib.pyplot import show
+from calibcant.analyze import calib_analyze_all
+
+
+def main(args):
+ usage = '%prog [options] filename'
+ p = OptionParser(usage)
+
+ p.add_option('-g', '--group', dest='group', default='/',
+ help='HDF5 group containing calibration data (%default).')
+ p.add_option('-d', '--dry-run', dest='dry_run', action='store_true',
+ help="Don't change the original file.")
+
+ options,args = p.parse_args(args)
+ filename = args[0]
+
+ k,k_s = calib_analyze_all(
+ filename=filename, group=options.group, dry_run=options.dry_run)
+ print "New spring constant:"
+ print "%g +/- %g N/m" % (k, k_s)
+
+
+if __name__ == '__main__':
+ import sys
+
+ sys.exit(main(sys.argv[1:]))
--- /dev/null
+#!/usr/bin/env python
+# calibcant - tools for thermally calibrating AFM cantilevers
+#
+# Copyright (C) 2008-2011 W. Trevor King <wking@drexel.edu>
+#
+# This file is part of calibcant.
+#
+# calibcant is free software: you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation, either
+# version 3 of the License, or (at your option) any later version.
+#
+# calibcant is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with calibcant. If not, see
+# <http://www.gnu.org/licenses/>.
+
+"""Load a calibration file and plot everything interesting.
+"""
+
+from optparse import OptionParser
+
+from matplotlib.pyplot import close, get_fignums, figure, show
+from calibcant.analyze import calib_load_all, calib_plot_all
+
+
+def main(args):
+ usage = '%prog [options] filename'
+ p = OptionParser(usage)
+
+ p.add_option('-g', '--group', dest='group', default='/',
+ help='HDF5 group containing calibration data (%default).')
+ p.add_option('-b', '--no-bumps', dest='bumps', default=True,
+ action='store_false',
+ help="Don't display bump details.")
+ p.add_option('-v', '--no-vibrations', dest='vibrations', default=True,
+ action='store_false',
+ help="Don't display vibration details.")
+ p.add_option('-s', '--save-figures', dest='save', action='store_true',
+ help="Save plots (instead of displaying them.")
+
+ options,args = p.parse_args(args)
+ filename = args[0]
+
+ stuff = calib_load_all(filename=filename, group=options.group)
+ if not options.bumps:
+ stuff['bump_details'] = []
+ if not options.vibrations:
+ stuff['vibration_details'] = []
+ calib_plot_all(**stuff)
+ if options.save:
+ for i in get_fignums():
+ fig = figure(i)
+ fig.savefig('%i.png' % i, dpi=300)
+ close(i)
+ else:
+ show()
+
+
+if __name__ == '__main__':
+ import sys
+
+ sys.exit(main(sys.argv[1:]))
--- /dev/null
+# T family.
+# Fairly stubby, since a one shot Temp measurement is a common thing.
+# We just wrap that to provide a consistent interface.
+
+from . import LOG as _LOG
+from . import base_config as _base_config
+from .T_analyze import T_analyze as _T_analyze
+from .T_analyze import T_save as _T_save
+
+
+def T_acquire(get_T=None):
+ """Measure the current temperature of the sample,
+
+ If `get_T` is `None`, fake it by returning
+ `base_config['temperature']`.
+ """
+ if get_T:
+ _LOG.info('measure temperature')
+ T = get_T()
+ else:
+ T = None
+ if not T:
+ _LOG.info('fake temperature %g' % _base_config['temperature'])
+ return _base_config['temperature']
+
+def T(get_T, temperature_config, filename, group='/'):
+ """Wrapper around T_acquire(), T_analyze(), T_save().
+
+ >>> import os
+ >>> import tempfile
+ >>> from pypiezo.config import pprint_HDF5
+ >>> from .config import HDF5_TemperatureConfig
+
+ >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
+ >>> os.close(fd)
+
+ >>> temperature_config = HDF5_TemperatureConfig(
+ ... filename=filename, group='/T/config/')
+ >>> def get_T():
+ ... return 19.2
+
+ >>> t = T(get_T=get_T, temperature_config=temperature_config,
+ ... filename=filename, group='/T/')
+ >>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
+ /
+ /T
+ /T/config
+ <HDF5 dataset "default": shape (), type "|S2">
+ no
+ <HDF5 dataset "units": shape (), type "|S7">
+ Celsius
+ <HDF5 dataset "processed": shape (), type "<f8">
+ 292.35
+ <HDF5 dataset "raw": shape (), type "<f8">
+ 19.2
+ >>> t = T(get_T=None, temperature_config=temperature_config,
+ ... filename=filename, group='/T/')
+ >>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
+ /
+ /T
+ /T/config
+ <HDF5 dataset "default": shape (), type "|S3">
+ yes
+ <HDF5 dataset "units": shape (), type "|S7">
+ Celsius
+ <HDF5 dataset "processed": shape (), type "<f8">
+ 295.15
+ <HDF5 dataset "raw": shape (), type "<i4">
+ 22
+
+ Cleanup our temporary config file.
+
+ >>> os.remove(filename)
+ """
+ T_raw = T_acquire(get_T)
+ T_ret = _T_analyze(T_raw, temperature_config)
+ temperature_config['default'] = not get_T
+ _T_save(filename, group=group, raw_T=T_raw,
+ temperature_config=temperature_config, processed_T=T_ret)
+ return T_ret
# License along with calibcant. If not, see
# <http://www.gnu.org/licenses/>.
+"""Temperature analysis.
+
+Separate the more general `T_analyze()` from the other `T_*()`
+functions in calibcant.
+
+The relevant physical quantities are:
+
+* `T` Temperature at which thermal vibration measurements were acquired
+
+>>> import os
+>>> import tempfile
+>>> import numpy
+>>> from .config import HDF5_TemperatureConfig
+>>> from pypiezo.config import pprint_HDF5
+
+>>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
+>>> os.close(fd)
+
+>>> temperature_config = HDF5_TemperatureConfig(
+... filename=filename, group='/T/config/')
+
+>>> raw_T = numpy.array([22, 23.5, 24])
+>>> processed_T = T_analyze(raw_T, temperature_config)
+>>> T_plot(raw_T=raw_T, processed_T=processed_T)
+>>> T_save(filename=filename, group='/T/', raw_T=raw_T,
+... temperature_config=temperature_config, processed_T=processed_T)
+
+>>> pprint_HDF5(filename) # doctest: +REPORT_UDIFF
+/
+ /T
+ /T/config
+ <HDF5 dataset "default": shape (), type "|S3">
+ yes
+ <HDF5 dataset "units": shape (), type "|S7">
+ Celsius
+ <HDF5 dataset "processed": shape (3,), type "<f8">
+ [ 295.15 296.65 297.15]
+ <HDF5 dataset "raw": shape (3,), type "<f8">
+ [ 22. 23.5 24. ]
+
+>>> raw_T,temperature_config,processed_T = T_load(
+... filename=filename, group='/T/')
+>>> print temperature_config.dump()
+units: Celsius
+default: yes
+>>> raw_T
+array([ 22. , 23.5, 24. ])
+>>> type(raw_T)
+<type 'numpy.ndarray'>
+>>> processed_T
+array([ 295.15, 296.65, 297.15])
+
+>>> os.remove(filename)
"""
-Separate the more general T_analyze() from the other T_*()
-functions in calibcant. Also provide a command line interface
-for analyzing data acquired through other workflows.
-The relevant physical quantities are :
- T Temperature at which thermal vibration measurements were aquired
-"""
-
-import numpy
+import h5py as _h5py
+from scipy.constants import C2K as _C2K
-import data_logger
-from splittable_kwargs import splittableKwargsFunction, \
- make_splittable_kwargs_function
+try:
+ import matplotlib as _matplotlib
+ import matplotlib.pyplot as _matplotlib_pyplot
+ import time as _time # for timestamping lines on plots
+except (ImportError, RuntimeError), e:
+ _matplotlib = None
+ _matplotlib_import_error = e
-from . import common
-from . import config
+from pypiezo.config import h5_create_group as _h5_create_group
+from . import LOG as _LOG
+from . import base_config as _base_config
+from .config import Celsius as _Celsius
+from .config import Kelvin as _Kelvin
+from .config import HDF5_TemperatureConfig as _HDF5_TemperatureConfig
-def C_to_K(celsius) :
- "Convert Celsius -> Kelvin."
- return celsius + 273.15
-def K_to_K(kelvin) :
- "Convert Kelvin -> Kelvin."
- return kelvin
+def T_analyze(T, temperature_config):
+ """Convert measured temperature to Kelvin.
-@splittableKwargsFunction()
-def T_analyze(T, convert_to_K=C_to_K) :
- """
- Not much to do here, just convert to Kelvin.
- Uses C_to_K (defined above) by default.
- """
- try : # if T is an array, convert element by element
- for i in range(len(T)) :
- T[i] = convert_to_K(T[i])
- except TypeError : # otherwise, make an array from a single T
- T = [convert_to_K(T)]
- return T
-
-@splittableKwargsFunction()
-def T_save(T, log_dir=None) :
- """
- Save either a single T (if you are crazy :p),
- or an array of them (more reasonable)
- """
- T = numpy.array(T, dtype=numpy.float)
- if log_dir != None :
- log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
- log_name="T_float")
- log.write_binary(T.tostring())
- if config.LOG_DATA != None :
- log = data_logger.data_log(config.LOG_DIR, noclobber_logsubdir=False,
- log_name="T_float")
- log.write_binary(T.tostring())
-
-def T_load(datafile=None) :
- """
- Load the saved T array (possibly of length 1), and return it. If
- datafile == None, return an array of one config.DEFAULT_TEMP
- instead.
- """
- if datafile == None :
- return numpy.array([config.DEFAULT_TEMP], dtype=numpy.float)
- else :
- dl = data_logger.data_load()
- return dl.read_binary(datafile)
-
-@splittableKwargsFunction()
-def T_plot(T, plotVerbose=False) :
- """
- Umm, just print the temperature?
+ `T` should be a numpy ndarray or scalar. `temperature_config`
+ should be a `config._TemperatureConfig` instance.
"""
- if plotVerbose or config.PYLAB_VERBOSE or config.TEXT_VERBOSE :
- print "Temperature ", T
-
-@splittableKwargsFunction((T_analyze, 'T', 'convert_to_K'),
- (T_plot, 'T'))
-def T_load_analyze_tweaked(tweak_file, convert_to_K=C_to_K, textVerboseFile=None, **kwargs) :
- "Load all the T array files from a tweak file and return a single array"
- T_analyze_kwargs,T_plot_kwargs = \
- T_load_analyze_tweaked._splitargs(T_load_analyze_tweaked, kwargs)
- Ts = []
- for line in file(tweak_file, 'r') :
- parsed = line.split()
- path = parsed[0].strip()
- if path[0] == '#': # a comment
- continue
- if textVerboseFile != None :
- print >> textVerboseFile, "Reading data from %s" % (path)
- # read the data
- data = T_load(path)
- Ts.extend(data)
- T_analyze(Ts, convert_to_K=convert_to_K)
- return numpy.array(Ts, dtype=numpy.float)
-
-# commandline interface functions
-import scipy.io, sys
-
-def read_data(ifile):
- "ifile can be a filename string or open (seekable) file object"
- if ifile == None : ifile = sys.stdin
- unlabeled_data=scipy.io.read_array(ifile)
- return unlabeled_data
-
-if __name__ == '__main__' :
- # command line interface
- from optparse import OptionParser
-
- usage_string = ('%prog <input-file>\n'
- '2008, W. Trevor King.\n'
- '\n'
- 'There are two operation modes, one to analyze a single T (temperature) file,\n'
- 'and one to analyze tweak files.\n'
- '\n'
- 'Single file mode (the default) :\n'
- 'Reads in single column ASCII file of temperatures and... prints them back out.\n'
- 'No need to do this, but the option is available for consistency with the other\n'
- 'calibcant modules.\n'
- '\n'
- 'Tweak file mode:\n'
- 'Runs the same analysis as in single file mode for each T file in\n'
- 'a tweak file. Each line in the tweak file specifies a single T file.\n'
- 'Blank lines and those beginning with a pound sign (#) are ignored.\n'
- 'A T file contains a sequence of 32 bit floats representing temperature in K.\n'
- )
- parser = OptionParser(usage=usage_string, version='%prog 0.1')
- parser.add_option('-C', '--celsius', dest='celsius',
- help='Use Celsius input temperatures instead of Kelvin (default %default)\n',
- action='store_true', default=False)
- parser.add_option('-o', '--output-file', dest='ofilename',
- help='write output to FILE (default stdout)',
- type='string', metavar='FILE')
- parser.add_option('-c', '--comma-out', dest='comma_out', action='store_true',
- help='Output comma-seperated values (default %default)',
- default=False)
- parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
- help='Run in tweak-file mode',
- default=False)
- parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
- help='Print lots of debugging information',
- default=False)
-
- options,args = parser.parse_args()
- parser.destroy()
- assert len(args) >= 1, "Need an input file"
-
- ifilename = args[0]
-
- if options.ofilename != None :
- ofile = file(options.ofilename, 'w')
- else :
- ofile = sys.stdout
- if options.verbose == True :
- vfile = sys.stderr
- else :
- vfile = None
- config.TEXT_VERBOSE = options.verbose
- config.PYLAB_VERBOSE = False
- config.GNUPLOT_VERBOSE = False
- if options.celsius :
- convert_to_K = C_to_K
- else :
- convert_to_K = K_to_K
-
- if options.tweakmode == False :
- data = read_data(ifilename)
- Ts = T_analyze(data, convert_to_K)
- else : # tweak file mode
- Ts = T_load_analyze_tweaked(ifilename, convert_to_K, textVerboseFile=vfile)
-
- if options.comma_out :
- sep = ','
- else :
- sep = '\n'
- common.write_array(ofile, Ts, sep)
-
- if options.ofilename != None :
- ofile.close()
+ if temperature_config['units'] == _Celsius:
+ return _C2K(T)
+ else:
+ return T
+
+def T_save(filename, group='/', raw_T=None, temperature_config=None,
+ processed_T=None):
+ with _h5py.File(filename, 'a') as f:
+ cwg = _h5_create_group(f, group)
+ if raw_T is not None:
+ try:
+ del cwg['raw']
+ except KeyError:
+ pass
+ cwg['raw'] = raw_T
+ if temperature_config:
+ config_cwg = _h5_create_group(cwg, 'config')
+ temperature_config.save(group=config_cwg)
+ if processed_T is not None:
+ try:
+ del cwg['processed']
+ except KeyError:
+ pass
+ cwg['processed'] = processed_T
+
+def T_load(filename, group='/'):
+ assert group.endswith('/')
+ raw_T = processed_T = None
+ with _h5py.File(filename, 'a') as f:
+ try:
+ raw_T = f[group+'raw'][...]
+ except KeyError:
+ pass
+ temperature_config = _HDF5_TemperatureConfig(
+ filename=filename, group=group+'config/')
+ try:
+ processed_T = f[group+'processed'][...]
+ except KeyError:
+ pass
+ temperature_config.load()
+ return (raw_T, temperature_config, processed_T)
+
+def T_plot(raw_T=None, processed_T=None):
+ if not _matplotlib:
+ raise _matplotlib_import_error
+ figure = _matplotlib_pyplot.figure()
+ timestamp = _time.strftime('%H%M%S')
+ if raw_T is None:
+ if processed_T is None:
+ return # nothing to plot
+ axes1 = None
+ axes2 = figure.add_subplot(1, 1, 1)
+ elif processed_T is None:
+ axes1 = figure.add_subplot(1, 1, 1)
+ axes2 = None
+ else:
+ axes1 = figure.add_subplot(2, 1, 1)
+ axes2 = figure.add_subplot(2, 1, 2)
+ if axes1:
+ axes1.set_title('Raw Temperatures %s' % timestamp)
+ axes1.plot(raw_T, label='raw')
+ if axes2:
+ axes2.set_title('Processed Temperatures %s' % timestamp)
+ axes2.plot(processed_T, label='processed')
+ figure.show()
# License along with calibcant. If not, see
# <http://www.gnu.org/licenses/>.
-__version__ = '0.5'
+import logging as _logging
+import logging.handlers as _logging_handlers
+
+
+__version__ = '0.4'
+
+
+LOG = _logging.getLogger('calibcant')
+"Calibcant logger"
+
+LOG.setLevel(_logging.WARN)
+_formatter = _logging.Formatter(
+ '%(asctime)s - %(name)s - %(levelname)s - %(message)s')
+
+_stream_handler = _logging.StreamHandler()
+_stream_handler.setLevel(_logging.DEBUG)
+_stream_handler.setFormatter(_formatter)
+LOG.addHandler(_stream_handler)
+
+_syslog_handler = None
+
+
+from .config import _BaseConfig
+from .config import find_base_config as _find_base_config
+
+
+def setup_base_config(config):
+ global base_config, _syslog_handler
+ base_config = config
+
+ LOG.setLevel(base_config['log-level'])
+
+ if base_config['syslog']:
+ if not _syslog_handler:
+ _syslog_handler = _logging_handlers.SysLogHandler()
+ _syslog_handler.setLevel(_logging.DEBUG)
+ LOG.handlers = [_syslog_handler]
+ else:
+ LOG.handlers = [_stream_handler]
+
+ LOG.info('setup base_config:\n%s' % config.dump())
+
+def clear_base_config():
+ setup_base_config(_BaseConfig())
+
+base_config = _find_base_config()
+setup_base_config(base_config)
# License along with calibcant. If not, see
# <http://www.gnu.org/licenses/>.
-"""
-Separate the more general calib_analyze() from the other calib_*()
-functions in calibcant. Also provide a command line interface
-for analyzing data acquired through other workflows.
-
-The relevent physical quantities are :
- Vzp_out Output z-piezo voltage (what we generate)
- Vzp Applied z-piezo voltage (after external ZPGAIN)
- Zp The z-piezo position
- Zcant The cantilever vertical deflection
- Vphoto The photodiode vertical deflection voltage (what we measure)
- Fcant The force on the cantilever
- T The temperature of the cantilever and surrounding solution
- (another thing we measure)
- k_b Boltzmann's constant
-
-Which are related by the parameters :
- zpGain Vzp_out / Vzp
- zpSensitivity Zp / Vzp
- photoSensitivity Vphoto / Zcant
- k_cant Fcant / Zcant
+"""Calculate `k` from arrays of bumps, temperatures, and vibrations.
+
+Separate the more general `calib_analyze()` from the other `calib_*()`
+functions in calibcant.
+
+The relevent physical quantities are :
+ Vzp_out Output z-piezo voltage (what we generate)
+ Vzp Applied z-piezo voltage (after external ZPGAIN)
+ Zp The z-piezo position
+ Zcant The cantilever vertical deflection
+ Vphoto The photodiode vertical deflection voltage (what we measure)
+ Fcant The force on the cantilever
+ T The temperature of the cantilever and surrounding solution
+ (another thing we measure)
+ k_b Boltzmann's constant
+
+Which are related by the parameters:
+ zp_gain Vzp_out / Vzp
+ zp_sensitivity Zp / Vzp
+ photo_sensitivity Vphoto / Zcant
+ k_cant Fcant / Zcant
+
+
+>>> import os
+>>> import tempfile
+>>> import numpy
+>>> from pypiezo.config import pprint_HDF5
+>>> from .config import HDF5_CalibrationConfig
+
+>>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
+>>> os.close(fd)
+
+>>> calibration_config = HDF5_CalibrationConfig(filename, '/calib/config/')
+>>> bumps = numpy.array((15.9e6, 16.9e6, 16.3e6))
+>>> temperatures = numpy.array((295, 295.2, 294.8))
+>>> vibrations = numpy.array((2.20e-5, 2.22e-5, 2.21e-5))
+
+>>> k,k_s = calib_analyze(bumps=bumps, temperatures=temperatures,
+... vibrations=vibrations)
+>>> (k, k_s) # doctest: +ELLIPSIS
+(0.0493..., 0.00248...)
+
+Most of the error in this example comes from uncertainty in the
+photodiode sensitivity (bumps).
+
+>>> k_s/k # doctest: +ELLIPSIS
+0.0503...
+>>> bumps.std()/bumps.mean() # doctest: +ELLIPSIS
+0.0251...
+>>> temperatures.std()/temperatures.mean() # doctest: +ELLIPSIS
+0.000553...
+>>> vibrations.std()/vibrations.mean() # doctest: +ELLIPSIS
+0.00369...
+
+>>> calib_save(filename=filename, group='/calib/',
+... bumps=bumps, temperatures=temperatures, vibrations=vibrations,
+... calibration_config=calibration_config, k=k, k_s=k_s)
+>>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
+/
+ /calib
+ /calib/config
+ <HDF5 dataset "num-bumps": shape (), type "|S2">
+ 10
+ <HDF5 dataset "num-temperatures": shape (), type "|S2">
+ 10
+ <HDF5 dataset "num-vibrations": shape (), type "|S2">
+ 20
+ <HDF5 dataset "temperature-sleep": shape (), type "|S1">
+ 1
+ <HDF5 dataset "vibration-spacing": shape (), type "|S5">
+ 5e-05
+ /calib/processed
+ /calib/processed/spring-constant
+ <HDF5 dataset "data": shape (), type "<f8">
+ 0.0493...
+ <HDF5 dataset "standard-deviation": shape (), type "<f8">
+ 0.00248...
+ <HDF5 dataset "units": shape (), type "|S3">
+ N/m
+ /calib/raw
+ /calib/raw/photodiode-sensitivity
+ <HDF5 dataset "data": shape (3,), type "<f8">
+ [ 15900000. 16900000. 16300000.]
+ <HDF5 dataset "units": shape (), type "|S3">
+ V/m
+ /calib/raw/temperature
+ <HDF5 dataset "data": shape (3,), type "<f8">
+ [ 295. 295.2 294.8]
+ <HDF5 dataset "units": shape (), type "|S1">
+ K
+ /calib/raw/thermal-vibration-variance
+ <HDF5 dataset "data": shape (3,), type "<f8">
+ [ 2.20...e-05 2.220...e-05 2.210...e-05]
+ <HDF5 dataset "units": shape (), type "|S3">
+ V^2
+
+>>> bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load(
+... filename=filename, group='/calib/')
+>>> (k, k_s) # doctest: +ELLIPSIS
+(0.0493..., 0.00248...)
+
+>>> os.remove(filename)
"""
-import numpy
+import h5py as _h5py
+import numpy as _numpy
+try:
+ from scipy.constants import Boltzmann as _kB # in J/K
+except ImportError:
+ from scipy.constants import Bolzmann as _kB # in J/K
+# Bolzmann -> Boltzmann patch submitted:
+# http://projects.scipy.org/scipy/ticket/1417
+# Fixed in scipy commit 4716d91, Apr 2, 2011, during work after v0.9.0rc5.
-import data_logger
-from splittable_kwargs import splittableKwargsFunction, \
- make_splittable_kwargs_function
+try:
+ import matplotlib as _matplotlib
+ import matplotlib.pyplot as _matplotlib_pyplot
+ import time as _time # for timestamping lines on plots
+except (ImportError, RuntimeError), e:
+ _matplotlib = None
+ _matplotlib_import_error = e
-from . import common
-from . import config
-from . import T_analyze
+from pypiezo.config import h5_create_group as _h5_create_group
+from . import LOG as _LOG
+from . import base_config as _base_config
+from .bump_analyze import bump_analyze as _bump_analyze
+from .bump_analyze import bump_load as _bump_load
+from .bump_analyze import bump_save as _bump_save
+from .config import HDF5_CalibrationConfig as _HDF5_CalibrationConfig
+from .T_analyze import T_analyze as _temperature_analyze
+from .T_analyze import T_load as _temperature_load
+from .T_analyze import T_save as _temperature_save
+from .vib_analyze import vib_analyze as _vibration_analyze
+from .vib_analyze import vib_load as _vibration_load
+from .vib_analyze import vib_save as _vibration_save
-kb = 1.3806504e-23 # Boltzmann's constant in J/K
-#@splittableKwargsFunction((calib_plot, 'bumps', 'Ts', 'vibs'))
-# Some of the child functions aren't yet defined, so postpone
-# make-splittable until later in the module.
-def calib_analyze(bumps, Ts, vibs, **kwargs) :
- """
- Analyze data from get_calibration_data()
- return (k, k_s, !!!a2_r, T_r, one_o_Vp2_r)
- Inputs (all are arrays of recorded data) :
- bumps measured (V_photodiode / nm_tip) proportionality constant
- Ts measured temperature (K)
- vibs measured V_photodiode variance in free solution (V**2)
- Outputs :
- k cantilever spring constant (in N/m, or equivalently nN/nm)
- k_s standard deviation in our estimate of k
- !!!a2_r relative error in a**2
- !!!T_r relative error in T
- !!!one_o_Vp2_r relative error in 1/Vphotodiode_variance
- Notes :
- We're assuming vib is mostly from thermal cantilever vibrations
- (and then only from vibrations in the single vertical degree of freedom),
- and not from other noise sources.
- The various relative errors are returned to help you gauge the
- largest source of random error in your measurement of k.
- If one of them is small, don't bother repeating that measurment too often.
- If one is large, try repeating that measurement more.
- Remember that you need enough samples to have a valid error estimate in
- the first place, and that none of this addresses any systematic errors.
+def calib_analyze(bumps, temperatures, vibrations):
+ """Analyze data from `get_calibration_data()`
+
+ Inputs (all are arrays of recorded data):
+ bumps measured (V_photodiode / nm_tip) proportionality constant
+ temperatures measured temperature (K)
+ vibrations measured V_photodiode variance in free solution (V**2)
+ Outputs:
+ k cantilever spring constant (in N/m, or equivalently nN/nm)
+ k_s standard deviation in our estimate of k
+
+ Notes:
+
+ We're assuming vib is mostly from thermal cantilever vibrations
+ (and then only from vibrations in the single vertical degree of
+ freedom), and not from other noise sources.
+
+ If the error is large, check the relative errors
+ (`x.std()/x.mean()`)of your input arrays. If one of them is
+ small, don't bother repeating that measurment too often. If one
+ is large, try repeating that measurement more. Remember that you
+ need enough samples to have a valid error estimate in the first
+ place, and that none of this addresses any systematic errors.
"""
- calib_plot_kwargs, = calib_analyze._splitargs(calib_analyze, kwargs)
- photoSensitivity2 = bumps**2
- one_o_Vphoto2 = 1/vibs
-
- photoSensitivity2_m = photoSensitivity2.mean()
- T_m = Ts.mean()
- one_o_Vphoto2_m = one_o_Vphoto2.mean()
- # Vphoto / photoSensitivity = x
- # k = kb T / <x**2> = kb T photoSensitiviy**2 * (1e9nm/m)**2 /
- # <Vphoto_std**2>
+ ps_m = bumps.mean() # ps for photo-sensitivity
+ ps_s = bumps.std()
+ T_m = temperatures.mean()
+ T_s = temperatures.std()
+ v2_m = vibrations.mean() # average voltage variance
+ v2_s = vibrations.std()
+
+ # Vphoto / photo_sensitivity = x
+ # k = kB T / <x**2> = kB T photo_sensitivity**2 / Vphoto_var
#
- # units, photoSensitivity = Vphoto/(Zcant in nm),
- # so Vphoto/photoSensitivity = Zcant in nm
- # so k = J/K * K / nm^2 * (1e9nm/m)**2 = N/m
- k = kb * T_m * photoSensitivity2_m * one_o_Vphoto2_m * 1e18
-
- # propogation of errors !!!
- # first, get standard deviations
- photoSensitivity2_s = photoSensitivity2.std()
- T_s = Ts.std()
- one_o_Vphoto2_s = one_o_Vphoto2.std()
- # !!!!now, get relative errors
- photoSensitivity2_r = photoSensitivity2_s / photoSensitivity2_m
- T_r = T_s / T_m
- one_o_Vphoto2_r = one_o_Vphoto2_s / one_o_Vphoto2_m
-
- k_s = k*(photoSensitivity2_r**2 + T_r**2 + one_o_Vphoto2_r**2)**0.5
-
- calib_plot(bumps, Ts, vibs, **calib_plot_kwargs)
-
- return (k, k_s,
- photoSensitivity2_m, photoSensitivity2_s,
- T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s)
-
-@splittableKwargsFunction()
-def string_errors(k_m, k_s,
- photoSensitivity2_m, photoSensitivity2_s,
- T_m, T_s,
- one_o_Vphoto2_m, one_o_Vphoto2_s) :
- k_r = k_s / k_m
- photoSensitivity2_r = photoSensitivity2_s / photoSensitivity2_m
- T_r = T_s / T_m
- one_o_Vphoto2_r = one_o_Vphoto2_s / one_o_Vphoto2_m
- string = "Variable (units) : mean +/- std. dev. (relative error)\n"
- string += "Cantilever k (N/m) : %g +/- %g (%g)\n" \
- % (k_m, k_s, k_r)
- string += "photoSensitivity**2 (V/nm)**2 : %g +/- %g (%g)\n" \
- % (photoSensitivity2_m, photoSensitivity2_s, photoSensitivity2_r)
- string += "T (K) : %g +/- %g (%g)\n" \
- % (T_m, T_s, T_r)
- string += "1/Vphoto**2 (1/V)**2 : %g +/- %g (%g)" \
- % (one_o_Vphoto2_m, one_o_Vphoto2_s, one_o_Vphoto2_r)
- return string
-
-@splittableKwargsFunction()
-def calib_save(bumps, Ts, vibs, log_dir=None) :
- """
- Save a dictonary with the bump, T, and vib data.
- """
- if log_dir != None :
- data = {'bump':bumps, 'T':Ts, 'vib':vibs}
- log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
- log_name="calib")
- log.write_dict_of_arrays(data)
-
-def calib_load(datafile) :
- "Load the dictionary data, using data_logger.date_load()"
- dl = data_logger.data_load()
- data = dl.read_dict_of_arrays(path)
- return (data['bump'], data['T'], data['vib'])
-
-def calib_save_analysis(k, k_s,
- photoSensitivity2_m, photoSensitivity2_s,
- T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s,
- log_dir=None) :
- if log_dir != None :
- log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
- log_name="calib_analysis_text")
- log.write_binary(string_errors(k, k_s,
- photoSensitivity2_m, photoSensitivity2_s,
- T_m, T_s,
- one_o_Vphoto2_m, one_o_Vphoto2_s)
- +'\n')
-
-@splittableKwargsFunction()
-def calib_plot(bumps, Ts, vibs, plotVerbose=False) :
- if plotVerbose or config.PYLAB_VERBOSE :
- common._import_pylab()
- common._pylab.figure(config.BASE_FIGNUM+4)
- common._pylab.subplot(311)
- common._pylab.plot(bumps, 'g.-')
- common._pylab.title('Photodiode sensitivity (V/nm)')
- common._pylab.subplot(312)
- common._pylab.plot(Ts, 'r.-')
- common._pylab.title('Temperature (K)')
- common._pylab.subplot(313)
- common._pylab.plot(vibs, 'b.-')
- common._pylab.title('Thermal deflection variance (Volts**2)')
- common._flush_plot()
-
-make_splittable_kwargs_function(calib_analyze,
- (calib_plot, 'bumps', 'Ts', 'vibs'))
-
-@splittableKwargsFunction((calib_analyze, 'bumps', 'Ts', 'vibs'))
-def calib_load_analyze_tweaks(bump_tweaks, vib_tweaks, T_tweaks=None) :
- raise NotImplementedError
- a = read_tweaked_bumps(bump_tweaks)
- vib = V_photo_variance_from_file(vib_tweaks)
- if T_tweaks == None:
- pass
- return analyze_calibration_data(a, T, vib, log_dir=log_dir)
-
-# commandline interface functions
-import scipy.io, sys
-
-def array_from_string(string):
- ret = []
- for num in string.split(',') :
- ret.append(float(num))
- assert len(ret) > 0
- return numpy.array(ret)
-
-def read_data(ifile):
- "ifile can be a filename string or open (seekable) file object"
- unlabeled_data=scipy.io.read_array(file)
- return unlabeled_data
-
-def get_array(string, filename, name) :
- "get an array from supplied command line options"
- if string != None :
- array = array_from_string(string)
- elif filename != None :
- array = read_data(filename)
- else :
- raise Exception, "no %s information given" % (name)
- return array
-
-if __name__ == '__main__' :
- # command line interface
- from optparse import OptionParser
-
- usage_string = ('%prog <bumps> <temps> <vibs>\n'
- '2008, W. Trevor King.\n'
- '\n'
- 'Takes arrays of Vphotodiode sensitivity (V/nm), Temperature (K), \n'
- 'and Vibration variance (V**2) as comma seperated lists.\n'
- 'Returns the cantilever spring constant (pN/nm).\n'
- 'for example:\n'
- ' $ %prog -b 0.02,0.03,0.025 -t 298.2,300.1 -v 6e-9,5.5e-9\n'
- )
- parser = OptionParser(usage=usage_string, version='%prog '+common.VERSION)
- parser.add_option('-b', '--bump-string', dest='bump_string',
- help='comma seperated photodiode sensitivities (V/nm)',
- type='string', metavar='BUMPS')
- parser.add_option('-t', '--temp-string', dest='temp_string',
- help='comma seperated temperatures (K)',
- type='string', metavar='TEMPS')
- parser.add_option('-v', '--vib-string', dest='vib_string',
- help='comma seperated vibration variances (V**2)',
- type='string', metavar='VIBS')
- parser.add_option('-B', '--bump-file', dest='bump_file',
- help='comma seperated photodiode sensitivities (V/nm)',
- type='string', metavar='BUMPFILE')
- parser.add_option('-T', '--temp-file', dest='temp_file',
- help='comma seperated temperatures (K)',
- type='string', metavar='TEMPFILE')
- parser.add_option('-V', '--vib-file', dest='vib_file',
- help='comma seperated vibration variances (V**2)',
- type='string', metavar='VIBFILE')
- parser.add_option('-C', '--celsius', dest='celsius',
- help='Use Celsius input temperatures instead of Kelvin (defaul %default)\n',
- action='store_true', default=False)
- parser.add_option('-o', '--output-file', dest='ofilename',
- help='write output to FILE (default stdout)',
- type='string', metavar='FILE')
- parser.add_option('-p', '--plot-inputs', dest='plot',
- help='plot the input arrays to check their distribution',
- action='store_true', default=False)
- parser.add_option('', '--verbose', dest='verbose', action='store_true',
- help='print lots of debugging information',
- default=False)
-
- options,args = parser.parse_args()
- parser.destroy()
-
- config.TEXT_VERBOSE = options.verbose
- config.PYLAB_INTERACTIVE = False
- config.PYLAB_VERBOSE = options.plot
-
- bumps = get_array(options.bump_string, options.bump_file, "bump")
- temps = get_array(options.temp_string, options.temp_file, "temp")
- vibs = get_array(options.vib_string, options.vib_file, "vib")
-
- #if options.plot :
- # calib_plot(bumps, temps, vibs)
-
- if options.celsius :
- for i in range(len(temps)) :
- temps[i] = T_analyze.C_to_K(temps[i])
-
- km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s = \
- calib_analyze(bumps, temps, vibs, plotVerbose=options.plot)
-
- if options.verbose :
- print >> sys.stderr, \
- string_errors(km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s)
-
- if options.ofilename != None :
- print >> file(options.ofilename, 'w'), km
- else :
- print km
+ # units, photo_sensitivity = Vphoto/(Zcant in m),
+ # so Vphoto/photo_sensitivity = Zcant in m
+ # so k = J/K * K / m^2 = J / m^2 = N/m
+ k = _kB * T_m * ps_m**2 / v2_m
+
+ # propogation of errors
+ # dk/dT = k/T
+ dk_T = k/T_m * T_s
+ # dk/dps = 2k/ps
+ dk_ps = 2*k/ps_m * ps_s
+ # dk/dv2 = -k/v2
+ dk_v = -k/v2_m * v2_s
+
+ k_s = _numpy.sqrt(dk_T**2 + dk_ps**2 + dk_v**2)
+
+ _LOG.info('variable (units) : '
+ 'mean +/- std. dev. (relative error)')
+ _LOG.info('cantilever k (N/m) : %g +/- %g (%g)' % (k, k_s, k_s/k))
+ _LOG.info('photo sensitivity (V/m) : %g +/- %g (%g)'
+ % (ps_m, ps_s, ps_s/ps_m))
+ _LOG.info('T (K) : %g +/- %g (%g)'
+ % (T_m, T_s, T_s/T_m))
+ _LOG.info('vibration variance (V^2) : %g +/- %g (%g)'
+ % (v2_m, v2_s, v2_s/v2_m))
+
+ if _base_config['matplotlib']:
+ calib_plot(bumps, temperatures, vibrations)
+
+ return (k, k_s)
+
+def calib_save(filename, group='/', bumps=None, temperatures=None,
+ vibrations=None, calibration_config=None, k=None, k_s=None):
+ with _h5py.File(filename, 'a') as f:
+ cwg = _h5_create_group(f, group)
+ if calibration_config is not None:
+ config_cwg = _h5_create_group(cwg, 'config')
+ calibration_config.save(group=config_cwg)
+ if bumps is not None:
+ try:
+ del cwg['raw/photodiode-sensitivity/data']
+ except KeyError:
+ pass
+ try:
+ del cwg['raw/photodiode-sensitivity/units']
+ except KeyError:
+ pass
+ cwg['raw/photodiode-sensitivity/data'] = bumps
+ cwg['raw/photodiode-sensitivity/units'] = 'V/m'
+ if temperatures is not None:
+ try:
+ del cwg['raw/temperature/data']
+ except KeyError:
+ pass
+ try:
+ del cwg['raw/temperature/units']
+ except KeyError:
+ pass
+ cwg['raw/temperature/data'] = temperatures
+ cwg['raw/temperature/units'] = 'K'
+ if vibrations is not None:
+ try:
+ del cwg['raw/thermal-vibration-variance/data']
+ except KeyError:
+ pass
+ try:
+ del cwg['raw/thermal-vibration-variance/units']
+ except KeyError:
+ pass
+ cwg['raw/thermal-vibration-variance/data'] = vibrations
+ cwg['raw/thermal-vibration-variance/units'] = 'V^2'
+ if k is not None:
+ try:
+ del cwg['processed/spring-constant/data']
+ except KeyError:
+ pass
+ try:
+ del cwg['processed/spring-constant/units']
+ except KeyError:
+ pass
+ cwg['processed/spring-constant/data'] = k
+ cwg['processed/spring-constant/units'] = 'N/m'
+ if k_s is not None:
+ try:
+ del cwg['processed/spring-constant/standard-deviation']
+ except KeyError:
+ pass
+ cwg['processed/spring-constant/standard-deviation'] = k_s
+
+def calib_load(filename, group='/'):
+ assert group.endswith('/')
+ bumps = temperatures = vibrations = k = k_s = None
+ configs = []
+ with _h5py.File(filename, 'a') as f:
+ try:
+ bumps = f[group+'raw/photodiode-sensitivity/data'][...]
+ except KeyError:
+ pass
+ try:
+ temperatures = f[group+'raw/temperature/data'][...]
+ except KeyError:
+ pass
+ try:
+ vibrations = f[group+'raw/thermal-vibration-variance/data'][...]
+ except KeyError:
+ pass
+ try:
+ k = f[group+'processed/spring-constant/data'][...]
+ except KeyError:
+ pass
+ try:
+ k_s = f[group+'processed/spring-constant/standard-deviation'][...]
+ except KeyError:
+ pass
+ calibration_config = _HDF5_CalibrationConfig(
+ filename=filename, group=group+'config/')
+ calibration_config.load()
+ return (bumps, temperatures, vibrations, calibration_config, k, k_s)
+
+def calib_plot(bumps, temperatures, vibrations):
+ if not _matplotlib:
+ raise _matplotlib_import_error
+ figure = _matplotlib_pyplot.figure()
+
+ bump_axes = figure.add_subplot(3, 1, 1)
+ T_axes = figure.add_subplot(3, 1, 2)
+ vib_axes = figure.add_subplot(3, 1, 3)
+
+ timestamp = _time.strftime('%H%M%S')
+ bump_axes.set_title('cantilever calibration %s' % timestamp)
+
+ bump_axes.plot(bumps, 'g.-')
+ bump_axes.set_ylabel('photodiode sensitivity (V/m)')
+ T_axes.plot(temperatures, 'r.-')
+ T_axes.set_ylabel('temperature (K)')
+ vib_axes.plot(vibrations, 'b.-')
+ vib_axes.set_ylabel('thermal deflection variance (V^2)')
+
+ figure.show()
+
+
+def calib_load_all(filename, group='/'):
+ "Load all data from a `calib()` run."
+ assert group.endswith('/'), group
+ bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load(
+ filename, group+'calibration/')
+ bump_details = []
+ for i in range(calibration_config['num-bumps']):
+ (raw_bump,bump_config,z_channel_config,z_axis_config,
+ deflection_channel_config,processed_bump) = _bump_load(
+ filename=filename, group='%sbump/%d/' % (group, i))
+ bump_details.append({
+ 'raw_bump': raw_bump,
+ 'bump_config': bump_config,
+ 'z_channel_config': z_channel_config,
+ 'z_axis_config': z_axis_config,
+ 'deflection_channel_config': deflection_channel_config,
+ 'processed_bump': processed_bump,
+ })
+ temperature_details = []
+ for i in range(calibration_config['num-temperatures']):
+ (raw_temperature,temperature_config,processed_temperature
+ ) = _temperature_load(
+ filename=filename, group='%stemperature/%d/' % (group, i))
+ temperature_details.append({
+ 'raw_temperature': raw_temperature,
+ 'temperature_config': temperature_config,
+ 'processed_temperature': processed_temperature,
+ })
+ vibration_details = []
+ for i in range(calibration_config['num-vibrations']):
+ (raw_vibration,vibration_config,deflection_channel_config,
+ processed_vibration) = _vibration_load(
+ filename=filename, group='%svibration/%d/' % (group, i))
+ vibration_details.append({
+ 'raw_vibration': raw_vibration,
+ 'vibration_config': vibration_config,
+ 'deflection_channel_config': deflection_channel_config,
+ 'processed_vibration': processed_vibration,
+ })
+ return {
+ 'bumps': bumps,
+ 'bump_details': bump_details,
+ 'temperatures': temperatures,
+ 'temperature_details': temperature_details,
+ 'vibrations': vibrations,
+ 'vibration_details': vibration_details,
+ 'calibration_config': calibration_config,
+ 'k': k,
+ 'k_s': k_s,
+ }
+
+def calib_analyze_all(filename, group='/', maximum_relative_error=1e-5,
+ dry_run=False):
+ "(Re)analyze all data from a `calib()` run."
+ assert group.endswith('/'), group
+ bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load(
+ filename, group+'calibration/')
+ changed_bump = changed_temperature = changed_vibration = False
+ for i in range(calibration_config['num-bumps']):
+ bump_group = '%sbump/%d/' % (group, i)
+ (raw_bump,bump_config,z_channel_config,z_axis_config,
+ deflection_channel_config,processed_bump) = _bump_load(
+ filename=filename, group=bump_group)
+ sensitivity = _bump_analyze(
+ data=raw_bump, bump_config=bump_config,
+ z_channel_config=z_channel_config,
+ z_axis_config=z_axis_config,
+ deflection_channel_config=deflection_channel_config)
+ bumps[i] = sensitivity
+ rel_error = abs(sensitivity - processed_bump)/processed_bump
+ if rel_error > maximum_relative_error:
+ _LOG.warn(("new analysis doesn't match for bump %d: %g -> %g "
+ "(difference: %g, relative error: %g)")
+ % (i, processed_bump, sensitivity,
+ sensitivity-processed_bump, rel_error))
+ changed_bump = True
+ if not dry_run:
+ _bump_save(filename, bump_group, processed_bump=sensitivity)
+ for i in range(calibration_config['num-temperatures']):
+ temperature_group = '%stemperature/%d/' % (group, i)
+ (raw_temperature,temperature_config,processed_temperature
+ ) = _temperature_load(
+ filename=filename, group=temperature_group)
+ temperature = _temperature_analyze(
+ raw_temperature, temperature_config)
+ temperatures[i] = temperature
+ rel_error = abs(temperature - processed_temperature
+ )/processed_temperature
+ if rel_error > maximum_relative_error:
+ _LOG.warn(("new analysis doesn't match for temperature %d: "
+ "%g -> %g (difference: %g, relative error: %g)")
+ % (i, processed_temperature, temperature,
+ temperature-processed_temperature, rel_error))
+ changed_temperature = True
+ if not dry_run:
+ _temperature_save(
+ filename, temperature_group,
+ processed_T=temperature)
+ for i in range(calibration_config['num-vibrations']):
+ vibration_group = '%svibration/%d/' % (group, i)
+ (raw_vibration,vibration_config,deflection_channel_config,
+ processed_vibration) = _vibration_load(
+ filename=filename, group=vibration_group)
+ variance = _vibration_analyze(
+ deflection=raw_vibration, vibration_config=vibration_config,
+ deflection_channel_config=deflection_channel_config)
+ vibrations[i] = variance
+ rel_error = abs(variance - processed_vibration)/processed_vibration
+ if rel_error > maximum_relative_error:
+ _LOG.warn(("new analysis doesn't match for vibration %d: %g -> %g "
+ "(difference: %g, relative error: %g)")
+ % (i, processed_vibration, variance,
+ variance-processed_vibration, rel_error))
+ changed_vibration = True
+ if not dry_run:
+ _vibration_save(
+ filename, vibration_group, processed_vibration=variance)
+
+ calib_group = '%scalibration/' % group
+
+ if changed_bump and not dry_run:
+ calib_save(filename, calib_group, bumps=bumps)
+ if changed_temperature and not dry_run:
+ calib_save(filename, calib_group, temperatures=temperatures)
+ if changed_vibration and not dry_run:
+ calib_save(filename, calib_group, vibrations=vibrations)
+
+ new_k,new_k_s = calib_analyze(
+ bumps=bumps, temperatures=temperatures, vibrations=vibrations)
+ rel_error = abs(new_k-k)/k
+ if rel_error > maximum_relative_error:
+ _LOG.warn(("new analysis doesn't match for k: %g -> %g "
+ "(difference: %g, relative error: %g)")
+ % (k, new_k, new_k-k, rel_error))
+ if not dry_run:
+ calib_save(filename, calib_group, k=new_k)
+ rel_error = abs(new_k_s-k_s)/k_s
+ if rel_error > maximum_relative_error:
+ _LOG.warn(("new analysis doesn't match for k_s: %g -> %g "
+ "(difference: %g, relative error: %g)")
+ % (k_s, new_k_s, new_k_s-k_s, rel_error))
+ if not dry_run:
+ calib_save(filename, calib_group, k_s=new_k_s)
+ return (new_k, new_k_s)
+
+def calib_plot_all(bumps, bump_details, temperatures, temperature_details,
+ vibrations, vibration_details, calibration_config, k, k_s,
+ maximum_relative_error=1e-5):
+ calib_plot(bumps, temperatures, vibrations)
+ for i,bump in enumerate(bump_details):
+ sensitivity = _bump_analyze(
+ data=bump['raw_bump'], bump_config=bump['bump_config'],
+ z_channel_config=bump['z_channel_config'],
+ z_axis_config=bump['z_axis_config'],
+ deflection_channel_config=bump['deflection_channel_config'],
+ plot=True)
+ rel_error = abs(sensitivity - bump['processed_bump']
+ )/bump['processed_bump']
+ if rel_error > maximum_relative_error:
+ _LOG.warn(("new analysis doesn't match for bump %d: %g != %g "
+ "(difference: %g, relative error: %g)")
+ % (i, sensitivity, bump['processed_bump'],
+ sensitivity-bump['processed_bump'], rel_error))
+ # nothing interesting to plot for temperatures...
+ for i,vibration in enumerate(vibration_details):
+ variance = _vibration_analyze(
+ deflection=vibration['raw_vibration'],
+ vibration_config=vibration['vibration_config'],
+ deflection_channel_config=vibration['deflection_channel_config'],
+ plot=True)
+ rel_error = abs(variance - vibration['processed_vibration']
+ )/vibration['processed_vibration']
+ if rel_error > maximum_relative_error:
+ _LOG.warn(("new analysis doesn't match for vibration %d: %g != %g "
+ "(difference: %g, relative error: %g)")
+ % (i, variance, vibration['processed_vibration'],
+ variance-vibration['processed_vibration'], rel_error))
# License along with calibcant. If not, see
# <http://www.gnu.org/licenses/>.
-"""
-Aquire, save, and load cantilever calibration bump data.
+"""Acquire, save, and load cantilever calibration bump data.
+
For measuring photodiode sensitivity.
-W. Trevor King Dec. 2007 - Oct. 2008
+The relevent physical quantities are:
+ Vzp_out Output z-piezo voltage (what we generate)
+ Vzp Applied z-piezo voltage (after external ZPGAIN)
+ Zp The z-piezo position
+ Zcant The cantilever vertical deflection
+ Vphoto The photodiode vertical deflection voltage (what we measure)
-The relevent physical quantities are :
- Vzp_out Output z-piezo voltage (what we generate)
- Vzp Applied z-piezo voltage (after external ZPGAIN)
- Zp The z-piezo position
- Zcant The cantilever vertical deflection
- Vphoto The photodiode vertical deflection voltage (what we measure)
+Which are related by the parameters:
+ zp_gain Vzp_out / Vzp
+ zp_sensitivity Zp / Vzp
+ photo_sensitivity Vphoto / Zcant
-Which are related by the parameters :
- zpGain Vzp_out / Vzp
- zpSensitivity Zp / Vzp
- photoSensitivity Vphoto / Zcant
+Cantilever calibration assumes a pre-calibrated z-piezo
+(zp_sensitivity) and amplifier (zp_gain). In our lab, the z-piezo is
+calibrated by imaging a calibration sample, which has features with
+well defined sizes, and the gain is set with a knob on our modified
+NanoScope.
-Cantilever calibration assumes a pre-calibrated z-piezo (zpSensitivity) and
-amplifier (zpGain). In our lab, the z-piezo is calibrated by imaging a
-calibration sample, which has features with well defined sizes, and the gain
-is set with a knob on the Nanoscope.
+Photo-sensitivity is measured by bumping the cantilever against the
+surface, where `Zp = Zcant` (see the `bump_*()` family of functions).
+The measured slope Vphoto/Vout is converted to photo-sensitivity via
-photoSensitivity is measured by bumping the cantilever against the surface,
-where Zp = Zcant (see the bump_*() family of functions)
-The measured slope Vphoto/Vout is converted to photoSensitivity via
-Vphoto/Vzp_out * Vzp_out/Vzp * Vzp/Zp * Zp/Zcant = Vphoto/Zcant
- (measured) (1/zpGain) (1/zpSensitivity) (1) (photoSensitivity)
+ Vphoto/Vzp_out * Vzp_out/Vzp * Vzp/Zp * Zp/Zcant = Vphoto/Zcant
+ (measured) (1/zp_gain) (1/zp_sensitivity) (1) (photo_sensitivity)
-We do all these measurements a few times to estimate statistical errors.
+We do all these measurements a few times to estimate statistical
+errors.
The functions are layed out in the families:
- bump_*()
-For each family, * can be any of :
- aquire get real-world data
- save store real-world data to disk
- load get real-world data from disk
- analyze interperate the real-world data.
- plot show a nice graphic to convince people we're working :p
- load_analyze_tweaked
- read a file with a list of paths to previously saved real world data
- load each file using *_load(), analyze using *_analyze(), and
- optionally plot using *_plot().
- Intended for re-processing old data.
-A family name without any _* extension (e.g. bump()),
- runs *_aquire(), *_save(), *_analyze(), *_plot().
+ bump_*()
+For each family, * can be any of:
+ acquire get real-world data
+ save store real-world data to disk
+ load get real-world data from disk
+ analyze interperate the real-world data.
+ plot show a nice graphic to convince people we're working :p
+
+A family name without any _* extension (e.g. `bump()`), runs `*_acquire()`,
+`*_save()`, `*_analyze()`.
+
+If `base_config['matplotlib']` is `True`, `*_analyze()` will call
+`*_plot()` internally.
"""
-import numpy
-import time
+import numpy as _numpy
-import data_logger
-import FFT_tools
-import piezo.z_piezo_utils as z_piezo_utils
+from pypiezo.base import convert_meters_to_bits as _convert_meters_to_bits
+from pypiezo.base import convert_bits_to_meters as _convert_bits_to_meters
-from .bump_analyze import bump_analyze
+from . import LOG as _LOG
+from .bump_analyze import bump_analyze as _bump_analyze
+from .bump_analyze import bump_save as _bump_save
-LOG_DATA = True # quietly grab all real-world data and log to LOG_DIR
-LOG_DIR = '${DEFAULT}/calibrate_cantilever'
+def bump_acquire(afm, bump_config):
+ """Ramps `push_depth` closer and returns to the original position.
-TEXT_VERBOSE = True # for debugging
+ Inputs:
+ afm a pyafm.AFM instance
+ bump_config a .config._BumpConfig instance
+ Returns the acquired ramp data dictionary, with data in DAC/ADC bits.
+ """
+ afm.move_just_onto_surface(
+ depth=bump_config['initial-position'], far=bump_config['far-steps'])
+
+ _LOG.info('bump the surface to a depth of %g m'
+ % bump_config['push-depth'])
+
+ axis = afm.piezo.axis_by_name(afm.axis_name)
+
+ start_pos = afm.piezo.last_output[afm.axis_name]
+ start_pos_m = _convert_bits_to_meters(
+ axis.axis_channel_config, axis.axis_config, start_pos)
+ close_pos_m = start_pos_m + bump_config['push-depth']
+ close_pos = _convert_meters_to_bits(
+ axis.axis_channel_config, axis.axis_config, close_pos_m)
+
+ dtype = afm.piezo.channel_dtype(afm.axis_name, direction='output')
+ appr = _numpy.linspace(
+ start_pos, close_pos, bump_config['samples']).astype(dtype)
+ # switch numpy.append to numpy.concatenate with version 2.0+
+ out = _numpy.append(appr, appr[::-1])
+ out = out.reshape((len(out), 1))
+
+ # (samples) / (meters) * (meters/second) = (samples/second)
+ freq = (bump_config['samples'] / bump_config['push-depth']
+ * bump_config['push-speed'])
+
+ data = afm.piezo.ramp(out, freq, output_names=[afm.axis_name],
+ input_names=['deflection'])
+
+ out = out.reshape((len(out),))
+ data = data.reshape((data.size,))
+ return {afm.axis_name: out, 'deflection': data}
+
+def bump(afm, bump_config, filename, group='/'):
+ """Wrapper around bump_acquire(), bump_analyze(), bump_save().
+
+ >>> import os
+ >>> import tempfile
+ >>> from pycomedi.device import Device
+ >>> from pycomedi.subdevice import StreamingSubdevice
+ >>> from pycomedi.channel import AnalogChannel, DigitalChannel
+ >>> from pycomedi.constant import AREF, IO_DIRECTION, SUBDEVICE_TYPE, UNIT
+ >>> from pypiezo.afm import AFMPiezo
+ >>> from pypiezo.base import PiezoAxis, InputChannel
+ >>> from pypiezo.config import (HDF5_ChannelConfig, HDF5_AxisConfig,
+ ... pprint_HDF5)
+ >>> from stepper import Stepper
+ >>> from pyafm import AFM
+ >>> from .config import HDF5_BumpConfig
+
+ >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
+ >>> os.close(fd)
+
+ >>> d = Device('/dev/comedi0')
+ >>> d.open()
+
+ Setup an `AFMPiezo` instance.
+
+ >>> s_in = d.find_subdevice_by_type(SUBDEVICE_TYPE.ai,
+ ... factory=StreamingSubdevice)
+ >>> s_out = d.find_subdevice_by_type(SUBDEVICE_TYPE.ao,
+ ... factory=StreamingSubdevice)
+
+ >>> axis_channel = s_out.channel(
+ ... 0, factory=AnalogChannel, aref=AREF.ground)
+ >>> input_channel = s_in.channel(0, factory=AnalogChannel, aref=AREF.diff)
+ >>> for chan in [axis_channel, input_channel]:
+ ... chan.range = chan.find_range(unit=UNIT.volt, min=-10, max=10)
+
+ We set the minimum voltage for the `z` axis to -9 (a volt above
+ the minimum possible voltage) to help with testing
+ `.get_surface_position`. Without this minimum voltage, small
+ calibration errors could lead to a railed -10 V input for the
+ first few surface approaching steps, which could lead to an
+ `EdgeKink` error instead of a `FlatFit` error.
+
+ >>> axis_config = HDF5_AxisConfig(filename, '/bump/config/z/axis')
+ >>> axis_config.update(
+ ... {'gain':20, 'sensitivity':8e-9, 'minimum':-9})
+ >>> axis_channel_config = HDF5_ChannelConfig(
+ ... filename, '/bump/config/z/channel')
+ >>> input_channel_config = HDF5_ChannelConfig(
+ ... filename, '/bump/config/deflection/channel')
+
+ >>> a = PiezoAxis(axis_config=axis_config,
+ ... axis_channel_config=axis_channel_config,
+ ... axis_channel=axis_channel, name='z')
+ >>> a.setup_config()
+
+ >>> c = InputChannel(
+ ... channel_config=input_channel_config, channel=input_channel,
+ ... name='deflection')
+ >>> c.setup_config()
+
+ >>> piezo = AFMPiezo(axes=[a], input_channels=[c])
+
+ Setup a `stepper` instance.
+
+ >>> s_d = d.find_subdevice_by_type(SUBDEVICE_TYPE.dio)
+ >>> d_channels = [s_d.channel(i, factory=DigitalChannel)
+ ... for i in (0, 1, 2, 3)]
+ >>> for chan in d_channels:
+ ... chan.dio_config(IO_DIRECTION.output)
+
+ >>> def write(value):
+ ... s_d.dio_bitfield(bits=value, write_mask=2**4-1)
+
+ >>> stepper = Stepper(write=write)
+
+ Setup an `AFM` instance.
+
+ >>> afm = AFM(piezo, stepper)
+
+ Test a bump:
-# bump family
+ >>> bump_config = HDF5_BumpConfig(
+ ... filename=filename, group='/bump/config/bump')
+ >>> bump(afm, bump_config, filename, group='/bump')
+ TODO: replace skipped example data with real-world values
+ >>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
-def bump_aquire(zpiezo, push_depth, npoints, freq) :
- """
- Ramps closer push_depth and returns to the original position.
- Inputs:
- zpiezo an opened zpiezo.zpiezo instance
- push_depth distance to approach, in nm
- npoints number points during the approach and during the retreat
- freq rate at which data is aquired
- log_dir directory to log data to (see data_logger.py).
- None to turn off logging (see also the global LOG_DATA).
- Returns the aquired ramp data dictionary, with data in DAC/ADC bits.
- """
- # generate the bump output
- start_pos = zpiezo.curPos()
- pos_dist = zpiezo.pos_nm2out(push_depth) - zpiezo.pos_nm2out(0)
- close_pos = start_pos + pos_dist
- appr = linspace(start_pos, close_pos, npoints)
- retr = linspace(close_pos, start_pos, npoints)
- out = concatenate((appr, retr))
- # run the bump, and measure deflection
- if TEXT_VERBOSE :
- print "Bump %g nm" % push_depth
- data = zpiezo.ramp(out, freq)
- # default saving, so we have a log in-case the operator is lazy ;)
- if LOG_DATA == True :
- log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
- log_name="bump_surface")
- log.write_dict_of_arrays(data)
- return data
-
-def bump_save(data, log_dir) :
- "Save the dictionary data, using data_logger.data_log()"
- if log_dir != None :
- log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
- log_name="bump")
- log.write_dict_of_arrays(data)
-
-def bump_load(datafile) :
- "Load the dictionary data, using data_logger.date_load()"
- dl = data_logger.data_load()
- data = dl.read_dict_of_arrays(path)
- return data
-
-def bump_plot(data, plotVerbose) :
- "Plot the bump (Vphoto vs Vzp) if plotVerbose or PYLAB_VERBOSE == True"
- if plotVerbose or PYLAB_VERBOSE :
- _import_pylab()
- _pylab.figure(BASE_FIGNUM)
- _pylab.plot(data["Z piezo output"], data["Deflection input"],
- '.', label='bump')
- _pylab.title("bump surface")
- _pylab.legend(loc='upper left')
- _flush_plot()
-
-def bump(zpiezo, push_depth, npoints=1024, freq=100e3,
- log_dir=None,
- plotVerbose=False) :
- """
- Wrapper around bump_aquire(), bump_save(), bump_analyze(), bump_plot()
- """
- data = bump_aquire(zpiezo, push_depth, npoints, freq)
- bump_save(data, log_dir)
- photoSensitivity = bump_analyze(data, zpiezo.gain, zpiezo.sensitivity,
- zpiezo.pos_out2V, zpiezo.def_in2V)
- bump_plot(data, plotVerbose)
- return photoSensitivity
-
-def bump_load_analyze_tweaked(tweak_file, zpGain=_usual_zpGain,
- zpSensitivity=_usual_zpSensitivity,
- Vzp_out2V=_usual_Vzp_out2V,
- Vphoto_in2V=_usual_Vphoto_in2V,
- plotVerbose=False) :
- "Load the output file of tweak_calib_bump.sh, return an array of slopes"
- photoSensitivity = []
- for line in file(tweak_file, 'r') :
- parsed = line.split()
- path = parsed[0].split('\n')[0]
- # read the data
- full_data = bump_load(path)
- if len(parsed) == 1 :
- data = full_data # use whole bump
- else :
- # use the listed sections
- zp = []
- df = []
- for rng in parsed[1:] :
- p = rng.split(':')
- starti = int(p[0])
- stopi = int(p[1])
- zp.extend(full_data['Z piezo output'][starti:stopi])
- df.extend(full_data['Deflection input'][starti:stopi])
- data = {'Z piezo output': array(zp),
- 'Deflection input':array(df)}
- pSi = bump_analyze(data, zpGain, zpSensitivity,
- Vzp_out2V, Vphoto_in2V, plotVerbose)
- photoSensitivity.append(pSi)
- bump_plot(data, plotVervose)
- return array(photoSensitivity, dtype=numpy.float)
+ Close the Comedi device.
+ >>> d.close()
+
+ Cleanup our temporary config file.
+
+ >>> os.remove(filename)
+ """
+ deflection_channel = afm.piezo.input_channel_by_name('deflection')
+ axis = afm.piezo.axis_by_name(afm.axis_name)
+
+ data = bump_acquire(afm, bump_config)
+ photo_sensitivity = _bump_analyze(
+ data, bump_config, z_channel_config=axis.axis_channel_config,
+ z_axis_config=axis.axis_config,
+ deflection_channel_config=deflection_channel.channel_config)
+ _bump_save(
+ filename, group, data, bump_config,
+ z_channel_config=axis.axis_channel_config,
+ z_axis_config=axis.axis_config,
+ deflection_channel_config=deflection_channel.channel_config,
+ processed_bump=photo_sensitivity)
+ return photo_sensitivity
# License along with calibcant. If not, see
# <http://www.gnu.org/licenses/>.
-"""
-Separate the more general bump_analyze() from the other bump_*()
-functions in calibcant. Also provide a command line interface
-for analyzing data acquired through other workflows.
-
-The relevant physical quantities are :
- Vzp_out Output z-piezo voltage (what we generate)
- Vzp Applied z-piezo voltage (after external ZPGAIN)
- Zp The z-piezo position
- Zcant The cantilever vertical deflection
- Vphoto The photodiode vertical deflection voltage (what we measure)
-
-Which are related by the parameters :
- zpGain Vzp_out / Vzp
- zpSensitivity Zp / Vzp
- photoSensitivity Vphoto / Zcant
-
-photoSensitivity is measured by bumping the cantilever against the
-surface, where Zp = Zcant (see calibrate.bump_aquire()). The measured
-slope Vphoto/Vout is converted to photoSensitivity with bump_analyze().
+"""Surface bump analysis (measures photodiode sensitivity).
+
+Separate the more general `bump_analyze()` from the other `bump_*()`
+functions in calibcant.
+
+The relevant physical quantities are:
+
+* `Vzp_out` Output z-piezo voltage (what we generate)
+* `Vzp` Applied z-piezo voltage (after external amplification)
+* `Zp` The z-piezo position
+* `Zcant` The cantilever vertical deflection
+* `Vphoto` The photodiode vertical deflection voltage (what we measure)
+
+Which are related by the parameters:
+
+* `zp_gain` Vzp_out / Vzp
+* `zp_sensitivity` Zp / Vzp
+* `photo_sensitivity` Vphoto / Zcant
+
+`photo_sensitivity` is measured by bumping the cantilever against the
+surface, where `Zp = Zcant` (see `calibrate.bump_acquire()`). The
+measured slope `Vphoto/Vout` is converted to `photo_sensitivity` with
+`bump_analyze()`.
+
+>>> import os
+>>> from pprint import pprint
+>>> import tempfile
+>>> import numpy
+>>> from .config import HDF5_BumpConfig
+>>> from pypiezo.config import pprint_HDF5, HDF5_ChannelConfig, HDF5_AxisConfig
+
+>>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
+>>> os.close(fd)
+
+>>> bump_config = HDF5_BumpConfig(filename=filename, group='/bump/config/')
+>>> bump_config.save()
+>>> z_channel_config = HDF5_ChannelConfig(filename=None)
+>>> z_channel_config['maxdata'] = 200
+>>> z_channel_config['conversion-coefficients'] = (0,1)
+>>> z_channel_config['conversion-origin'] = 0
+>>> z_axis_config = HDF5_AxisConfig(filename=None)
+>>> deflection_channel_config = HDF5_ChannelConfig(filename=None)
+>>> deflection_channel_config['maxdata'] = 200
+>>> deflection_channel_config['conversion-coefficients'] = (0,1)
+>>> deflection_channel_config['conversion-origin'] = 0
+
+>>> raw_bump = {
+... 'z': numpy.arange(100, dtype=numpy.uint16),
+... 'deflection': numpy.arange(100, dtype=numpy.uint16),
+... }
+>>> raw_bump['deflection'][:50] = 50
+>>> processed_bump = bump_analyze(
+... raw_bump, bump_config, z_channel_config, z_axis_config,
+... deflection_channel_config)
+>>> bump_plot(data=raw_bump) # TODO: convert to V and m
+>>> bump_save(filename=filename, group='/bump/', raw_bump=raw_bump,
+... z_channel_config=z_channel_config, z_axis_config=z_axis_config,
+... deflection_channel_config=deflection_channel_config,
+... processed_bump=processed_bump)
+
+>>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
+/
+ /bump
+ /bump/config
+ /bump/config/deflection
+ /bump/config/deflection/channel
+ <HDF5 dataset "channel": shape (), type "|S1">
+ 0
+ <HDF5 dataset "conversion-coefficients": shape (), type "|S4">
+ 0, 1
+ <HDF5 dataset "conversion-origin": shape (), type "|S1">
+ 0
+ <HDF5 dataset "device": shape (), type "|S12">
+ /dev/comedi0
+ <HDF5 dataset "inverse-conversion-coefficients": shape (), type "|S1">
+<BLANKLINE>
+ <HDF5 dataset "inverse-conversion-origin": shape (), type "|S3">
+ 1.0
+ <HDF5 dataset "maxdata": shape (), type "|S3">
+ 200
+ <HDF5 dataset "range": shape (), type "|S1">
+ 1
+ <HDF5 dataset "subdevice": shape (), type "|S2">
+ -1
+ <HDF5 dataset "far-steps": shape (), type "|S3">
+ 200
+ <HDF5 dataset "initial-position": shape (), type "|S6">
+ -5e-08
+ <HDF5 dataset "model": shape (), type "|S9">
+ quadratic
+ <HDF5 dataset "push-depth": shape (), type "|S5">
+ 2e-07
+ <HDF5 dataset "push-speed": shape (), type "|S5">
+ 1e-06
+ <HDF5 dataset "samples": shape (), type "|S4">
+ 1024
+ <HDF5 dataset "setpoint": shape (), type "|S3">
+ 2.0
+ /bump/config/z
+ /bump/config/z/axis
+ <HDF5 dataset "gain": shape (), type "|S3">
+ 1.0
+ <HDF5 dataset "maximum": shape (), type "|S4">
+ None
+ <HDF5 dataset "minimum": shape (), type "|S4">
+ None
+ <HDF5 dataset "sensitivity": shape (), type "|S3">
+ 1.0
+ /bump/config/z/channel
+ <HDF5 dataset "channel": shape (), type "|S1">
+ 0
+ <HDF5 dataset "conversion-coefficients": shape (), type "|S4">
+ 0, 1
+ <HDF5 dataset "conversion-origin": shape (), type "|S1">
+ 0
+ <HDF5 dataset "device": shape (), type "|S12">
+ /dev/comedi0
+ <HDF5 dataset "inverse-conversion-coefficients": shape (), type "|S1">
+<BLANKLINE>
+ <HDF5 dataset "inverse-conversion-origin": shape (), type "|S3">
+ 1.0
+ <HDF5 dataset "maxdata": shape (), type "|S3">
+ 200
+ <HDF5 dataset "range": shape (), type "|S1">
+ 1
+ <HDF5 dataset "subdevice": shape (), type "|S2">
+ -1
+ <HDF5 dataset "processed": shape (), type "<f8">
+ 1.00000002115
+ /bump/raw
+ <HDF5 dataset "deflection": shape (100,), type "<u2">
+ [50 50 ... 50 51 52 ... 97 98 99]
+ <HDF5 dataset "z": shape (100,), type "<u2">
+ [ 0 1 2 3 ... 97 98 99]
+
+>>> (raw_bump,bump_config,z_channel_config,z_axis_config,
+... deflection_channel_config,processed_bump) = bump_load(
+... filename=filename, group='/bump/')
+
+>>> pprint(raw_bump) # doctest: +ELLIPSIS
+{'deflection': array([50, 50, ... 51, 52, 53, ..., 97, 98, 99], dtype=uint16),
+ 'z': array([ 0, 1, 2, ..., 97, 98, 99], dtype=uint16)}
+>>> processed_bump # doctest: +ELLIPSIS
+1.00...
+
+>>> os.remove(filename)
"""
-import numpy
-import scipy.optimize
-
-import data_logger
-from splittable_kwargs import splittableKwargsFunction, \
- make_splittable_kwargs_function
-
-from . import common
-from . import config
-
-
-@splittableKwargsFunction()
-def Vzp_bits2nm(data_bits, zpGain=config.zpGain,
- zpSensitivity=config.zpSensitivity,
- Vzp_out2V=config.Vzp_out2V):
- scale_Vzp_bits2V = Vzp_out2V(1) - Vzp_out2V(0)
- data_V = data_bits / scale_Vzp_bits2V
- # bits / (bits/V) = V
- data_nm = data_V * zpGain * zpSensitivity
- return data_nm
-
-@splittableKwargsFunction()
-def Vphoto_bits2V(data_bits, Vphoto_in2V=config.Vphoto_in2V):
- scale_Vphoto_bits2V = Vphoto_in2V(1) - Vphoto_in2V(0)
- Vphoto_V = data_bits / scale_Vphoto_bits2V
- # bits / (bits/V) = V
- return Vphoto_V
-
-@splittableKwargsFunction((Vzp_bits2nm, 'data_bits'),
- (Vphoto_bits2V, 'data_bits'))
-def slope_bitspbit2Vpnm(slope_bitspbit, **kwargs):
- zp_kwargs,photo_kwargs = slope_bitspbit2Vpnm._splitargs(slope_bitspbit2Vpnm, kwargs)
- Vzp_bits = 1.0
- Vphoto_bits = slope_bitspbit * Vzp_bits
- return Vphoto_bits2V(Vphoto_bits, **photo_kwargs)/Vzp_bits2nm(Vzp_bits, **zp_kwargs)
-
-#@splittableKwargsFunction((bump_fit, 'zpiezo_output_bits',
-# 'deflection_input_bits'),
-# (slope_bitspbit2Vpnm, 'slope_bitspbit'))
-# Some of the child functions aren't yet defined, so postpone
-# make-splittable until later in the module.
-def bump_analyze(data, **kwargs) :
- """
- Return the slope of the bump ;).
+import h5py as _h5py
+import numpy as _numpy
+from scipy.optimize import leastsq as _leastsq
+
+try:
+ import matplotlib as _matplotlib
+ import matplotlib.pyplot as _matplotlib_pyplot
+ import time as _time # for timestamping lines on plots
+except (ImportError, RuntimeError), e:
+ _matplotlib = None
+ _matplotlib_import_error = e
+
+from pypiezo.base import convert_bits_to_volts as _convert_bits_to_volts
+from pypiezo.base import convert_bits_to_meters as _convert_bits_to_meters
+from pypiezo.config import HDF5_ChannelConfig as _HDF5_ChannelConfig
+from pypiezo.config import HDF5_AxisConfig as _HDF5_AxisConfig
+from pypiezo.config import h5_create_group as _h5_create_group
+
+from . import LOG as _LOG
+from . import base_config as _base_config
+from .config import Linear as _Linear
+from .config import Quadratic as _Quadratic
+from .config import HDF5_BumpConfig as _HDF5_BumpConfig
+
+
+def bump_analyze(data, bump_config, z_channel_config, z_axis_config,
+ deflection_channel_config, plot=False):
+ """Return the slope of the bump.
+
Inputs:
- data dictionary of data in DAC/ADC bits
- Vzp_out2V function that converts output DAC bits to Volts
- Vphoto_in2V function that converts input ADC bits to Volts
- zpGain zpiezo applied voltage per output Volt
- zpSensitivity nm zpiezo response per applied Volt
+ data dictionary of data in DAC/ADC bits
+ bump_config `.config._BumpConfig` instance
+ z_channel_config z `pypiezo.config._ChannelConfig` instance
+ z_axis_config z `pypiezo.config._AxisConfig` instance
+ deflection_channel_config
+ deflection `pypiezo.config._ChannelConfig` instance
+ plot boolean overriding matplotlib config setting.
Returns:
- photoSensitivity (Vphoto/Zcant) in Volts/nm
- Checks for strong correlation (r-value) and low randomness chance (p-value)
-
- With the current implementation, the data is regressed in DAC/ADC bits
- and THEN converted, so we're assuming that both conversions are LINEAR.
- If they aren't, rewrite to convert before the regression.
+ photo_sensitivity (Vphoto/Zcant) in Volts/m
+
+ Checks for strong correlation (r-value) and low randomness chance
+ (p-value).
"""
- bump_fit_kwargs,slope_bitspbit2Vpnm_kwargs = \
- bump_analyze._splitargs(bump_analyze, kwargs)
- Vphoto2Vzp_out_bit = bump_fit(data['Z piezo output'],
- data['Deflection input'],
- **bump_fit_kwargs)
- return slope_bitspbit2Vpnm(Vphoto2Vzp_out_bit, **slope_bitspbit2Vpnm_kwargs)
+ z = _convert_bits_to_meters(z_channel_config, z_axis_config, data['z'])
+ deflection = _convert_bits_to_volts(
+ deflection_channel_config, data['deflection'])
+ if bump_config['model'] == _Linear:
+ kwargs = {
+ 'param_guesser': limited_linear_param_guess,
+ 'model': limited_linear,
+ 'sensitivity_from_fit_params': limited_linear_sensitivity,
+ }
+ else: # _Quadratic
+ kwargs = {
+ 'param_guesser': limited_quadratic_param_guess,
+ 'model': limited_quadratic,
+ 'sensitivity_from_fit_params': limited_quadratic_sensitivity,
+ }
+ photo_sensitivity = bump_fit(z, deflection, plot=plot, **kwargs)
+ return photo_sensitivity
def limited_linear(x, params):
"""
high_voltage_rail = 2**16 - 1 # bits
x_contact,y_contact,slope = params
y = slope*(x-x_contact) + y_contact
- y = numpy.clip(y, y_contact, high_voltage_rail)
+ y = _numpy.clip(y, y_contact, high_voltage_rail)
return y
-def limited_linear_param_guess(x, y) :
+def limited_linear_param_guess(x, y):
"""
Guess rough parameters for a limited_linear model. Assumes the
bump approaches (raising the deflection as it does so) first.
i = 0
y_low = y_contact + 0.3 * (y_max-y_contact)
y_high = y_contact + 0.7 * (y_max-y_contact)
- while y[i] < y_low :
+ while y[i] < y_low:
i += 1
i_low = i
- while y[i] < y_high :
+ while y[i] < y_high:
i += 1
i_high = i
x_contact = float(x[i_low])
high_voltage_rail = 2**16 - 1 # bits
x_contact,y_contact,slope,quad = params
y = slope*(x-x_contact) + quad*(x-x_contact)**2+ y_contact
- y = numpy.clip(y, y_contact, high_voltage_rail)
+ y = _numpy.clip(y, y_contact, high_voltage_rail)
return y
-def limited_quadratic_param_guess(x, y) :
+def limited_quadratic_param_guess(x, y):
"""
Guess rough parameters for a limited_quadratic model. Assumes the
bump approaches (raising the deflection as it does so) first.
of thresholds 0.3 and 0.7 of the y value's total range. Not the
most efficient algorithm, but it seems fairly robust.
"""
- x_contact,y_contact,slope = limited_linear_param_guess(x,y)
- quad = 0
+ x_contact,y_contact,linear_slope = limited_linear_param_guess(x,y)
+ contact_depth = x.max() - x_contact
+ slope = linear_slope / 2
+ quad = slope / contact_depth
+ # The above slope and quad were chosen so
+ # x = contact_depth
+ # x*slope + x**2 * slope == x * linear_slope
return (x_contact, y_contact, slope, quad)
def limited_quadratic_sensitivity(params):
slope = params[2]
return slope
-@splittableKwargsFunction()
-def bump_fit(zpiezo_output_bits, deflection_input_bits,
+def bump_fit(z, deflection,
param_guesser=limited_quadratic_param_guess,
model=limited_quadratic,
sensitivity_from_fit_params=limited_quadratic_sensitivity,
- plotVerbose=False) :
- x = zpiezo_output_bits
- y = deflection_input_bits
- def residual(p, y, x) :
- return model(x, p) - y
- param_guess = param_guesser(x, y)
- p,cov,info,mesg,ier = \
- scipy.optimize.leastsq(residual, param_guess, args=(y, x),
- full_output=True, maxfev=int(10e3))
- if config.TEXT_VERBOSE :
- print "Fitted params:",p
- print "Covariance mx:",cov
- print "Info:", info
- print "mesg:", mesg
- if ier == 1 :
- print "Solution converged"
- else :
- print "Solution did not converge"
- if plotVerbose or config.PYLAB_VERBOSE :
- yguess = model(x, param_guess)
- #yguess = None # Don't print the guess, since I'm convinced it's ok ;).
- yfit = model(x, p)
- bump_plot(data={"Z piezo output":x, "Deflection input":y},
- yguess=yguess, yfit=yfit, plotVerbose=plotVerbose)
+ plot=False):
+ """Fit a aurface bump and return the photodiode sensitivitiy.
+
+ Parameters:
+ z piezo position in meters
+ deflection photodiode deflection in volts
+ param_guesser function that guesses initial model parameters
+ model parametric model of deflection vs. z
+ sensitivity_from_fit_params given fit params, return the sensitivity
+ plot boolean overriding matplotlib config setting.
+ Returns:
+ photodiode_sensitivity photodiode volts per piezo meter
+ """
+ _LOG.debug('fit bump data with model %s' % model)
+ def residual(p, deflection, z):
+ return model(z, p) - deflection
+ param_guess = param_guesser(z, deflection)
+ p,cov,info,mesg,ier = _leastsq(
+ residual, param_guess, args=(deflection, z), full_output=True,
+ maxfev=int(10e3))
+ _LOG.debug('fitted params: %s' % p)
+ _LOG.debug('covariance matrix: %s' % cov)
+ #_LOG.debug('info: %s' % info)
+ _LOG.debug('message: %s' % mesg)
+ if ier == 1:
+ _LOG.debug('solution converged')
+ else:
+ _LOG.debug('solution did not converge')
+ if plot or _base_config['matplotlib']:
+ yguess = model(z, param_guess)
+ yfit = model(z, p)
+ bump_plot({'z': z, 'deflection': deflection}, yguess=yguess, yfit=yfit)
return sensitivity_from_fit_params(p)
-@splittableKwargsFunction()
-def bump_save(data, log_dir=None) :
- "Save the dictionary data, using data_logger.data_log()"
- if log_dir != None :
- log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
- log_name="bump")
- log.write_dict_of_arrays(data)
-
-def bump_load(datafile) :
- "Load the dictionary data, using data_logger.date_load()"
- dl = data_logger.data_load()
- data = dl.read_dict_of_arrays(datafile)
- return data
-
-@splittableKwargsFunction()
-def bump_plot(data, yguess=None, yfit=None, plotVerbose=False) :
- "Plot the bump (Vphoto vs Vzp) if plotVerbose or PYLAB_VERBOSE == True"
- if plotVerbose or config.PYLAB_VERBOSE :
- common._import_pylab()
- common._pylab.figure(config.BASE_FIGNUM)
- if yfit != None: # two subplot figure
- common._pylab.subplot(211)
- common._pylab.hold(False)
- common._pylab.plot(data["Z piezo output"], data["Deflection input"],
- '.', label='bump')
- common._pylab.hold(True)
- if yguess != None:
- common._pylab.plot(data["Z piezo output"], yguess,
- 'g-', label='guess')
- if yfit != None:
- common._pylab.plot(data["Z piezo output"], yfit,
- 'r-', label='fit')
- common._pylab.hold(False)
- common._pylab.title("bump surface")
- common._pylab.legend(loc='upper left')
- common._pylab.xlabel("Z piezo output voltage (bits)")
- common._pylab.ylabel("Photodiode input voltage (bits)")
- if yfit != None:
- # second subplot for residual
- common._pylab.subplot(212)
- common._pylab.plot(data["Z piezo output"],
- data["Deflection input"] - yfit,
- 'r-', label='residual')
- common._pylab.legend(loc='upper right')
- common._pylab.xlabel("Z piezo output voltage (bits)")
- common._pylab.ylabel("Photodiode input voltage (bits)")
- common._flush_plot()
-
-make_splittable_kwargs_function(bump_analyze,
- (bump_fit, 'zpiezo_output_bits',
- 'deflection_input_bits'),
- (slope_bitspbit2Vpnm, 'slope_bitspbit'))
-
-@splittableKwargsFunction((bump_analyze, 'data'))
-def bump_load_analyze_tweaked(tweak_file, **kwargs):
- "Load the output file of tweak_calib_bump.sh, return an array of slopes"
- bump_analyze_kwargs, = \
- bump_load_analyze_tweaked._splitargs(bump_load_analyze_tweaked, kwargs)
- photoSensitivity = []
- for line in file(tweak_file, 'r') :
- parsed = line.split()
- path = parsed[0].strip()
- if path[0] == '#' : # a comment
- continue
- if config.TEXT_VERBOSE :
- print "Reading data from %s with ranges %s" % (path, parsed[1:])
- # read the data
- full_data = bump_load(path)
- if len(parsed) == 1 :
- data = full_data # use whole bump
- else :
- # use the listed sections
- zp = []
- df = []
- for rng in parsed[1:] :
- p = rng.split(':')
- starti = int(p[0])
- stopi = int(p[1])
- zp.extend(full_data['Z piezo output'][starti:stopi])
- df.extend(full_data['Deflection input'][starti:stopi])
- data = {'Z piezo output': numpy.array(zp),
- 'Deflection input': numpy.array(df)}
- pSi = bump_analyze(data, **bump_analyze_kwargs)
- photoSensitivity.append(pSi)
- return numpy.array(photoSensitivity, dtype=numpy.float)
-
-# commandline interface functions
-import scipy.io, sys
-
-def read_data(ifile):
- "ifile can be a filename string or open (seekable) file object"
- if ifile == None : ifile = sys.stdin
- unlabeled_data=scipy.io.read_array(ifile)
- data = {}
- data['Z piezo output'] = unlabeled_data[:,0]
- data['Deflection input'] = unlabeled_data[:,1]
- return data
-
-def remove_further_than(data, zp_crit) :
- ndata = {}
- ndata['Z piezo output'] = []
- ndata['Deflection input'] = []
- for zp,df in zip(data['Z piezo output'],data['Deflection input']) :
- if zp > zp_crit :
- ndata['Z piezo output'].append(zp)
- ndata['Deflection input'].append(df)
- return ndata
-
-if __name__ == '__main__' :
- # command line interface
- from optparse import OptionParser
-
- usage_string = ('%prog <input-file>\n'
- '2008, W. Trevor King.\n'
- '\n'
- 'There are two operation modes, one to analyze a single bump file,\n'
- 'and one to analyze tweak files.\n'
- '\n'
- 'Single file mode (the default) :\n'
- 'Scales raw DAC/ADC bit data and fits a bounded quadratic.\n'
- 'Returns photodiode sensitivity Vphotodiode/Zcantilever in V/nm, determined by.\n'
- 'the slope at the kink between the non-contact region and the contact region.\n'
- '<input-file> should be whitespace-delimited, 2 column ASCII\n'
- 'without a header line. e.g: "<zp_DAC>\\t<deflection_ADC>\\n"\n'
- '\n'
- 'Tweak file mode:\n'
- 'Runs the same analysis as in single file mode for each bump in\n'
- 'a tweak file. Each line in the tweak file specifies a single bump.\n'
- 'Blank lines and those beginning with a pound sign (#) are ignored.\n'
- 'The format of a line is a series of whitespace-separated fields--\n'
- 'a base file path followed by optional point index ranges, e.g.:\n'
- '20080919/20080919132500_bump_surface 10:651 1413:2047\n'
- 'which only discards all points outside the index ranges [10,651)\n'
- 'and [1413,2047) (indexing starts at 0).\n'
- )
- parser = OptionParser(usage=usage_string, version='%prog '+common.VERSION)
- parser.add_option('-o', '--output-file', dest='ofilename',
- help='write output to FILE (default stdout)',
- type='string', metavar='FILE')
- parser.add_option('-c', '--comma-out', dest='comma_out', action='store_true',
- help='Output comma-seperated values (default %default)',
- default=False)
- parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
- help='Produce pylab fit checks during execution',
- default=False)
- parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
- help='Run in tweak-file mode',
- default=False)
- parser.add_option('-d', '--datalogger-mode', dest='datalogger_mode', action='store_true',
- help='Run input files with datalogger.read_dict_of_arrays(). This is useful, for example, to test a single line from a tweakfile.',
- default=False)
- parser.add_option('-q', '--disable-quadratic', dest='quadratic', action='store_false',
- help='Disable quadratic term in fitting (i.e. use bounded linear fits).',
- default=True)
- parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
- help='Print lots of debugging information',
- default=False)
-
- options,args = parser.parse_args()
- parser.destroy()
- assert len(args) >= 1, "Need an input file"
-
- ifilename = args[0]
-
- if options.ofilename != None :
- ofile = file(options.ofilename, 'w')
- else :
- ofile = sys.stdout
- config.TEXT_VERBOSE = options.verbose
- config.PYLAB_INTERACTIVE = False
- config.PYLAB_VERBOSE = options.pylab
- config.GNUPLOT_VERBOSE = False
- if options.quadratic == True:
- param_guesser = limited_quadratic_param_guess
- model = limited_quadratic
- sensitivity_from_fit_params = limited_quadratic_sensitivity
+def bump_save(filename, group='/', raw_bump=None, bump_config=None,
+ z_channel_config=None, z_axis_config=None,
+ deflection_channel_config=None, processed_bump=None):
+ with _h5py.File(filename, 'a') as f:
+ cwg = _h5_create_group(f, group)
+ if raw_bump is not None:
+ try:
+ del cwg['raw/z']
+ except KeyError:
+ pass
+ try:
+ del cwg['raw/deflection']
+ except KeyError:
+ pass
+ cwg['raw/z'] = raw_bump['z']
+ cwg['raw/deflection'] = raw_bump['deflection']
+ for config,key in [(bump_config, 'config/bump'),
+ (z_channel_config, 'config/z/channel'),
+ (z_axis_config, 'config/z/axis'),
+ (deflection_channel_config,
+ 'config/deflection/channel')]:
+ if config is None:
+ continue
+ config_cwg = _h5_create_group(cwg, key)
+ config.save(group=config_cwg)
+ if processed_bump is not None:
+ try:
+ del cwg['processed']
+ except KeyError:
+ pass
+ cwg['processed'] = processed_bump
+
+def bump_load(filename, group='/'):
+ assert group.endswith('/')
+ raw_bump = processed_bump = None
+ configs = []
+ with _h5py.File(filename, 'a') as f:
+ try:
+ raw_bump = {
+ 'z': f[group+'raw/z'][...],
+ 'deflection': f[group+'raw/deflection'][...],
+ }
+ except KeyError:
+ pass
+ for Config,key in [(_HDF5_BumpConfig, 'config/bump'),
+ (_HDF5_ChannelConfig, 'config/z/channel'),
+ (_HDF5_AxisConfig, 'config/z/axis'),
+ (_HDF5_ChannelConfig, 'config/deflection/channel')]:
+ config = Config(filename=filename, group=group+key)
+ configs.append(config)
+ try:
+ processed_bump = f[group+'processed'][...]
+ except KeyError:
+ pass
+ ret = [raw_bump]
+ ret.extend(configs)
+ ret.append(processed_bump)
+ for config in configs:
+ config.load()
+ return tuple(ret)
+
+def bump_plot(data, yguess=None, yfit=None):
+ "Plot the bump (Vphoto vs Vzp)"
+ if not _matplotlib:
+ raise _matplotlib_import_error
+ figure = _matplotlib_pyplot.figure()
+ if yfit is None:
+ axes = figure.add_subplot(1, 1, 1)
else:
- param_guesser = limited_linear_param_guess
- model = limited_linear
- sensitivity_from_fit_params = limited_linear_sensitivity
-
- if options.tweakmode == False :
- if options.datalogger_mode:
- data = bump_load(ifilename)
- else:
- data = read_data(ifilename)
- photoSensitivity = bump_analyze(data,
- param_guesser=param_guesser,
- model=model,
- sensitivity_from_fit_params=sensitivity_from_fit_params)
-
- print >> ofile, photoSensitivity
- else : # tweak file mode
- slopes = bump_load_analyze_tweaked(ifilename,
- param_guesser=param_guesser,
- model=model,
- sensitivity_from_fit_params=sensitivity_from_fit_params)
- if options.comma_out :
- sep = ','
- else :
- sep = '\n'
- common.write_array(ofile, slopes, sep)
-
- if options.ofilename != None :
- ofile.close()
+ axes = figure.add_subplot(2, 1, 1)
+ residual_axes = figure.add_subplot(2, 1, 2)
+ timestamp = _time.strftime('%H%M%S')
+
+ axes.hold(True)
+ axes.plot(data['z'], data['deflection'], '.', label='bump')
+ if yguess != None:
+ axes.plot(data['z'], yguess, 'g-', label='guess')
+ if yfit != None:
+ axes.plot(data['z'], yfit, 'r-', label='fit')
+
+ axes.set_title('bump surface %s' % timestamp)
+ #axes.legend(loc='upper left')
+ axes.set_xlabel('Z piezo (meters)')
+ axes.set_ylabel('Photodiode (Volts)')
+ if yfit is not None:
+ # second subplot for residual
+ residual_axes.plot(data['z'], data['deflection'] - yfit,
+ 'r-', label='residual')
+ #residual_axes.legend(loc='upper right')
+ residual_axes.set_xlabel('Z piezo (meters)')
+ residual_axes.set_ylabel('Photodiode (Volts)')
+ figure.show()
# License along with calibcant. If not, see
# <http://www.gnu.org/licenses/>.
-"""
-Aquire and analyze cantilever calibration data.
-
-W. Trevor King Dec. 2007-Jan. 2008
+"""Acquire and analyze cantilever calibration data.
-GPL BOILERPLATE
+The relevent physical quantities are:
+* Vzp_out Output z-piezo voltage (what we generate)
+* Vzp Applied z-piezo voltage (after external ZPGAIN)
+* Zp The z-piezo position
+* Zcant The cantilever vertical deflection
+* Vphoto The photodiode vertical deflection voltage (what we measure)
+* Fcant The force on the cantilever
+* T The temperature of the cantilever and surrounding solution
+* (another thing we measure or guess)
+* k_b Boltzmann's constant
-The relevent physical quantities are :
- Vzp_out Output z-piezo voltage (what we generate)
- Vzp Applied z-piezo voltage (after external ZPGAIN)
- Zp The z-piezo position
- Zcant The cantilever vertical deflection
- Vphoto The photodiode vertical deflection voltage (what we measure)
- Fcant The force on the cantilever
- T The temperature of the cantilever and surrounding solution
- (another thing we measure or guess)
- k_b Boltzmann's constant
+Which are related by the parameters:
-Which are related by the parameters :
- zpGain Vzp_out / Vzp
- zpSensitivity Zp / Vzp
- photoSensitivity Vphoto / Zcant
- k_cant Fcant / Zcant
+* zpGain Vzp_out / Vzp
+* zpSensitivity Zp / Vzp
+* photoSensitivity Vphoto / Zcant
+* k_cant Fcant / Zcant
Cantilever calibration assumes a pre-calibrated z-piezo
(zpSensitivity) and a amplifier (zpGain). In our lab, the z-piezo is
well defined sizes, and the gain is set with a knob on the Nanoscope.
photoSensitivity is measured by bumping the cantilever against the
-surface, where Zp = Zcant (see bump_aquire() and the bump_analyze
+surface, where Zp = Zcant (see bump_acquire() and the bump_analyze
submodule).
k_cant is measured by watching the cantilever vibrate in free solution
-(see the vib_aquire() and the vib_analyze submodule). The average
+(see the vib_acquire() and the vib_analyze submodule). The average
energy of the cantilever in the vertical direction is given by the
equipartition theorem.
- 1/2 k_b T = 1/2 k_cant <Zcant**2>
- so k_cant = k_b T / Zcant**2
- but Zcant = Vphoto / photoSensitivity
- so k_cant = k_b T * photoSensitivty**2 / <Vphoto**2>
+
+.. math:: \frac{1}{2} k_b T = \frac{1}{2} k_cant <Zcant**2>
+
+so
+
+.. math:: k_cant = \frac{k_b T}{Zcant**2}
+
+but
+
+.. math:: Zcant = \frac{Vphoto}{photoSensitivity}
+
+so
+
+.. math:: k_cant = \frac{k_b T * photoSensitivty^2}{<Vphoto**2>}
We measured photoSensitivity with the surface bumps. We can either
measure T using an external function (see temperature.py), or just
-estimate it (see T_aquire() and the T_analyze submodule). Guessing
+estimate it (see T_acquire() and the T_analyze submodule). Guessing
room temp ~22 deg C is actually fairly reasonable. Assuming the
actual fluid temperature is within +/- 5 deg, the error in the spring
constant k_cant is within 5/273.15 ~= 2%. A time series of Vphoto
We do all these measurements a few times to estimate statistical
errors.
-The functions are layed out in the families:
- bump_*(), vib_*(), T_*(), and calib_*()
-where calib_{save|load|analyze}() deal with derived data, not
-real-world data.
-
-For each family, * can be any of :
- aquire get real-world data
- save store real-world data to disk
- load get real-world data from disk
- analyze interperate the real-world data.
- plot show a nice graphic to convince people we're working :p
- load_analyze_tweaked
- read a file with a list of paths to previously saved
- real world data load each file using *_load(), analyze
- using *_analyze(), and optionally plot using *_plot().
- Intended for re-processing old data.
-A family name without any _* extension (e.g. bump()), runs *_aquire(),
- *_save(), *_analyze(), *_plot().
-
-We also define the two positioning functions:
- move_just_onto_surface() and move_far_from_surface()
-which make automating the calibration procedure more straightforward.
-"""
+The functions are layed out in the families::
-import numpy
-import time
+ bump_*(), vib_*(), T_*(), and calib_*()
-import FFT_tools
-import piezo.z_piezo_utils as z_piezo_utils
-from splittable_kwargs import splittableKwargsFunction, \
- make_splittable_kwargs_function
+For each family, * can be any of:
-from . import common
-from . import config
-from . import bump_analyze
-from . import T_analyze
-from . import vib_analyze
-from . import analyze
+* acquire get real-world data
+* save store real-world data to disk
+* load get real-world data from disk
+* analyze interperate the real-world data.
+* plot show a nice graphic to convince people we're working :p
+A family name without any `_*` extension (e.g. `bump()`), runs
+`*_acquire()`, `*_analyze()`, and `*_save()`. `*_analyze()` will run
+`*_plot()` if `matplotlib` is set in `calibcant.base_config`.
+"""
-# bump family
+from numpy import zeros as _zeros
+from numpy import float as _float
+from time import sleep as _sleep
-@splittableKwargsFunction()
-def bump_aquire(zpiezo, push_depth=200, npoints=1024, push_speed=1000) :
- """
- Ramps closer push_depth and returns to the original position.
- Inputs:
- zpiezo an opened zpiezo.zpiezo instance
- push_depth distance to approach, in nm
- npoints number points during the approach and during the retreat
- push_speed piezo speed during approach and retreat, in nm/s
- Returns the aquired ramp data dictionary, with data in DAC/ADC bits.
- """
- # generate the bump output
- nm_per_step = float(push_depth) / npoints
- freq = push_speed / nm_per_step # freq is sample frequency in Hz
- start_pos = zpiezo.curPos()
- pos_dist = zpiezo.pos_nm2out(push_depth) - zpiezo.pos_nm2out(0)
- close_pos = start_pos + pos_dist
- appr = numpy.linspace(start_pos, close_pos, npoints)
- retr = numpy.linspace(close_pos, start_pos, npoints)
- out = numpy.concatenate((appr, retr))
- # run the bump, and measure deflection
- if config.TEXT_VERBOSE :
- print "Bump %g nm at %g nm/s" % (push_depth, push_speed)
- data = zpiezo.ramp(out, freq)
- return data
-
-@splittableKwargsFunction(bump_aquire,
- (bump_analyze.bump_save, 'data'),
- (bump_analyze.bump_analyze, 'data'))
-def bump(**kwargs):
- """
- Wrapper around bump_aquire(), bump_save(), bump_analyze()
- """
- bump_aquire_kwargs,bump_save_kwargs,bump_analyze_kwargs = \
- bump._splitargs(bump, kwargs)
- data = bump_aquire(**bump_aquire_kwargs)
- bump_analyze.bump_save(data, **bump_save_kwargs)
- photoSensitivity = bump_analyze.bump_analyze(data, **bump_analyze_kwargs)
- return photoSensitivity
-
-# T family.
-# Fairly stubby, since a one shot Temp measurement is a common thing.
-# We just wrap that to provide a consistent interface.
-
-@splittableKwargsFunction()
-def T_aquire(get_T=None) :
- """
- Measure the current temperature of the sample,
- or, if get_T == None, fake it by returning config.DEFAULT_TEMP
- """
- if get_T == None :
- if config.TEXT_VERBOSE :
- print "Fake temperature %g" % config.DEFAULT_TEMP
- return config.DEFAULT_TEMP
- else :
- if config.TEXT_VERBOSE :
- print "Measure temperature"
- return get_T()
-
-@splittableKwargsFunction(T_aquire,
- (T_analyze.T_save, 'T'),
- (T_analyze.T_analyze, 'T'))
-def T(**kwargs):
- """
- Wrapper around T_aquire(), T_save(), T_analyze(), T_plot()
- """
- T_aquire_kwargs,T_save_kwargs,T_analyze_kwargs = \
- T._splitargs(T, kwargs)
- T_raw = T_aquire(**T_aquire_kwargs)
- T_analyze.T_save(T_raw, **T_save_kwargs)
- T_ret = T_analyze.T_analyze(T_raw, **T_analyze_kwargs) # returns array
- return T_ret[0]
+from . import LOG as _LOG
-# vib family
+from .bump import bump as _bump
+from .T import T as _T
+from .vib import vib as _vib
+from .analyze import calib_analyze as _calib_analyze
+from .analyze import calib_save as _calib_save
-@splittableKwargsFunction()
-def vib_aquire(zpiezo, time=1, freq=50e3) :
- """
- Record data for TIME seconds at FREQ Hz from ZPIEZO at it's current position.
- """
- # round up to the nearest power of two, for efficient FFT-ing
- nsamps = FFT_tools.ceil_pow_of_two(time*freq)
- time = nsamps / freq
- # take some data, keeping the position voltage at it's current value
- out = numpy.ones((nsamps,), dtype=numpy.uint16) * zpiezo.curPos()
- if config.TEXT_VERBOSE :
- print "get %g seconds of data" % time
- data = zpiezo.ramp(out, freq)
- data['sample frequency Hz'] = numpy.array([freq])
- return data
-
-@splittableKwargsFunction(vib_aquire,
- (vib_analyze.vib_save, 'data'),
- (vib_analyze.vib_analyze, 'deflection_bits', 'freq'))
-def vib(**kwargs) :
- """
- Wrapper around vib_aquire(), vib_save(), vib_analyze()
- """
- vib_aquire_kwargs,vib_save_kwargs,vib_analyze_kwargs = \
- vib._splitargs(vib, kwargs)
- data = vib_aquire(**vib_aquire_kwargs)
- vib_analyze.vib_save(data, **vib_save_kwargs)
- freq = data['sample frequency Hz']
- deflection_bits = data['Deflection input']
- Vphoto_var = vib_analyze.vib_analyze(deflection_bits=deflection_bits,
- freq=freq, **vib_analyze_kwargs)
- return Vphoto_var
-
-# A few positioning functions, so we can run bump_aquire() and vib_aquire()
-# with proper spacing relative to the surface.
-
-@splittableKwargsFunction()
-def move_just_onto_surface(stepper, zpiezo, Depth_nm=-50, setpoint=2) :
- """
- Uses z_piezo_utils.getSurfPos() to pinpoint the position of the
- surface. Adjusts the stepper position as required to get within
- stepper_tol nm of the surface. Then set Vzp to place the
- cantilever Depth_nm onto the surface. Negative Depth_nm values
- will place the cantilever that many nm _off_ the surface.
-
- If getSurfPos() fails to find the surface, backs off (for safety)
- and steps in (without moving the zpiezo) until Vphoto > setpoint.
- """
- stepper_tol = 250 # nm, generous estimate of the fullstep stepsize
-
- if config.TEXT_VERBOSE :
- print "moving just onto surface"
- # Zero the piezo
- if config.TEXT_VERBOSE :
- print "zero the z piezo output"
- zpiezo.jumpToPos(zpiezo.pos_nm2out(0))
- # See if we're near the surface already
- if config.TEXT_VERBOSE :
- print "See if we're starting near the surface"
- try :
- dist = zpiezo.pos_out2nm( \
- z_piezo_utils.getSurfPos(zpiezo, zpiezo.def_V2in(setpoint))
- )
- except (z_piezo_utils.tooClose, z_piezo_utils.poorFit), string :
- if config.TEXT_VERBOSE :
- print "distance failed with: ", string
- print "Back off 200 half steps"
- # Back away 200 steps
- stepper.step_rel(-400)
- stepper.step_rel(200)
- sp = zpiezo.def_V2in(setpoint) # sp = setpoint in bits
- zpiezo.updateInputs()
- cd = zpiezo.curDef() # cd = current deflection in bits
- if config.TEXT_VERBOSE :
- print "Single stepping approach"
- while cd < sp :
- if config.TEXT_VERBOSE :
- print "deflection %g < setpoint %g. step closer" % (cd, sp)
- stepper.step_rel(2) # Full step in
- zpiezo.updateInputs()
- cd = zpiezo.curDef()
- # Back off two steps (protecting against backlash)
- if config.TEXT_VERBOSE :
- print "Step back 4 half steps to get off the setpoint"
- stepper.step_rel(-200)
- stepper.step_rel(196)
- # get the distance to the surface
- zpiezo.updateInputs()
- if config.TEXT_VERBOSE :
- print "get surf pos, with setpoint %g (%d)" % (setpoint, zpiezo.def_V2in(setpoint))
- for i in range(20) : # HACK, keep stepping back until we get a distance
- try :
- dist = zpiezo.pos_out2nm( \
- z_piezo_utils.getSurfPos(zpiezo,zpiezo.def_V2in(setpoint)))
- except (z_piezo_utils.tooClose, z_piezo_utils.poorFit), string :
- stepper.step_rel(-200)
- stepper.step_rel(198)
- continue
- break
- if i >= 19 :
- print "tried %d times, still too close! bailing" % i
- print "probably an invalid setpoint."
- raise Exception, "weirdness"
- if config.TEXT_VERBOSE :
- print 'distance to surface ', dist, ' nm'
- # fine tune the stepper position
- while dist < -stepper_tol : # step back if we need to
- stepper.step_rel(-200)
- stepper.step_rel(198)
- dist = zpiezo.pos_out2nm( \
- z_piezo_utils.getSurfPos(zpiezo, zpiezo.def_V2in(setpoint)))
- if config.TEXT_VERBOSE :
- print 'distance to surface ', dist, ' nm, step back'
- while dist > stepper_tol : # and step forward if we need to
- stepper.step_rel(2)
- dist = zpiezo.pos_out2nm( \
- z_piezo_utils.getSurfPos(zpiezo, zpiezo.def_V2in(setpoint)))
- if config.TEXT_VERBOSE :
- print 'distance to surface ', dist, ' nm, step closer'
- # now adjust the zpiezo to place us just onto the surface
- target = dist + Depth_nm
- zpiezo.jumpToPos(zpiezo.pos_nm2out(target))
- # and we're there :)
- if config.TEXT_VERBOSE :
- print "We're %g nm into the surface" % Depth_nm
-
-@splittableKwargsFunction()
-def move_far_from_surface(stepper, um_back=50) :
- """
- Step back a specified number of microns.
- (uses very rough estimate of step distance at the moment)
- """
- step_nm = 100
- steps = int(um_back*1000/step_nm)
- print "step back %d steps" % steps
- stepper.step_rel(-steps)
-
-
-# and finally, the calib family
-
-@splittableKwargsFunction((move_just_onto_surface, 'stepper', 'zpiezo'),
- (bump, 'zpiezo', 'log_dir', 'Vphoto_in2V'),
- (move_far_from_surface, 'stepper'),
- (T, 'log_dir'),
- (vib, 'zpiezo', 'log_dir', 'Vphoto_in2V'),
- (analyze.calib_save, 'bumps','Ts','vibs','log_dir'))
-def calib_aquire(stepper, zpiezo, num_bumps=10, num_Ts=10, num_vibs=20,
- log_dir=config.LOG_DIR, Vphoto_in2V=config.Vphoto_in2V,
- **kwargs):
+
+def move_far_from_surface(stepper, distance):
+ """Step back approximately `distance` meters.
"""
- Aquire data for calibrating a cantilever in one function.
- return (bump, T, vib), each of which is an array.
- Inputs :
- stepper a stepper.stepper_obj for coarse Z positioning
- zpiezo a z_piezo.z_piezo for fine positioning and deflection readin
- num_bumps number of 'bumps' (see Outputs)
- num_temps number of 'Ts' (see Outputs)
- num_vibs number of 'vib's (see Outputs)
- log_dir directory to log data to. Default 'None' disables logging.
- Vphoto_in2V function to convert photodiode input bits to Volts
-
- + other kwargs. Run calib_aquire._kwargs(calib_aquire) to see
- all options. Run calib_aquire._childSplittables to see a list
- of kwarg functions that this function calls.
-
- Outputs (all are arrays of recorded data) :
- bumps measured (V_photodiode / nm_tip) proportionality constant
- Ts measured temperature (K)
- vibs measured V_photodiode variance in free solution
+ steps = int(distance/stepper.step_size)
+ _LOG.info('step back %d steps (~%g m)' % (steps, distance))
+ stepper.step_relative(-steps)
+
+def calib_acquire(afm, calibration_config, bump_config, temperature_config,
+ vibration_config, filename=None, group='/'):
+ """Acquire data for calibrating a cantilever in one function.
+
+ Inputs:
+ afm a pyafm.AFM instance
+ calibration_config a .config._CalibrationConfig instance
+ bump_config a .config._BumpConfig instance
+ temperature_config a .config._TConfig instance
+ vibration_config a .config._VibrationConfig instance
+
+ Outputs (all are arrays of recorded data):
+ bumps measured (V_photodiode / nm_tip) proportionality constant
+ Ts measured temperature (K)
+ vibs measured V_photodiode variance (Volts**2) in free solution
+
+ The temperatures are collected after moving far from the surface
+ but before and vibrations are measured to give everything time to
+ settle after the big move.
"""
- move_just_onto_surface_kwargs,bump_kwargs,move_far_from_surface_kwargs, \
- T_kwargs,vib_kwargs,calib_save_kwargs = \
- calib_aquire._splitargs(calib_aquire, kwargs)
- # get bumps
- bumps = numpy.zeros((num_bumps,), dtype=numpy.float)
- for i in range(num_bumps) :
- move_just_onto_surface(stepper, zpiezo, **move_just_onto_surface_kwargs)
- bumps[i] = bump(zpiezo=zpiezo, log_dir=log_dir,
- Vphoto_in2V=Vphoto_in2V, **bump_kwargs)
- if config.TEXT_VERBOSE :
- print bumps
-
- move_far_from_surface(stepper, **move_far_from_surface_kwargs)
-
- # get Ts
- Ts = numpy.zeros((num_Ts,), dtype=numpy.float)
- for i in range(num_Ts) :
- Ts[i] = T(**T_kwargs)
- time.sleep(1) # wait a bit to get an independent temperature measure
- print Ts
+ assert group.endswith('/'), group
+
+ bumps = _zeros((calibration_config['num-bumps'],), dtype=_float)
+ for i in range(calibration_config['num-bumps']):
+ _LOG.info('acquire bump %d of %d' % (i, calibration_config['num-bumps']))
+ bumps[i] = _bump(afm=afm, bump_config=bump_config,
+ filename=filename, group='%sbump/%d/' % (group, i))
+ _LOG.debug('bumps: %s' % bumps)
+
+ move_far_from_surface(
+ afm.stepper, distance=calibration_config['vibration-spacing'])
+
+ Ts = _zeros((calibration_config['num-temperatures'],), dtype=_float)
+ for i in range(calibration_config['num-temperatures']):
+ _LOG.info('acquire T %d of %d'
+ % (i, calibration_config['num-temperatures']))
+ Ts[i] = _T(
+ get_T=afm.get_temperature, temperature_config=temperature_config,
+ filename=filename, group='%stemperature/%d/' % (group, i))
+ _sleep(calibration_config['temperature-sleep'])
+ _LOG.debug('temperatures: %s' % Ts)
# get vibs
- vibs = numpy.zeros((num_vibs,), dtype=numpy.float)
- for i in range(num_vibs) :
- vibs[i] = vib(zpiezo=zpiezo, log_dir=log_dir, Vphoto_in2V=Vphoto_in2V,
- **vib_kwargs)
- print vibs
-
- analyze.calib_save(bumps, Ts, vibs, log_dir, **calib_save_kwargs)
-
+ vibs = _zeros((calibration_config['num-vibrations'],), dtype=_float)
+ for i in range(calibration_config['num-vibrations']):
+ vibs[i] = _vib(
+ piezo=afm.piezo, vibration_config=vibration_config,
+ filename=filename, group='%svibration/%d/' % (group, i))
+ _LOG.debug('vibrations: %s' % vibs)
+
return (bumps, Ts, vibs)
+def calib(afm, calibration_config, bump_config, temperature_config,
+ vibration_config, filename=None, group='/'):
+ """Calibrate a cantilever in one function.
-@splittableKwargsFunction( \
- (calib_aquire, 'log_dir'),
- (analyze.calib_analyze, 'bumps','Ts','vibs'))
-def calib(log_dir=config.LOG_DIR, **kwargs) :
- """
- Calibrate a cantilever in one function.
- The I-don't-care-about-the-details black box version :p.
- return (k, k_s)
Inputs:
- (see calib_aquire())
- Outputs :
- k cantilever spring constant (in N/m, or equivalently nN/nm)
- k_s standard deviation in our estimate of k
- Notes :
- See get_calibration_data() for the data aquisition code
- See analyze_calibration_data() for the analysis code
+ (see `calib_acquire()`)
+
+ Outputs:
+ k cantilever spring constant (in N/m, or equivalently nN/nm)
+ k_s standard deviation in our estimate of k
+
+ >>> import os
+ >>> from pprint import pprint
+ >>> import tempfile
+ >>> from pycomedi.device import Device
+ >>> from pycomedi.subdevice import StreamingSubdevice
+ >>> from pycomedi.channel import AnalogChannel, DigitalChannel
+ >>> from pycomedi.constant import AREF, IO_DIRECTION, SUBDEVICE_TYPE, UNIT
+ >>> from pypiezo.afm import AFMPiezo
+ >>> from pypiezo.base import PiezoAxis, InputChannel
+ >>> from pypiezo.config import (HDF5_ChannelConfig, HDF5_AxisConfig,
+ ... pprint_HDF5)
+ >>> from stepper import Stepper
+ >>> from pyafm import AFM
+ >>> from .config import (HDF5_CalibrationConfig, HDF5_BumpConfig,
+ ... HDF5_TemperatureConfig, HDF5_VibrationConfig)
+ >>> from .analyze import calib_load_all
+
+ >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
+ >>> os.close(fd)
+
+ >>> d = Device('/dev/comedi0')
+ >>> d.open()
+
+ Setup an `AFMPiezo` instance.
+
+ >>> s_in = d.find_subdevice_by_type(SUBDEVICE_TYPE.ai,
+ ... factory=StreamingSubdevice)
+ >>> s_out = d.find_subdevice_by_type(SUBDEVICE_TYPE.ao,
+ ... factory=StreamingSubdevice)
+
+ >>> axis_channel = s_out.channel(
+ ... 0, factory=AnalogChannel, aref=AREF.ground)
+ >>> input_channel = s_in.channel(0, factory=AnalogChannel, aref=AREF.diff)
+ >>> for chan in [axis_channel, input_channel]:
+ ... chan.range = chan.find_range(unit=UNIT.volt, min=-10, max=10)
+
+ We set the minimum voltage for the `z` axis to -9 (a volt above
+ the minimum possible voltage) to help with testing
+ `.get_surface_position`. Without this minimum voltage, small
+ calibration errors could lead to a railed -10 V input for the
+ first few surface approaching steps, which could lead to an
+ `EdgeKink` error instead of a `FlatFit` error.
+
+ >>> axis_config = HDF5_AxisConfig(filename, '/bump/config/z/axis')
+ >>> axis_config.update(
+ ... {'gain':20, 'sensitivity':8e-9, 'minimum':-9})
+ >>> axis_channel_config = HDF5_ChannelConfig(
+ ... filename, '/bump/config/z/channel')
+ >>> input_channel_config = HDF5_ChannelConfig(
+ ... filename, '/bump/config/deflection/channel')
+
+ >>> a = PiezoAxis(axis_config=axis_config,
+ ... axis_channel_config=axis_channel_config,
+ ... axis_channel=axis_channel, name='z')
+ >>> a.setup_config()
+
+ >>> c = InputChannel(
+ ... channel_config=input_channel_config, channel=input_channel,
+ ... name='deflection')
+ >>> c.setup_config()
+
+ >>> piezo = AFMPiezo(axes=[a], input_channels=[c])
+
+ Setup a `stepper` instance.
+
+ >>> s_d = d.find_subdevice_by_type(SUBDEVICE_TYPE.dio)
+ >>> d_channels = [s_d.channel(i, factory=DigitalChannel)
+ ... for i in (0, 1, 2, 3)]
+ >>> for chan in d_channels:
+ ... chan.dio_config(IO_DIRECTION.output)
+
+ >>> def write(value):
+ ... s_d.dio_bitfield(bits=value, write_mask=2**4-1)
+
+ >>> stepper = Stepper(write=write)
+
+ Setup an `AFM` instance.
+
+ >>> afm = AFM(piezo, stepper)
+
+ Test calibration:
+
+ >>> calibration_config = HDF5_CalibrationConfig(
+ ... filename=filename, group='/bump/config/calibration/')
+ >>> bump_config = HDF5_BumpConfig(
+ ... filename=filename, group='/bump/config/bump/')
+ >>> temperature_config = HDF5_TemperatureConfig(
+ ... filename=filename, group='/bump/config/temperature/')
+ >>> vibration_config = HDF5_VibrationConfig(
+ ... filename=filename, group='/bump/config/vibration')
+ >>> calib(afm, calibration_config, bump_config, temperature_config,
+ ... vibration_config, filename=filename, group='/')
+ TODO: replace skipped example data with real-world values
+ >>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
+ >>> everything = calib_load_all(filename, '/')
+ >>> pprint(everything)
+
+ Close the Comedi device.
+
+ >>> d.close()
+
+ Cleanup our temporary config file.
+
+ os.remove(filename)
"""
- calib_aquire_kwargs,calib_analyze_kwargs = \
- calib._splitargs(calib, kwargs)
- a, T, vib = calib_aquire(**calib_aquire_kwargs)
- k,k_s,ps2_m, ps2_s,T_m,T_s,one_o_Vp2_m,one_o_Vp2_s = \
- analyze.calib_analyze(a, T, vib, **calib_analyze_kwargs)
- analyze.calib_save_analysis(k, k_s, ps2_m, ps2_s, T_m, T_s,
- one_o_Vp2_m, one_o_Vp2_s, log_dir)
+ bumps, Ts, vibs = calib_acquire(
+ afm, calibration_config, bump_config, temperature_config,
+ vibration_config, filename=filename, group=group)
+ # TODO: convert vib units?
+ k,k_s = _calib_analyze(bumps, Ts, vibs)
+ _calib_save(filename, group=group+'calibration/', bumps=bumps, Ts=Ts,
+ vibs=vibs, calibration_config=calibration_config, k=k, k_s=k_s)
return (k, k_s)
-
-
-
+++ /dev/null
-# calibcant - tools for thermally calibrating AFM cantilevers
-#
-# Copyright (C) 2008-2011 W. Trevor King <wking@drexel.edu>
-#
-# This file is part of calibcant.
-#
-# calibcant is free software: you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation, either
-# version 3 of the License, or (at your option) any later version.
-#
-# calibcant is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with calibcant. If not, see
-# <http://www.gnu.org/licenses/>.
-
-from . import config
-
-
-# handle extra verbose input modules, only imported if we need them
-_flush_plot = None
-_final_flush_plot = None
-_pylab = None
-def _dummy_fn(): pass
-
-def _import_pylab() : # TODO: auto detect no DISPLAY and abort
- """Import pylab plotting functions for when we need to plot. This
- function can be called multiple times, and ensures that the pylab
- setup is only imported once. It defines the functions
- _flush_plot() to be called after incremental changes,
- _final_flush_plot(), to be called at the end of any graphical processing, and
- _pylab, a reference to the pylab module."""
- global _pylab
- global _flush_plot
- global _final_flush_plot
- if _pylab == None :
- import pylab as _pylab
- if _flush_plot == None :
- if config.PYLAB_INTERACTIVE :
- _flush_plot = _pylab.draw
- else :
- _flush_plot = _pylab.show #_dummy_fn
- if _final_flush_plot == None :
- if config.PYLAB_INTERACTIVE :
- _final_flush_plot = _pylab.draw
- else :
- _final_flush_plot = _dummy_fn # _pylab.show
-
-def write_array(ofile, array, seperator):
- """Erite an iterable array seperated by seperator.
- Terminate the output with an endline."""
- strings = [repr(x) for x in array]
- ofile.write(seperator.join(strings))
- ofile.write('\n')
"""Define some variables to configure the package for a particular lab
and workflow."""
-import piezo.z_piezo as z_piezo
-
-
-DEFAULT_TEMP = 22 # assume 22 deg C
-LOG_DATA = True # quietly grab all real-world data and log to LOG_DIR
-LOG_DIR = '${DEFAULT}/calibrate_cantilever'
-GNUFIT_DATA_BASE='./calibrate_cantilever_fitdata'
-TEXT_VERBOSE = True # for debugging
-GNUPLOT_VERBOSE = True # turn on fit check plotting
-PYLAB_VERBOSE = True # turn on plotting
-PYLAB_INTERACTIVE = True # select between draw() and show() for flushing plots
-BASE_FIGNUM = 20 # to avoid writing to already existing figures
-
-
-# HACK
-# make sure you make a system note (add_system_note) if you change
-# these in case you don't have access to a z_piezo for conversion
-# functions
-
-# zpGain zpiezo applied voltage per output Volt
-zpGain = z_piezo.DEFAULT_GAIN
-# zpSensitivity nm zpiezo response per applied Volt
-zpSensitivity = z_piezo.DEFAULT_SENSITIVITY
-# Vzp_out2V function that converts output DAC bits to Volts
-Vzp_out2V = z_piezo.DEFAULT_VZP_OUT_2_VOLTS
-# Vphoto_in2V function that converts input ADC bits to Volts
-Vphoto_in2V = z_piezo.DEFAULT_VPHOTO_IN_2_VOLTS
-# zeroVphoto_bits ADC bit output for a 0V input
-zeroVphoto_bits = z_piezo.DEFAULT_ZERO_PHOTODIODE_BITS
+import logging as _logging
+import os.path as _os_path
+import sys as _sys
+
+from FFT_tools import window_hann as _window_hann
+from pypiezo.config import (
+ _ChoiceSetting, _BooleanSetting, _IntegerSetting, _FloatSetting,
+ _FloatListSetting, _Config, _BackedConfig, _HDF5Config, _YAMLConfig)
+
+from . import LOG as _LOG
+
+
+class _BaseConfig (_Config):
+ "Configure `calibcant` module operation"
+ settings = [
+ _ChoiceSetting(
+ name='log-level',
+ help='Module logging level.',
+ default=_logging.WARN,
+ choices=[
+ ('critical', _logging.CRITICAL),
+ ('error', _logging.ERROR),
+ ('warn', _logging.WARN),
+ ('info', _logging.INFO),
+ ('debug', _logging.DEBUG),
+ ]),
+ _BooleanSetting(
+ name='syslog',
+ help='Log to syslog (otherwise log to stderr).',
+ default=False),
+ _BooleanSetting(
+ name='matplotlib',
+ help='Plot piezo motion using `matplotlib`.',
+ default=False),
+ _FloatSetting(
+ name='temperature',
+ help=('Default temperature for thermal calibration in degrees '
+ 'Celsius.'),
+ default=22),
+ ]
+
+class _TemperatureUnit (object):
+ pass
+class Celsius (_TemperatureUnit):
+ pass
+class Kelvin (_TemperatureUnit):
+ pass
+
+class _TemperatureConfig (_Config):
+ "Configure `calibcant` temperature operation"
+ settings = [
+ _ChoiceSetting(
+ name='units',
+ help='Units of raw temperature measurements.',
+ default=Celsius,
+ choices=[
+ ('Celsius', Celsius),
+ ('Kelvin', Kelvin),
+ ]),
+ _BooleanSetting(
+ name='default',
+ help=('The temperature values are defaults (vs. real '
+ 'measurements).'),
+ default=True),
+ ]
+
+class _BumpModel (object):
+ pass
+class Linear (_BumpModel):
+ pass
+class Quadratic (_BumpModel):
+ pass
+
+class _BumpConfig (_Config):
+ "Configure `calibcant` bump operation"
+ settings = [
+ _FloatSetting(
+ name='initial-position',
+ help=('Position relative to surface for start of bump in meters. '
+ 'Should be less than zero to ensure non-contact region '
+ 'before you hit the surface.'),
+ default=-50e-9),
+ _FloatSetting(
+ name='setpoint',
+ help=('Maximum deflection in volts in case of stepper positioning '
+ 'to achieve the initial position.'),
+ default=2.0),
+ _IntegerSetting(
+ name='far-steps',
+ help=('Number of stepper steps to move "far" away from the '
+ 'surface. For possible stepper adjustments while initially '
+ 'locating the surface.'),
+ default=200),
+ _FloatSetting(
+ name='push-depth',
+ help='Distance to approach in meters.',
+ default=200e-9),
+ _FloatSetting(
+ name='push-speed',
+ help='Approach/retract speed in meters/second.',
+ default=1e-6),
+ _FloatSetting(
+ name='samples',
+ help='Number of samples during approach and during retreat.',
+ default=1024),
+ _ChoiceSetting(
+ name='model',
+ help='Bump deflection model.',
+ default=Quadratic,
+ choices=[
+ ('linear', Linear),
+ ('quadratic', Quadratic),
+ ]),
+ ]
+
+class _VibrationModel (object):
+ pass
+class Variance (_VibrationModel):
+ pass
+class BreitWigner (_VibrationModel):
+ pass
+class OffsetBreitWigner (_VibrationModel):
+ pass
+
+class _VibrationConfig (_Config):
+ "Configure `calibcant` vibration operation"
+ settings = [
+ _FloatSetting(
+ name='frequency',
+ help='Sampling frequency in Hz.',
+ default=50e3),
+ _FloatSetting(
+ name='sample-time',
+ help=('Aquisition time in seconds. This is rounded up as required '
+ 'so the number of samples will be an integer power of two.'),
+ default=1),
+ _ChoiceSetting(
+ name='model',
+ help='Vibration model.',
+ default=BreitWigner,
+ choices=[
+ ('variance', Variance),
+ ('Breit-Wigner', BreitWigner),
+ ('offset Breit-Wigner', OffsetBreitWigner),
+ ]),
+ _IntegerSetting(
+ name='chunk-size',
+ help='FFT chunk size (for PSD fits).',
+ default=2048),
+ _BooleanSetting(
+ name='overlap',
+ help='Overlap FFT chunks (for PSD fits).'),
+ _ChoiceSetting(
+ name='window',
+ help='FFT chunk window (for PSD fits).',
+ default=_window_hann,
+ choices=[
+ ('Hann', _window_hann),
+ ]),
+ _FloatSetting(
+ name='minimum-fit-frequency',
+ help='Lower bound of Lorentzian fitting region.',
+ default=500.),
+ _FloatSetting(
+ name='maximum-fit-frequency',
+ help='Upper bound of Lorentzian fitting region.',
+ default=25e3),
+ ]
+
+
+class _CalibrationConfig (_Config):
+ "Configure a full `calibcant` calibration run"
+ settings = [
+ _IntegerSetting(
+ name='num-bumps',
+ help='Number of surface bumps.',
+ default=10),
+ _IntegerSetting(
+ name='num-temperatures',
+ help='Number of temperature measurements.',
+ default=10),
+ _IntegerSetting(
+ name='num-vibrations',
+ help='Number of thermal vibration measurements.',
+ default=20),
+ _FloatSetting(
+ name='temperature-sleep',
+ help=('Time between temperature measurements (in seconds) to get '
+ 'independent measurements when reading from slow sensors.'),
+ default=1),
+ _FloatSetting(
+ name='vibration-spacing',
+ help=('Approximate distance from the surface in meters for the '
+ 'vibration measurements. This should be large enough that '
+ 'surface effects are negligable.'),
+ default=50e-6),
+ ]
+
+
+# Define HDF5- and YAML-backed subclasses of the basic _Config types.
+for name,obj in locals().items():
+ if (obj != _Config and
+ type(obj) == type and
+ issubclass(obj, _Config) and
+ not issubclass(obj, _BackedConfig)):
+ for prefix,base in [('HDF5', _HDF5Config), ('YAML', _YAMLConfig)]:
+ _name = '%s%s' % (prefix, name)
+ _bases = (base, obj)
+ _dict = {}
+ _class = type(_name, _bases, _dict)
+ setattr(_sys.modules[__name__], _name, _class)
+
+del name, obj, prefix, base, _name, _bases, _dict, _class
+
+
+def find_base_config():
+ "Return the best `_BaseConfig` match after scanning the filesystem"
+ _LOG.info('looking for base_config file')
+ user_basepath = _os_path.join(_os_path.expanduser('~'), '.calibcantrc')
+ system_basepath = _os_path.join('/etc', 'calibcant', 'config')
+ distributed_basepath = _os_path.join(
+ '/usr', 'share', 'calibcant', 'config')
+ for basepath in [user_basepath, system_basepath, distributed_basepath]:
+ for (extension, config) in [('.h5', HDF5_BaseConfig),
+ ('.yaml', YAML_BaseConfig)]:
+ filename = basepath + extension
+ if _os_path.exists(filename):
+ _LOG.info('base_config file found at %s' % filename)
+ base_config = config(filename)
+ base_config.load()
+ return base_config
+ else:
+ _LOG.debug('no base_config file at %s' % filename)
+ _LOG.info('new base_config file at %s' % filename)
+ basepath = user_basepath
+ filename = basepath + extension
+ return config(filename)
+++ /dev/null
-"""
-windowed_filter(array, winfunc, s)
-with winfuncs :
- windowed_mean_point
- windowed_median_point
-
-Command line bindings for filtering 2 column ACSII data.
-"""
-
-import numpy
-
-def windowed_mean_point(array, i, width=50) :
- """
- Filter data with a windowed average (mean).
- The ith point is replaced by the average of the points with indices in
- the range i +/- width (inclusive).
-
- This function returns the ith point in the filtered array,
- given the raw array and window half-width.
-
- The weights on the average are uniform at the moment.
- """
- x=0
- n=0
- for j in range(i-width, i+width) :
- if j >= 0 and j < len(array) :
- x+=array[j]
- n+=1
- x /= float(n)
- return x
-
-def windowed_median_point(array, i, width=9) :
- """
- Filter data with a windowed median.
- The ith point is replaced by the median of the points with indices in
- the range i +/- width (inclusive).
-
- This function returns the ith point in the filtered array,
- given the raw array and window half-width.
- """
- imin = i-width
- if imin < 0 : imin = 0
- imax = i+width
- if imax >= len(array) : imax = len(array-1)
- slice = array[imin:imax+1].copy()
- slice.sort()
- imid = numpy.floor((imax-imin)/2)
- return slice[imid]
-
-def windowed_filter(array, winfunc, width=None) :
- """
- Filter data with a windowing function winfunc.
- The ith point is replaced by the winfunc(array, i, s).
-
- See the windowed_* functions for possible winfunc options.
- """
- out = array.copy()
- if width == None :
- win = lambda i : winfunc(array, i) # user winfunc's default s
- else :
- win = lambda i : winfunc(array, i, width)
- for i in range(len(out)) :
- out[i] = win(i)
- return out
-
-
-# commandline interface functions
-import scipy.io, sys
-
-def read_data(ifile):
- """
- ifile can be a filename string or open (seekable) file object.
- returns (column 1 array, column 2 array)
- """
- if ifile == None : ifile = sys.stdin
- data=scipy.io.read_array(ifile)
- return (data[:,0], data[:,1])
-
-def write_data(ofile, x, y):
- """
- ofile can be a filename string or open (seekable) file object.
- """
- data = numpy.zeros((len(x),2))
- data[:,0] = x
- data[:,1] = y
- if ofile == None : ofile = sys.stdout
- scipy.io.write_array(ofile, data, separator='\t')
-
-if __name__ == '__main__' :
- # command line interface
- from optparse import OptionParser
-
- usage_string = ('%prog <input-file>\n'
- '2008, W. Trevor King.\n'
- '\n'
- 'Apply various filters to data y(x).'
- '<input-file> should be whitespace-delimited, 2 column ASCII\n'
- 'without a header line.\n'
- 'e.g: "<x>\\t<y>\\n"')
- parser = OptionParser(usage=usage_string, version='%prog 0.1')
- parser.add_option('-o', '--output-file', dest='ofilename',
- help='write output to FILE (default stdout)',
- type='string', metavar='FILE')
- parser.add_option('-w', '--width', dest='width', default=None,
- help='window width (i +/- width, inclusive)',
- type='int', metavar='WIDTH')
- parser.add_option('-t', '--type', dest='type', default='mean',
- help='filter type (default %default)',
- type='string', metavar='TYPE')
- parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
- help='Print lots of debugging information',
- default=False)
-
- options,args = parser.parse_args()
- parser.destroy()
- assert len(args) >= 1, "Need an input file"
-
- ifilename = args[0]
-
- x,y = read_data(ifilename)
-
- if options.type == 'mean' :
- winfunc = windowed_mean_point
- elif options.type == 'median' :
- winfunc = windowed_median_point
- else :
- raise Exception, "unrecognized window type '%s'" % options.type
-
- y = windowed_filter(y, winfunc, width=options.width)
-
- write_data(options.ofilename, x, y)
+++ /dev/null
-if __name__ == "__main__" :
-
- import sys
- import filter, calibcant_vib_analyze
- from optparse import OptionParser
-
- usage_string = ('%prog <input-file>\n'
- '2008, W. Trevor King.\n'
- '\n'
- 'Apply various filters to data y(x).'
- '<input-file> should be whitespace-delimited, 2 column ASCII\n'
- 'without a header line.\n'
- 'e.g: "<x>\\t<y>\\n"')
- parser = OptionParser(usage=usage_string, version='%prog 0.1')
- parser.add_option('-o', '--output-file', dest='ofilename',
- help='write output to FILE (default stdout)',
- type='string', metavar='FILE')
- parser.add_option('-w', '--width', dest='width', default=None,
- help='window width (i +/- width, inclusive)',
- type='int', metavar='WIDTH')
- parser.add_option('-t', '--type', dest='type', default='mean',
- help='filter type (default %default)',
- type='string', metavar='TYPE')
- parser.add_option('-m', '--min-freq', dest='min_freq', default=500,
- help='minimum frequency in Hz (default %default)',
- type='float', metavar='FREQ')
- parser.add_option('-M', '--max-freq', dest='max_freq', default=7000,
- help='maximum frequency in Hz (default %default)',
- type='float', metavar='FREQ')
- parser.add_option('-g', '--gnuplot', dest='gnuplot', action='store_true',
- help='Print gnuplot fit check script to stderr',
- default=False)
- parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
- help='Print lots of debugging information',
- default=False)
-
- options,args = parser.parse_args()
- parser.destroy()
- assert len(args) >= 1, "Need an input file"
-
- ifilename = args[0]
-
- # mirror filter behavior:
- freq,psd = filter.read_data(ifilename)
-
- if options.type == 'mean' :
- winfunc = filter.windowed_mean_point
- elif options.type == 'median' :
- winfunc = filter.windowed_median_point
- else :
- raise Exception, "unrecognized window type '%s'" % options.type
-
- psd = filter.windowed_filter(psd, winfunc, width=options.width)
-
- # mirror calibcant_vib_analyze behavior:
- (A,B,C) = calibcant_vib_analyze.fit_psd(freq, psd,
- minFreq=options.min_freq,
- maxFreq=options.max_freq)
- if options.verbose :
- print >> sys.stderr, "Fit f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with"
- print >> sys.stderr, "A = %g, \t B = %g, \t C = %g" % (A, B, C)
- if options.gnuplot :
- print >> sys.stderr, \
- calibcant_vib_analyze.gnuplot_check_fit(ifilename, A, B, C)
-
- area = calibcant_vib_analyze.lorentzian_area(A,B,C)
-
- if options.ofilename != None :
- print >> file(options.ofilename, 'w'), area
- else :
- print area
--- /dev/null
+# calibcant - tools for thermally calibrating AFM cantilevers
+#
+# Copyright (C) 2008-2011 W. Trevor King <wking@drexel.edu>
+#
+# This file is part of calibcant.
+#
+# calibcant is free software: you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation, either
+# version 3 of the License, or (at your option) any later version.
+#
+# calibcant is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with calibcant. If not, see
+# <http://www.gnu.org/licenses/>.
+
+"""Acquire, save, and load cantilever thermal vibration data.
+
+For measuring the cantilever's spring constant.
+
+The relevent physical quantity is:
+ Vphoto The photodiode vertical deflection voltage (what we measure)
+"""
+
+from numpy import zeros as _zeros
+from FFT_tools import ceil_pow_of_two as _ceil_pow_of_two
+from pycomedi.constant import TRIG_SRC as _TRIG_SRC
+from pycomedi.utility import inttrig_insn as _inttrig_insn
+from pycomedi.utility import Reader as _Reader
+
+from . import LOG as _LOG
+from .vib_analyze import vib_analyze as _vib_analyze
+from .vib_analyze import vib_save as _vib_save
+
+
+def vib_acquire(piezo, vibration_config):
+ """Record thermal vibration data for `piezo` at its current position.
+
+ Inputs:
+ piezo a pypiezo.afm.AFMPiezo instance
+ vibration_config a .config._VibrationConfig instance
+ """
+ _LOG.debug('prepare vibration aquisition command')
+
+ # round up to the nearest power of two, for efficient FFT-ing
+ n_samps = _ceil_pow_of_two(
+ vibration_config['sample-time']*vibration_config['frequency'])
+ time = n_samps / vibration_config['frequency']
+ scan_period_ns = int(1e9 / vibration_config['frequency'])
+
+ input_channel = piezo.input_channel_by_name('deflection')
+ channel = input_channel.channel
+
+ channels = [channel]
+
+ dtype = piezo.deflection_dtype()
+ data = _zeros((n_samps, len(channels)), dtype=dtype)
+
+ cmd = channel.subdevice.get_cmd_generic_timed(
+ len(channels), scan_period_ns)
+ cmd.start_src = _TRIG_SRC.int
+ cmd.start_arg = 0
+ cmd.stop_src = _TRIG_SRC.count
+ cmd.stop_arg = n_samps
+ cmd.chanlist = channels
+ channel.subdevice.cmd = cmd
+ for i in range(3):
+ rc = channel.subdevice.command_test()
+ if rc == None: break
+ _LOG.debug('command test %d: %s' % (i,rc))
+
+ _LOG.info('get %g seconds of vibration data at %g Hz'
+ % (vibration_config['sample-time'],
+ vibration_config['frequency']))
+ channel.subdevice.command()
+ reader = _Reader(channel.subdevice, data)
+ reader.start()
+ channel.subdevice.device.do_insn(_inttrig_insn(channel.subdevice))
+ reader.join()
+ data = data.reshape((data.size,))
+ return data
+
+def vib(piezo, vibration_config, filename, group='/'):
+ """Wrapper around vib_acquire(), vib_analyze(), vib_save().
+
+ >>> import os
+ >>> import tempfile
+ >>> from pycomedi.device import Device
+ >>> from pycomedi.subdevice import StreamingSubdevice
+ >>> from pycomedi.channel import AnalogChannel
+ >>> from pycomedi.constant import AREF, SUBDEVICE_TYPE, UNIT
+ >>> from pypiezo.afm import AFMPiezo
+ >>> from pypiezo.base import InputChannel
+ >>> from pypiezo.config import HDF5_ChannelConfig, pprint_HDF5
+ >>> from .config import HDF5_VibrationConfig
+
+ Setup an `AFMPiezo` instance.
+
+ >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
+ >>> os.close(fd)
+
+ >>> d = Device('/dev/comedi0')
+ >>> d.open()
+
+ >>> s_in = d.find_subdevice_by_type(SUBDEVICE_TYPE.ai,
+ ... factory=StreamingSubdevice)
+
+ >>> channel = s_in.channel(0, factory=AnalogChannel, aref=AREF.diff)
+ >>> channel.range = channel.find_range(
+ ... unit=UNIT.volt, min=-10, max=10)
+ >>> channel_config = HDF5_ChannelConfig(
+ ... filename, group='/vibration/config/deflection/channel')
+
+ >>> c = InputChannel(
+ ... channel_config=channel_config, channel=channel,
+ ... name='deflection')
+ >>> c.setup_config()
+
+ >>> piezo = AFMPiezo(axes=[], input_channels=[c])
+
+ Test a vibration:
+
+ >>> vibration_config = HDF5_VibrationConfig(
+ ... filename=filename, group='/vibration/config/vibration')
+ >>> vib(piezo, vibration_config, filename, group='/vibration')
+ TODO: replace skipped example data with real-world values
+ >>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
+ /
+ /vibration
+ /vibration/config
+ /vibration/config/deflection
+ /vibration/config/deflection/channel
+ <HDF5 dataset "channel": shape (), type "|S1">
+ 0
+ <HDF5 dataset "conversion-coefficients": shape (), type "|S24">
+ -10.0, 0.000305180437934
+ <HDF5 dataset "conversion-origin": shape (), type "|S3">
+ 0.0
+ <HDF5 dataset "device": shape (), type "|S12">
+ /dev/comedi0
+ <HDF5 dataset "inverse-conversion-coefficients": shape (), type "|S12">
+ 0.0, 3276.75
+ <HDF5 dataset "inverse-conversion-origin": shape (), type "|S5">
+ -10.0
+ <HDF5 dataset "maxdata": shape (), type "|S5">
+ 65535
+ <HDF5 dataset "range": shape (), type "|S1">
+ 0
+ <HDF5 dataset "subdevice": shape (), type "|S1">
+ 0
+ /vibration/config/vibration
+ <HDF5 dataset "chunk-size": shape (), type "|S4">
+ 2048
+ <HDF5 dataset "frequency": shape (), type "|S7">
+ 50000.0
+ <HDF5 dataset "maximum-fit-frequency": shape (), type "|S7">
+ 25000.0
+ <HDF5 dataset "minimum-fit-frequency": shape (), type "|S5">
+ 500.0
+ <HDF5 dataset "model": shape (), type "|S12">
+ Breit-Wigner
+ <HDF5 dataset "overlap": shape (), type "|S2">
+ no
+ <HDF5 dataset "sample-time": shape (), type "|S1">
+ 1
+ <HDF5 dataset "window": shape (), type "|S4">
+ Hann
+ <HDF5 dataset "processed": shape (), type "<f8">
+ ...
+ /vibration/raw
+ <HDF5 dataset "deflection": shape (65536,), type "<u2">
+ [...]
+
+ Close the Comedi device.
+
+ >>> d.close()
+
+ Cleanup our temporary config file.
+
+ >>> os.remove(filename)
+ """
+ deflection_input_channel = piezo.input_channel_by_name('deflection')
+
+ deflection_channel_config = deflection_input_channel.channel_config
+
+ deflection = vib_acquire(piezo, vibration_config)
+ variance = _vib_analyze(
+ deflection, vibration_config, deflection_channel_config)
+ _vib_save(
+ filename, group, raw_vibration=deflection,
+ vibration_config=vibration_config,
+ deflection_channel_config=deflection_channel_config,
+ processed_vibration=variance)
+ return variance
# License along with calibcant. If not, see
# <http://www.gnu.org/licenses/>.
-"""
-Separate the more general vib_analyze() from the other vib_*()
-functions in calibcant. Also provide a command line interface
-for analyzing data acquired through other workflows.
+"""Thermal vibration analysis.
+
+Separate the more general `vib_analyze()` from the other `vib_*()`
+functions in calibcant.
The relevent physical quantities are :
- Vphoto The photodiode vertical deflection voltage (what we measure)
+ Vphoto The photodiode vertical deflection voltage (what we measure)
+
+>>> import os
+>>> from pprint import pprint
+>>> import random
+>>> import tempfile
+>>> import numpy
+>>> from .config import HDF5_VibrationConfig
+>>> from pypiezo.config import pprint_HDF5, HDF5_ChannelConfig
+>>> from pypiezo.base import convert_volts_to_bits
+
+>>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
+>>> os.close(fd)
+
+>>> vibration_config = HDF5_VibrationConfig(
+... filename=filename, group='/vibration/config/vibration')
+>>> vibration_config['frequency'] = 50e3
+>>> vibration_config.save()
+>>> deflection_channel_config = HDF5_ChannelConfig(filename=None)
+>>> deflection_channel_config['maxdata'] = 200
+>>> deflection_channel_config['conversion-coefficients'] = (0,1)
+>>> deflection_channel_config['conversion-origin'] = 0
+>>> deflection_channel_config['inverse-conversion-coefficients'] = (0,1)
+>>> deflection_channel_config['inverse-conversion-origin'] = 0
+
+We'll be generating a test vibration time series with the following
+parameters. Make sure these are all floats to avoid accidental
+integer division in later steps.
+
+>>> m = 5e-11 # kg
+>>> gamma = 1.6e-6 # N*s/m
+>>> k = 0.05 # N/m
+>>> T = 1/vibration_config['frequency']
+>>> T # doctest: +ELLIPSIS
+2.00...e-05
+>>> N = int(2**15) # count
+>>> F_sigma = 1e3 # N
+
+where `T` is the sampling period, `N` is the number of samples, and
+`F_sigma` is the standard deviation of the white-noise external force.
+Note that the resonant frequency is less than the Nyquist frequency so
+we don't have to worry too much about aliasing.
+
+>>> w0 = numpy.sqrt(k/m)
+>>> f0 = w0/(2*numpy.pi)
+>>> f0 # doctest: +ELLIPSIS
+5032.9...
+>>> f_nyquist = vibration_config['frequency']/2
+>>> f_nyquist # doctest: +ELLIPSIS
+25000.0...
+
+The damping ratio is
+
+>>> damping = gamma / (2*m*w0)
+>>> damping # doctest: +ELLIPSIS
+0.505...
+
+The quality factor is
+
+>>> Q = m*w0 / gamma
+>>> Q # doctest: +ELLIPSIS
+0.988...
+>>> (1 / (2*damping)) / Q # doctest: +ELLIPSIS
+1.000...
+
+We expect the white-noise power spectral density (PSD) to be a flat
+line at
+
+>>> F0 = F_sigma**2 * 2 * T
+
+because the integral from `0` `1/2T` should be `F_sigma**2`.
+
+The expected time series PSD parameters are
+
+>>> A = f0
+>>> B = gamma/(m*2*numpy.pi)
+>>> C = F0/(m**2*(2*numpy.pi)**4)
+
+Simulate a time series with the proper PSD using center-differencing.
+
+ m\ddot{x} + \gamma \dot{x} + kx = F
+
+ m \frac{x_{i+1} - 2x_i + x_{i-1}}{T**2}
+ + \gamma \frac{x_{i+1}-x_{i-1}}{T}
+ + kx_i = F_i
+
+ a x_{i+1} + b x_{i} + c x_{i-1} = F_i
+
+where `T` is the sampling period, `i=t/T` is the measurement index,
+`a=m/T**2+gamma/2T`, `b=-2m/T**2+k`, and `c=m/T**2-gamma/2T`.
+Rearranging and shifting to `j=i-1`
+
+ x_j = \frac{F_{i-1} - bx_{i-1} - cx_{i-2}}{a}
+
+>>> a = m/T**2 + gamma/(2*T)
+>>> b = -2*m/T**2 + k
+>>> c = m/T**2 - gamma/(2*T)
+>>> x = numpy.zeros((N+2,), dtype=numpy.float) # two extra initial points
+>>> F = numpy.zeros((N,), dtype=numpy.float)
+>>> for i in range(2, x.size):
+... Fp = random.gauss(mu=0, sigma=F_sigma)
+... F[i-2] = Fp
+... xp = x[i-1]
+... xpp = x[i-2]
+... x[i] = (Fp - b*xp - c*xpp)/a
+>>> x = x[2:] # drop the initial points
+
+Convert the simulated data to bits.
+
+>>> deflection = x
+>>> deflection_bits = convert_volts_to_bits(deflection_channel_config, x)
+
+Analyze the simulated data.
+
+>>> naive_vibration = vib_analyze_naive(deflection)
+>>> naive_vibration # doctest: +SKIP
+136871517.43486854
+>>> abs(naive_vibration / 136.9e6 - 1) < 0.1
+True
+
+>>> processed_vibration = vib_analyze(
+... deflection_bits, vibration_config, deflection_channel_config)
+>>> processed_vibration # doctest: +SKIP
+136457906.25574699
+
+>>> vib_plot(deflection=deflection_bits, vibration_config=vibration_config)
+>>> vib_save(filename=filename, group='/vibration/',
+... raw_vibration=deflection_bits, vibration_config=vibration_config,
+... deflection_channel_config=deflection_channel_config,
+... processed_vibration=processed_vibration)
+
+>>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
+/
+ /vibration
+ /vibration/config
+ /vibration/config/deflection
+ <HDF5 dataset "channel": shape (), type "|S1">
+ 0
+ <HDF5 dataset "conversion-coefficients": shape (), type "|S4">
+ 0, 1
+ <HDF5 dataset "conversion-origin": shape (), type "|S1">
+ 0
+ <HDF5 dataset "device": shape (), type "|S12">
+ /dev/comedi0
+ <HDF5 dataset "inverse-conversion-coefficients": shape (), type "|S4">
+ 0, 1
+ <HDF5 dataset "inverse-conversion-origin": shape (), type "|S1">
+ 0
+ <HDF5 dataset "maxdata": shape (), type "|S3">
+ 200
+ <HDF5 dataset "range": shape (), type "|S1">
+ 1
+ <HDF5 dataset "subdevice": shape (), type "|S2">
+ -1
+ /vibration/config/vibration
+ <HDF5 dataset "chunk-size": shape (), type "|S4">
+ 2048
+ <HDF5 dataset "frequency": shape (), type "|S7">
+ 50000.0
+ <HDF5 dataset "maximum-fit-frequency": shape (), type "|S7">
+ 25000.0
+ <HDF5 dataset "minimum-fit-frequency": shape (), type "|S5">
+ 500.0
+ <HDF5 dataset "model": shape (), type "|S12">
+ Breit-Wigner
+ <HDF5 dataset "overlap": shape (), type "|S2">
+ no
+ <HDF5 dataset "sample-time": shape (), type "|S1">
+ 1
+ <HDF5 dataset "window": shape (), type "|S4">
+ Hann
+ <HDF5 dataset "processed": shape (), type "<f8">
+ ...
+ /vibration/raw
+ <HDF5 dataset "deflection": shape (32768,), type "<f8">
+ [...]
+
+>>> (raw_vibration,vibration_config,deflection_channel_config,
+... processed_vibration) = vib_load(
+... filename=filename, group='/vibration/')
+
+>>> processed_vibration # doctest: +SKIP
+136457906.25574699
+>>> abs(processed_vibration / 136.5e6 - 1) < 0.1
+True
+
+>>> os.remove(filename)
"""
-import os, time, numpy
-import GnuplotBiDir # can be used for fitting lorentzian
-import scipy.optimize # can be used for fitting lorentzian
+import os as _os
+import time as _time
-import data_logger
-import FFT_tools
-from splittable_kwargs import splittableKwargsFunction, \
- make_splittable_kwargs_function
+import h5py as _h5py
+import numpy as _numpy
+from scipy.optimize import leastsq as _leastsq
-from . import common
-from . import config
+try:
+ import matplotlib as _matplotlib
+ import matplotlib.pyplot as _matplotlib_pyplot
+ import time as _time # for timestamping lines on plots
+except (ImportError, RuntimeError), e:
+ _matplotlib = None
+ _matplotlib_import_error = e
+import FFT_tools as _FFT_tools
+from pypiezo.base import convert_bits_to_volts as _convert_bits_to_volts
+from pypiezo.config import HDF5_ChannelConfig as _HDF5_ChannelConfig
+from pypiezo.config import h5_create_group as _h5_create_group
-PLOT_GUESSED_LORENTZIAN=False
+from . import LOG as _LOG
+from . import base_config as _base_config
+from .config import Variance as _Variance
+from .config import BreitWigner as _BreitWigner
+from .config import OffsetBreitWigner as _OffsetBreitWigner
+from .config import HDF5_VibrationConfig as _HDF5_VibrationConfig
-@splittableKwargsFunction()
-def vib_analyze_naive(deflection_bits, Vphoto_in2V=config.Vphoto_in2V,
- zeroVphoto_bits=config.zeroVphoto_bits) :
- """
- Calculate the variance of the raw data, and convert it to Volts**2.
- This method is simple and easy to understand, but it highly succeptible
- to noise, drift, etc.
+
+def vib_analyze_naive(deflection):
+ """Calculate the deflection variance in Volts**2.
+
+ This method is simple and easy to understand, but it highly
+ succeptible to noise, drift, etc.
- Vphoto_in2V is a function converting Vphoto values from bits to Volts.
- This function is assumed linear, see bump_analyze().
- """
- if config.TEXT_VERBOSE :
- print "Calculating the naive (raw) variance"
- Vphoto_std_bits = deflection_bits.std()
- Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zeroVphoto_bits)
- if config.TEXT_VERBOSE :
- print "Average deflection (Volts):", Vphoto_in2V(deflection_bits.mean())
- return Vphoto_std**2
-
-#@splittableKwargsFunction((fit_psd, 'freq_axis', 'psd_data'),
-# (vib_plot, 'deflection_bits', 'freq_axis', 'power',
-# 'A', 'B', 'C', 'minFreq', 'maxFreq'))
-# Some of the child functions aren't yet defined, so postpone
-# make-splittable until later in the module.
-def vib_analyze(deflection_bits, freq,
- chunk_size=4096, overlap=True, window=FFT_tools.window_hann,
- Vphoto_in2V=config.Vphoto_in2V, **kwargs) :
+ Inputs
+ deflection : numpy array with deflection timeseries in Volts.
"""
- Calculate the variance in the raw data due to the cantilever's
- thermal oscillation and convert it to Volts**2.
+ std = deflection.std()
+ var = std**2
+ _LOG.debug('naive deflection variance: %g V**2' % var)
+ return var
- Improves on vib_analyze_naive() by first converting Vphoto(t) to
- frequency space, and fitting a Lorentzian in the relevant
+def vib_analyze(deflection, vibration_config, deflection_channel_config,
+ plot=False):
+ """Calculate the deflection variance in Volts**2.
+
+ Improves on `vib_analyze_naive()` by first converting `Vphoto(t)`
+ to frequency space, and fitting a Breit-Wigner in the relevant
frequency range (see cantilever_calib.pdf for derivation).
However, there may be cases where the fit is thrown off by noise
spikes in frequency space. To protect from errors, the fitted
variance is compared to the naive variance (where all noise is
included), and the minimum variance is returned.
- Vphoto_in2V is a function converting Vphoto values from bits to Volts.
- This function is assumed linear, see bump_analyze().
+ Inputs:
+ deflection Vphoto deflection input in bits.
+ vibration_config `.config._VibrationConfig` instance
+ deflection_channel_config
+ deflection `pypiezo.config._ChannelConfig` instance
+ plot boolean overriding matplotlib config setting.
+
+ The conversion to frequency space generates an average power
+ spectrum by breaking the data into windowed chunks and averaging
+ the power spectrums for the chunks together. See
+ `FFT_tools.unitary_avg_power_spectrum()` for details.
"""
- fit_psd_kwargs,vib_plot_kwargs = vib_analyze._splitargs(vib_analyze,kwargs)
- if 'minFreq' in fit_psd_kwargs:
- vib_plot_kwargs['minFreq'] = fit_psd_kwargs['minFreq']
- if 'maxFreq' in fit_psd_kwargs:
- vib_plot_kwargs['maxFreq'] = fit_psd_kwargs['maxFreq']
-
- Vphoto_t = numpy.zeros((len(deflection_bits),),
- dtype=numpy.float)
# convert the data from bits to volts
- if config.TEXT_VERBOSE :
- print "Converting %d bit values to voltages" % len(Vphoto_t)
- for i in range(len(deflection_bits)) :
- Vphoto_t[i] = Vphoto_in2V(deflection_bits[i])
- if config.TEXT_VERBOSE :
- print "Average deflection (Volts):", Vphoto_t.mean()
+ deflection_v = _convert_bits_to_volts(
+ deflection_channel_config, deflection)
+ mean = deflection_v.mean()
+ _LOG.debug('average thermal deflection (Volts): %g' % mean)
+
+ naive_variance = vib_analyze_naive(deflection_v)
+ if vibration_config['model'] == _Variance:
+ return naive_variance
# Compute the average power spectral density per unit time (in uV**2/Hz)
- if config.TEXT_VERBOSE :
- print "Compute the averaged power spectral density in uV**2/Hz"
- (freq_axis, power) = \
- FFT_tools.unitary_avg_power_spectrum((Vphoto_t - Vphoto_t.mean())*1e6,
- freq, chunk_size,overlap, window)
+ _LOG.debug('compute the averaged power spectral density in uV**2/Hz')
+ freq_axis,power = _FFT_tools.unitary_avg_power_spectrum(
+ (deflection_v - mean)*1e6, vibration_config['frequency'],
+ vibration_config['chunk-size'], vibration_config['overlap'],
+ vibration_config['window'])
- A,B,C,D = fit_psd(freq_axis, power, **kwargs)
+ A,B,C,D = fit_psd(
+ freq_axis, power,
+ min_frequency=vibration_config['minimum-fit-frequency'],
+ max_frequency=vibration_config['maximum-fit-frequency'],
+ offset=vibration_config['model'] == _OffsetBreitWigner)
- if config.TEXT_VERBOSE :
- print "Fitted f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with"
- print "A = %g, \t B = %g, \t C = %g, \t D = %g" % (A, B, C, D)
+ _LOG.debug('fit PSD(f) = C / ((A**2-f**2)**2 + (f*B)**2) with '
+ 'A = %g, B = %g, C = %g, D = %g' % (A, B, C, D))
- vib_plot(deflection_bits, freq_axis, power, A, B, C, D, **kwargs)
+ if plot or _base_config['matplotlib']:
+ vib_plot(deflection, freq_axis, power, A, B, C, D,
+ vibration_config=vibration_config)
# Our A is in uV**2, so convert back to Volts**2
- fitted_variance = lorentzian_area(A,B,C) * 1e-12
+ fitted_variance = breit_wigner_area(A,B,C) * 1e-12
+
+ _LOG.debug('fitted deflection variance: %g V**2' % fitted_variance)
+
+ if _base_config['matplotlib']:
+ vib_plot(deflection, freq_axis, power, A, B, C, D,
+ vibration_config=vibration_config)
- naive_variance = vib_analyze_naive(deflection_bits,
- Vphoto_in2V=Vphoto_in2V)
- if config.TEXT_VERBOSE:
- print "Fitted variance: %g V**2" % fitted_variance
- print "Naive variance : %g V**2" % naive_variance
-
return min(fitted_variance, naive_variance)
-@splittableKwargsFunction()
-def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=25000) :
- """
- freq_axis : array of frequencies in Hz
- psd_data : array of PSD amplitudes for the frequencies in freq_axis
+def breit_wigner(f, A, B, C, D=0):
+ """Breit-Wigner (sortof).
- Calculate the variance in the raw data due to the cantilever's
- thermal oscillation and convert it to Volts**2.
+ Inputs
+ f Frequency
+ A Resonant frequency
+ B Quality Q = A/B
+ C Scaling factor
+ D Optional white-noise offset
- The conversion to frequency space generates an average power spectrum
- by breaking the data into windowed chunks and averaging the power
- spectrums for the chunks together.
- See FFT_tools.unitary_avg_power_spectrum().
+ All parameters must be postive.
+ """
+ return abs(C) / ((A**2-f**2)**2 + (B*f)**2) + abs(D)
+
+def fit_psd(freq_axis, psd_data, min_frequency=500, max_frequency=25000,
+ offset=False):
+ """Fit the FFTed vibration data to a Breit-Wigner.
+
+ Inputs
+ freq_axis array of frequencies in Hz
+ psd_data array of PSD amplitudes for the frequencies in freq_axis
+ min_frequency lower bound of Breit-Wigner fitting region
+ max_frequency upper bound of Breit-Wigner fitting region
+ offset add a white-noise offset to the Breit-Wigner model
+ Output
+ Breit-Wigner model fit parameters `A`, `B`, `C`, and `D`.
"""
# cut out the relevent frequency range
- if config.TEXT_VERBOSE :
- print "Cut the frequency range %g to %g Hz" % (minFreq, maxFreq)
+ _LOG.debug('cut the frequency range %g to %g Hz'
+ % (min_frequency, max_frequency))
imin = 0
- while freq_axis[imin] < minFreq : imin += 1
+ while freq_axis[imin] < min_frequency : imin += 1
imax = imin
- while freq_axis[imax] < maxFreq : imax += 1
- assert imax >= imin + 10 , \
- "Less than 10 points in freq range (%g,%g)" % (minFreq, maxFreq)
+ while freq_axis[imax] < max_frequency : imax += 1
+ assert imax >= imin + 10 , 'less than 10 points in freq range (%g,%g)' % (
+ min_frequency, max_frequency)
- # generate guesses for Lorentzian parameters A,B,C
- max_psd_index = numpy.argmax(psd_data[imin:imax]) + imin
+ # generate guesses for Breit-Wigner parameters A, B, C, and D
+ max_psd_index = _numpy.argmax(psd_data[imin:imax]) + imin
max_psd = psd_data[max_psd_index]
- half_max_index = imin
- while psd_data[half_max_index] < max_psd/2 :
- half_max_index += 1
res_freq = freq_axis[max_psd_index]
- half_width = freq_axis[max_psd_index] - freq_axis[half_max_index]
- # Lorentzian L(x) = C / ((A**2-x**2)**2 + (B*x)**2)
+
+ # Breit-Wigner L(x) = C / ((A**2-x**2)**2 + (B*x)**2)
# is expected power spectrum for
# x'' + B x' + A^2 x'' = F_external(t)/m
# (A = omega_0)
# C = (2 k_B T B) / (pi m)
#
# A = resonant frequency
- # peak at x_res = sqrt(A^2 - B^2/2)
+ # peak at x_res = sqrt(A^2 - B^2/2) (by differentiating)
# which could be complex if there isn't a peak (overdamped)
# peak height = C / (3 x_res^4 + A^4)
+ #
# Q = A/B
#
- # guess Q = 1, and adjust from there
+ # Guessing Q = 1 is pretty safe.
+
Q_guess = 1
+
# so x_res^2 = B^2 Q^2 - B^2/2 = (Q^2-1/2)B^2
# B = x_res / sqrt(Q^2-1/2)
- if config.TEXT_VERBOSE :
- print "params : %g, %g" % (res_freq, max_psd)
- B_guess = res_freq / numpy.sqrt(Q_guess**2-0.5)
+ B_guess = res_freq / _numpy.sqrt(Q_guess**2-0.5)
A_guess = Q_guess*B_guess
C_guess = max_psd * (-res_freq**4 + A_guess**4)
- D_guess = 0 # allow a noise floor offset (DISABLED in fitting below)
- #
+ if offset:
+ D_guess = psd_data[max_psd_index]
+ C_guess -= D_guess
+ else:
+ D_guess = 0
+ _LOG.debug(('guessed params: resonant freq %g, max psd %g, Q %g, A %g, '
+ 'B %g, C %g, D %g') % (
+ res_freq, max_psd, Q_guess, A_guess, B_guess, C_guess, D_guess))
# Half width w on lower side when L(a-w) = L(a)/2
# (a**2 - (a-w)**2)**2 + (b*(a-w))**2 = 2*(b*a)**2
# Let W=(a-w)**2, A=a**2, and B=b**2
# The mass m is given by m = k_B T / (pi A)
# The spring constant k is given by k = m (omega_0)**2
# The quality factor Q is given by Q = omega_0 m / gamma
- if config.TEXT_VERBOSE :
- print "guesses : %g, %g, %g" % (A_guess, B_guess, C_guess)
- if PLOT_GUESSED_LORENTZIAN:
- vib_plot(None, freq_axis, psd_data, A_guess, B_guess, C_guess,
- minFreq, maxFreq, plotVerbose=True)
# Fitting the PSD of V = photoSensitivity*x just rescales the parameters
- if False: # Gnuplot fit worse than scipy.optimize.leastsq fit.
- # fit Lorentzian using Gnuplot's 'fit' command
-
- # save about-to-be-fitted data to a temporary file
- datafilename = "%s_%d" % (config.GNUFIT_DATA_BASE, time.time())
- fd = file(datafilename, 'w')
- for i in range(imin, imax) :
- fd.write("%g\t%g\n" % (freq_axis[i], psd_data[i]))
- fd.close()
-
- g = GnuplotBiDir.Gnuplot()
- g("f(x) = C / ((A**2-x**2)**2 + (B*x)**2)")
- g("A=%g; B=%g; C=%g" % (A_guess, B_guess, C_guess))
- g("fit f(x) '%s' via A,B,C" % datafilename)
- A=abs(float(g.getvar('A'))) # A and B only show up as squares in f(x)
- B=abs(float(g.getvar('B'))) # so ensure we get positive values
- C=float(g.getvar('C'))
-
- os.remove(datafilename)
+ # fit Breit-Wigner using scipy.optimize.leastsq
+ def residual(p, y, x):
+ return breit_wigner(x, *p) - y
+ if offset:
+ guess = _numpy.array((A_guess, B_guess, C_guess, D_guess))
+ else:
+ guess = _numpy.array((A_guess, B_guess, C_guess))
+
+ p,cov,info,mesg,ier = _leastsq(
+ residual, guess,
+ args=(psd_data[imin:imax], freq_axis[imin:imax]),
+ full_output=True, maxfev=10000)
+ _LOG.debug('fitted params: %s' % p)
+ _LOG.debug('covariance matrix: %s' % cov)
+ #_LOG.debug('info: %s' % info)
+ _LOG.debug('message: %s' % mesg)
+ if ier == 1:
+ _LOG.debug('solution converged')
+ else:
+ _LOG.debug('solution did not converge')
+ if offset:
+ A,B,C,D = p
else:
- # fit Lorenzian using scipy.optimize.leastsq
- def residual(p, y, x):
- A = p[0]
- B = p[1]
- C = p[2]
- Y = C / ((A**2-x**2)**2 + (B*x)**2)
- return Y-y
- leastsq = scipy.optimize.leastsq
- p,cov,info,mesg,ier = leastsq(residual, (A_guess, B_guess, C_guess),
- args=(psd_data[imin:imax],
- freq_axis[imin:imax]),
- full_output=True, maxfev=10000)
- if config.TEXT_VERBOSE:
- print "Fitted params:",p
- print "Covariance mx:",cov
- #print "Info:", info # too much information :p
- print "mesg:", mesg
- if ier == 1 :
- print "Solution converged"
- else :
- print "Solution did not converge"
A,B,C = p
- D = 0 # do not fit the noise floor (fit doesn't look very convincing)
- A=abs(A) # A and B only show up as squares in f(x)
- B=abs(B) # so ensure we get positive values
+ D = 0
+ A=abs(A) # A and B only show up as squares in f(x)
+ B=abs(B) # so ensure we get positive values
return (A, B, C, D)
-def lorentzian_area(A, B, C) :
+def breit_wigner_area(A, B, C):
# Integrating the the power spectral density per unit time (power) over the
# frequency axis [0, fN] returns the total signal power per unit time
# int_0^fN power(f)df = 1/T int_0^T |x(t)**2| dt
# = <V(t)**2>, the variance for our AC signal.
- # The variance from our fitted Lorentzian is the area under the Lorentzian
+ # The variance from our fitted Breit-Wigner is the area under the Breit-Wigner
# <V(t)**2> = (pi*C) / (2*B*A**2)
- return (numpy.pi * C) / (2 * B * A**2)
+ return (_numpy.pi * C) / (2 * B * A**2)
+
+def vib_save(filename, group='/', raw_vibration=None, vibration_config=None,
+ deflection_channel_config=None, processed_vibration=None):
+ with _h5py.File(filename, 'a') as f:
+ cwg = _h5_create_group(f, group)
+ if raw_vibration is not None:
+ try:
+ del cwg['raw/deflection']
+ except KeyError:
+ pass
+ cwg['raw/deflection'] = raw_vibration
+ for config,key in [(vibration_config, 'config/vibration'),
+ (deflection_channel_config,
+ 'config/deflection')]:
+ if config is None:
+ continue
+ config_cwg = _h5_create_group(cwg, key)
+ config.save(group=config_cwg)
+ if processed_vibration is not None:
+ try:
+ del cwg['processed']
+ except KeyError:
+ pass
+ cwg['processed'] = processed_vibration
+
+def vib_load(filename, group='/'):
+ assert group.endswith('/')
+ raw_vibration = processed_vibration = None
+ configs = []
+ with _h5py.File(filename, 'a') as f:
+ try:
+ raw_vibration = f[group+'raw/deflection'][...]
+ except KeyError:
+ pass
+ for Config,key in [(_HDF5_VibrationConfig, 'config/vibration'),
+ (_HDF5_ChannelConfig, 'config/deflection')]:
+ config = Config(filename=filename, group=group+key)
+ configs.append(config)
+ try:
+ processed_vibration = f[group+'processed'][...]
+ except KeyError:
+ pass
+ ret = [raw_vibration]
+ ret.extend(configs)
+ ret.append(processed_vibration)
+ for config in configs:
+ config.load()
+ return tuple(ret)
+
+def vib_plot(deflection=None, freq_axis=None, power=None, A=None, B=None,
+ C=None, D=0, vibration_config=None, analyze=False):
+ """Plot 3 subfigures displaying vibration data and analysis.
-def gnuplot_check_fit(datafilename, A, B, C) :
- """
- return a string containing a gnuplot script displaying the fit.
- """
- string = 'f(x) = C / ((A**2-x**2)**2 + (x*B)**2)\n'
- string += 'A= %g ; B= %g ; C= %g\n' % (A, B, C)
- string += 'set logscale y\n'
- string += "plot '%s', f(x)\n" % datafilename
- return string
-
-@splittableKwargsFunction()
-def vib_save(data, log_dir=None) :
- """Save the dictionary data, using data_logger.data_log()
- data should be a dictionary of arrays with the fields
- 'Z piezo output'
- 'Z piezo input'
- 'Deflection input'
- """
- if log_dir :
- freq = data['sample frequency Hz']
- log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
- log_name="vib_%gHz" % freq)
- log.write_dict_of_arrays(data)
-
-def vib_load(datafile=None) :
- """Load the dictionary data, using data_logger.date_load()"
- data should be a dictionary of arrays with the fields
- 'Z piezo output'
- 'Z piezo input'
- 'Deflection input'
- """
- dl = data_logger.data_load()
- data = dl.read_dict_of_arrays(datafile)
- return data
-
-@splittableKwargsFunction()
-def vib_plot(deflection_bits, freq_axis, power, A, B, C, D,
- minFreq=None, maxFreq=None, plotVerbose=False) :
- """
- If plotVerbose or config.PYLAB_VERBOSE == True, plot 3 subfigures:
Time series (Vphoto vs sample index) (show any major drift events),
A histogram of Vphoto, with a gaussian fit (demonstrate randomness), and
FFTed Vphoto data (Vphoto vs Frequency) (show major noise components).
"""
- if plotVerbose or config.PYLAB_VERBOSE :
- common._import_pylab()
- common._pylab.figure(config.BASE_FIGNUM+2)
-
- if deflection_bits != None:
- # plot time series
- common._pylab.subplot(311)
- common._pylab.hold(False)
- common._pylab.plot(deflection_bits, 'r.')
- common._pylab.title("free oscillation")
-
- # plot histogram distribution and gaussian fit
- common._pylab.subplot(312)
- common._pylab.hold(False)
- n, bins, patches = \
- common._pylab.hist(deflection_bits, bins=30,
- normed=1, align='center')
- gaus = numpy.zeros((len(bins),), dtype=numpy.float)
- mean = deflection_bits.mean()
- std = deflection_bits.std()
- pi = numpy.pi
- exp = numpy.exp
- for i in range(len(bins)) :
- gaus[i] = (2*pi)**(-.5)/std * exp(-0.5*((bins[i]-mean)/std)**2)
- common._pylab.hold(True)
- common._pylab.plot(bins, gaus, 'r-');
- common._pylab.hold(False)
-
- # plot FFTed data
- axes = common._pylab.subplot(313)
- else:
- # use a nice big subplot just for the FFTed data
- axes = common._pylab.subplot(111)
- common._pylab.hold(False)
- common._pylab.semilogy(freq_axis, power, 'r.-')
- common._pylab.hold(True)
- xmin,xmax = axes.get_xbound()
- ymin,ymax = axes.get_ybound()
-
- if minFreq is not None and maxFreq is not None:
- # highlight the region we're fitting in
- common._pylab.axvspan(minFreq, maxFreq, facecolor='g', alpha=0.1,
- zorder = -1)
-
- fitdata = C / ((A**2-freq_axis**2)**2 + (B*freq_axis)**2) + D
- common._pylab.plot(freq_axis, fitdata, 'b-')
- noisefloor = D + 0*freq_axis;
- common._pylab.plot(freq_axis, noisefloor)
- common._pylab.hold(False)
- axes.axis([xmin,xmax,ymin,ymax])
-
- common._flush_plot()
- if (plotVerbose or config.GNUPLOT_VERBOSE) and False: # TODO: cleanup and test
- # write all the ft data now
- fd = file(datafilename, 'w')
- for i in range(len(freq_axis)) :
- fd.write("%g\t%g\n" % (freq_axis[i], power[i]))
- fd.write("\n") # blank line to separate fit data.
- # write the fit data
- for i in range(imin, imax) :
- x = freq_axis[i]
- fd.write("%g\t%g\n" % (freq_axis[i],
- C / ((A**2-x**2)**2 + (x*B)**2) ) )
- fd.close()
- print Vphoto_t.var(), A, B, C, ((numpy.pi * C)/(2 * B * A**2) * 1e-12)
- g("set terminal epslatex color solid")
- g("set output 'calibration.tex'")
- g("set size 2,2") # multipliers on default 5"x3.5"
- #g("set title 'Thermal calibration'")
- g("set logscale y")
- g("set xlabel 'Frequency (Hz)'")
- g("set ylabel 'Power ($\mu$V$^2$/Hz)'")
- g("set xrange [0:20000]")
- g("g(x) = x < %g || x > %g ? 1/0 : f(x)" % (minFreq, maxFreq)) # only plot on fitted region
- g("plot '%s' title 'Vibration data', g(x) title 'Lorentzian fit'" % datafilename)
-
- g("set mouse")
- g("pause mouse 'click with mouse'")
- g.getvar("MOUSE_BUTTON")
- del(g)
-
-# Make vib_analyze splittable (was postponed until the child functions
-# were defined).
-make_splittable_kwargs_function(vib_analyze,
- (fit_psd, 'freq_axis', 'psd_data'),
- (vib_plot, 'deflection_bits', 'freq_axis', 'power', 'A', 'B', 'C', 'D',
- 'minFreq', 'maxFreq'))
-
-@splittableKwargsFunction((vib_analyze, 'deflection_bits', 'freq'))
-def vib_load_analyze_tweaked(tweak_file, analyze=vib_analyze, **kwargs) :
- """
- Read the vib files listed in tweak_file.
- Return an array of Vphoto variances (due to the cantilever) in Volts**2
- """
- analyze_kwargs, = \
- vib_load_analyze_tweaked._splitargs(vib_load_analyze_tweaked, kwargs)
- dl = data_logger.data_load()
- Vphoto_var = []
- for path in file(tweak_file, 'r') :
- path = path.strip()
- if path[0] == '#': # a comment
- continue
- # read the data
- data = vib_load(path)
- deflection_bits = data['Deflection input']
- freq = float(path.split('_')[-1].split('Hz')[0])
- if config.TEXT_VERBOSE :
- print "Analyzing '%s' at %g Hz" % (path, freq)
- # analyze
- if 'freq' in analyze._kwargs(analyze):
- var = analyze(deflection_bits, freq, **analyze_kwargs)
- else:
- var = analyze(deflection_bits, **analyze_kwargs)
- Vphoto_var.append(var)
- return numpy.array(Vphoto_var, dtype=numpy.float)
-
-# commandline interface functions
-import scipy.io, sys
-
-def read_data(ifile):
- """
- ifile can be a filename string or open (seekable) file object.
- returns an array of data values (1 column)
- """
- #returns (column 1 array, column 2 array)
- #"""
- if ifile == None : ifile = sys.stdin
- unlabeled_data=scipy.io.read_array(ifile)
- return unlabeled_data #(unlabeled_data[:,0], unlabeled_data[:,1])
-
-if __name__ == '__main__' :
- # command line interface
- from optparse import OptionParser
-
- usage_string = ('%prog <input-file>\n'
- '2007-2009 W. Trevor King.\n'
- '\n'
- 'There are several operation modes, each handling a different <input-file>\n'
- 'type. In each case, the output is the fitted variance (in square Volts).\n'
- '\n'
- 'Single file mode without datalogger mode (the default) :\n'
- '<input-file> should be a 1 column ASCII without a header line of cantilever\n'
- 'deflection (in bits)\n'
- 'e.g: "<deflection bits>\\n..."\n'
- '\n'
- 'Single file mode with datalogger mode :\n'
- 'Same as without datalogger mode, except the input should be a datalogger file\n'
- 'with data stored in bits (e.g. as saved by the unfold module).\n'
- '\n'
- 'Tweak file mode:\n'
- 'Runs the same analysis as in single file mode for each vib file in a tweak\n'
- 'file. Each line in the tweak file specifies a single vib file. Blank lines\n'
- 'and those beginning with a pound sign (#) are ignored. Vib files are valid\n'
- 'datalogger files as per single-file-with-datalogger-mode.\n'
- )
- parser = OptionParser(usage=usage_string, version='%prog 0.1')
- parser.add_option('-f', '--sample-freq', dest='freq', default=100e3,
- help='sample frequency in Hz (default %default)',
- type='float', metavar='FREQ')
- parser.add_option('-m', '--min-freq', dest='min_freq', default=500,
- help='minimum frequency in Hz (default %default)',
- type='float', metavar='FREQ')
- parser.add_option('-M', '--max-freq', dest='max_freq', default=7000,
- help='maximum frequency in Hz (default %default)',
- type='float', metavar='FREQ')
- parser.add_option('-o', '--output-file', dest='ofilename',
- help='write output to FILE (default stdout)',
- type='string', metavar='FILE')
- parser.add_option('-c', '--comma-out', dest='comma_out', action='store_true',
- help='Output comma-seperated values (default %default)',
- default=False)
- parser.add_option('-g', '--gnuplot', dest='gnuplot', action='store_true',
- help='Print gnuplot fit check script to stderr',
- default=False)
- parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
- help='Produce pylab fit checks during execution',
- default=False)
- parser.add_option('-G', '--plot-guess', dest='plot_guess', action='store_true',
- help='Produce plots of the pre-fit Lorentzian',
- default=False)
- parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
- help='Run in tweak-file mode',
- default=False)
- parser.add_option('-d', '--datalogger-mode', dest='datalogger_mode', action='store_true',
- help='Read input files with datalogger.read_dict_of_arrays(). This is useful, for example, to test a single line from a tweakfile.',
- default=False)
- parser.add_option('-s', '--disable-spectrum-fitting', dest='spectrum_fitting',
- action='store_false',
- help='Disable PSD fitting, just use the raw variance',
- default=True)
- parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
- help='Print lots of debugging information',
- default=False)
-
- options,args = parser.parse_args()
- parser.destroy()
- assert len(args) >= 1, "Need an input file"
-
- ifilename = args[0]
-
- if options.ofilename != None :
- ofile = file(options.ofilename, 'w')
- else :
- ofile = sys.stdout
- config.TEXT_VERBOSE = options.verbose
- config.PYLAB_INTERACTIVE = False
- config.PYLAB_VERBOSE = options.pylab
- config.GNUPLOT_VERBOSE = options.gnuplot
- PLOT_GUESSED_LORENTZIAN = options.plot_guess
- if options.spectrum_fitting == True:
- analyze = vib_analyze
- analyze_args = {'freq':options.freq,
- 'minFreq':options.min_freq,
- 'maxFreq':options.max_freq}
+ if not _matplotlib:
+ raise _matplotlib_import_error
+ figure = _matplotlib_pyplot.figure()
+
+ if power is None:
+ assert deflection != None, (
+ 'must set at least one of `deflection` or `power`.')
+ time_axes = figure.add_subplot(2, 1, 1)
+ hist_axes = figure.add_subplot(2, 1, 2)
+ freq_axes = None
+ elif deflection is None:
+ time_axes = hist_axes = None
+ freq_axes = figure.add_subplot(1, 1, 1)
else:
- analyze = vib_analyze_naive
- analyze_args = dict()
- make_splittable_kwargs_function(vib_load_analyze_tweaked,
- (vib_analyze, 'deflection_bits'))
-
- if options.tweakmode == False : # single file mode
- if options.datalogger_mode:
- data = vib_load(ifilename)
- deflection_bits = data['Deflection input']
- else:
- deflection_bits = read_data(ifilename)
- Vphoto_var = analyze(deflection_bits, **analyze_args)
- print >> ofile, Vphoto_var
- else : # tweak mode
- if 'freq' in analyze_args:
- analyze_args.pop('freq') # freq extracted from vib file names
- Vphoto_vars = vib_load_analyze_tweaked(ifilename, analyze=analyze,
- **analyze_args)
- if options.comma_out :
- sep = ','
- else :
- sep = '\n'
- common.write_array(ofile, Vphoto_vars, sep)
-
- if common._final_flush_plot != None:
- common._final_flush_plot()
+ time_axes = figure.add_subplot(3, 1, 1)
+ hist_axes = figure.add_subplot(3, 1, 2)
+ freq_axes = figure.add_subplot(3, 1, 3)
+
+ timestamp = _time.strftime('%H%M%S')
+
+ if deflection is not None:
+ time_axes.plot(deflection, 'r.')
+ time_axes.set_title('free oscillation')
+
+ # plot histogram distribution and gaussian fit
+ hist_axes.hold(True)
+ n,bins,patches = hist_axes.hist(
+ deflection, bins=30, normed=True, align='mid')
+ gauss = _numpy.zeros((len(bins),), dtype=_numpy.float)
+ mean = deflection.mean()
+ std = deflection.std()
+ pi = _numpy.pi
+ exp = _numpy.exp
+ gauss = _numpy.sqrt(2*pi)/std * exp(-0.5*((bins-mean)/std)**2)
+ # Matplotlib's normed histogram uses bin heights of n/(len(x)*dbin)
+ dbin = bins[1]-bins[0]
+ hist_axes.plot(bins, gauss/dbin, 'r-')
+ if power is not None:
+ freq_axes.hold(True)
+ freq_axes.set_yscale('log')
+ freq_axes.plot(freq_axis, power, 'r.-')
+ xmin,xmax = freq_axes.get_xbound()
+ ymin,ymax = freq_axes.get_ybound()
- if options.ofilename != None :
- ofile.close()
-
-# LocalWords: calibcant AFM
+ # highlight the region we're fitting
+ freq_axes.axvspan(
+ vibration_config['minimum-fit-frequency'],
+ vibration_config['maximum-fit-frequency'],
+ facecolor='g', alpha=0.1, zorder=-2)
+
+ if A is not None:
+ fitdata = breit_wigner(freq_axis, A, B, C, D)
+ freq_axes.plot(freq_axis, fitdata, 'b-')
+ noisefloor = D + 0*freq_axis;
+ freq_axes.plot(freq_axis, noisefloor)
+
+ if B**2 < 2*A**2:
+ res_freq = _numpy.sqrt(A**2 - B**2/2)
+ freq_axes.axvline(res_freq, color='b', zorder=-1)
+
+ freq_axes.set_title('power spectral density %s' % timestamp)
+ freq_axes.axis([xmin,xmax,ymin,ymax])
+ freq_axes.set_xlabel('frequency (Hz)')
+ freq_axes.set_ylabel('PSD')
+
+ figure.show()
"calibcant: tools for thermally calibrating AFM cantilevers"
+package_name = 'calibcant'
classifiers = """\
Development Status :: 2 - Pre-Alpha
Intended Audience :: Developers
packages = find_packages()
-setup(name='calibcant',
+setup(name=package_name,
version=__version__,
maintainer='W. Trevor King',
maintainer_email='wking@drexel.edu',
- url = 'http://www.physics.drexel.edu/~wking/code/python/',
- download_url = 'http://www.physics.drexel.edu/~wking/code/python/calibcant-%s.tar.gz' % __version__,
- license = 'GNU General Public License (GPL)',
- platforms = ['all'],
- description = __doc__,
- long_description = open('README', 'r').read(),
- classifiers = filter(None, classifiers.split('\n')),
- packages = packages,
- provides = ['calibcant (%s)' % __version__],
- requires = ['piezo (>= 0.3)'],
+ url='http://www.physics.drexel.edu/~wking/unfolding-disasters/posts/%s/' % package_name,
+ download_url='http://www.physics.drexel.edu/~wking/code/python/%s-%s.tar.gz' % (package_name, __version__),
+ license='GNU General Public License (GPL)',
+ platforms=['all'],
+ description=__doc__,
+ long_description=open('README', 'r').read(),
+ classifiers=filter(None, classifiers.split('\n')),
+ packages=packages,
+ provides=['calibcant (%s)' % __version__],
+ requires=['pypiezo (>= 0.4)'],
)
projects / calibcant.git / commitdiff