# calibcant - tools for thermally calibrating AFM cantilevers
#
-# Copyright (C) 2008-2011 W. Trevor King <wking@drexel.edu>
+# Copyright (C) 2008-2012 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 free software: you can redistribute it and/or modify it under
+# the terms of the GNU 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.
+# 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 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/>.
+# You should have received a copy of the GNU General Public 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 h5config.storage.hdf5 import HDF5_Storage, pprint_HDF5
+>>> from .config import CalibrationConfig
+
+>>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
+>>> os.close(fd)
+
+>>> calibration_config = CalibrationConfig(storage=HDF5_Storage(
+... filename=filename, group='/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 "bump": shape (), type "|S1">
+<BLANKLINE>
+ <HDF5 dataset "num-bumps": shape (), type "<i4">
+ 10
+ <HDF5 dataset "num-temperatures": shape (), type "<i4">
+ 10
+ <HDF5 dataset "num-vibrations": shape (), type "<i4">
+ 20
+ <HDF5 dataset "temperature": shape (), type "|S1">
+<BLANKLINE>
+ <HDF5 dataset "temperature-sleep": shape (), type "<i4">
+ 1
+ <HDF5 dataset "vibration": shape (), type "|S1">
+<BLANKLINE>
+ <HDF5 dataset "vibration-spacing": shape (), type "<f8">
+ 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 h5config.storage.hdf5 import HDF5_Storage as _HDF5_Storage
+from h5config.storage.hdf5 import h5_create_group as _h5_create_group
+from . import LOG as _LOG
+from . import package_config as _package_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 CalibrationConfig as _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 _package_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')
+ storage = _HDF5_Storage()
+ storage.save(config=calibration_config, 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 = float(f[group+'processed/spring-constant/data'][...])
+ except KeyError:
+ pass
+ try:
+ k_s = float(
+ f[group+'processed/spring-constant/standard-deviation'][...])
+ except KeyError:
+ pass
+ calibration_config = _CalibrationConfig(storage=_HDF5_Storage(
+ 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)')
+
+ if hasattr(figure, 'show'):
+ 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_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_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/')
+ if bumps is None:
+ bumps = _numpy.zeros(
+ (calibration_config['num-bumps'],), dtype=float)
+ if temperatures is None:
+ temperatures = _numpy.zeros(
+ (calibration_config['num-temperatures'],), dtype=float)
+ if vibrations is None:
+ vibrations = _numpy.zeros(
+ (calibration_config['num-vibrations'],), dtype=float)
+ 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_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_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_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))
projects / calibcant.git / blobdiff