calibcant/T_analyze.py

   1 # calibcant - tools for thermally calibrating AFM cantilevers
   2 #
   3 # Copyright (C) 2008-2011 W. Trevor King <wking@drexel.edu>
   4 #
   5 # This file is part of calibcant.
   6 #
   7 # calibcant is free software: you can redistribute it and/or
   8 # modify it under the terms of the GNU Lesser General Public
   9 # License as published by the Free Software Foundation, either
  10 # version 3 of the License, or (at your option) any later version.
  11 #
  12 # calibcant is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU Lesser General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU Lesser General Public
  18 # License along with calibcant.  If not, see
  19 # <http://www.gnu.org/licenses/>.
  20
  21 """Temperature analysis.
  22
  23 Separate the more general `T_analyze()` from the other `T_*()`
  24 functions in calibcant.
  25
  26 The relevant physical quantities are:
  27
  28 * `T` Temperature at which thermal vibration measurements were acquired
  29
  30 >>> import os
  31 >>> import tempfile
  32 >>> import numpy
  33 >>> from .config import TemperatureConfig
  34 >>> from h5config.storage.hdf5 import pprint_HDF5, HDF5_Storage
  35
  36 >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
  37 >>> os.close(fd)
  38
  39 >>> temperature_config = TemperatureConfig(storage=HDF5_Storage(
  40 ...         filename=filename, group='/T/config/'))
  41
  42 >>> raw_T = numpy.array([22, 23.5, 24])
  43 >>> processed_T = T_analyze(raw_T, temperature_config)
  44 >>> T_plot(raw_T=raw_T, processed_T=processed_T)
  45 >>> T_save(filename=filename, group='/T/', raw_T=raw_T,
  46 ...     temperature_config=temperature_config, processed_T=processed_T)
  47
  48 >>> pprint_HDF5(filename)  # doctest: +REPORT_UDIFF
  49 /
  50   /T
  51     /T/config
  52       <HDF5 dataset "default": shape (), type "|b1">
  53         True
  54       <HDF5 dataset "units": shape (), type "|S7">
  55         Celsius
  56     <HDF5 dataset "processed": shape (3,), type "<f8">
  57       [ 295.15  296.65  297.15]
  58     <HDF5 dataset "raw": shape (3,), type "<f8">
  59       [ 22.   23.5  24. ]
  60
  61 >>> raw_T,temperature_config,processed_T = T_load(
  62 ...     filename=filename, group='/T/')
  63 >>> print temperature_config.dump()
  64 units: Celsius
  65 default: yes
  66 >>> raw_T
  67 array([ 22. ,  23.5,  24. ])
  68 >>> type(raw_T)
  69 <type 'numpy.ndarray'>
  70 >>> processed_T
  71 array([ 295.15,  296.65,  297.15])
  72
  73 >>> os.remove(filename)
  74 """
  75
  76 import h5py as _h5py
  77 from scipy.constants import C2K as _C2K
  78
  79 try:
  80     import matplotlib as _matplotlib
  81     import matplotlib.pyplot as _matplotlib_pyplot
  82     import time as _time  # for timestamping lines on plots
  83 except (ImportError, RuntimeError), e:
  84     _matplotlib = None
  85     _matplotlib_import_error = e
  86
  87 from h5config.storage.hdf5 import HDF5_Storage as _HDF5_Storage
  88 from h5config.storage.hdf5 import h5_create_group as _h5_create_group
  89
  90 from . import LOG as _LOG
  91 from . import package_config as _package_config
  92 from .config import Celsius as _Celsius
  93 from .config import Kelvin as _Kelvin
  94 from .config import TemperatureConfig as _TemperatureConfig
  95
  96
  97 def T_analyze(T, temperature_config):
  98     """Convert measured temperature to Kelvin.
  99
 100     `T` should be a numpy ndarray or scalar.  `temperature_config`
 101     should be a `config._TemperatureConfig` instance.
 102     """
 103     if temperature_config['units'] == _Celsius:
 104         return _C2K(T)
 105     else:
 106         return T
 107
 108 def T_save(filename, group='/', raw_T=None, temperature_config=None,
 109            processed_T=None):
 110     with _h5py.File(filename, 'a') as f:
 111         cwg = _h5_create_group(f, group)
 112         if raw_T is not None:
 113             try:
 114                 del cwg['raw']
 115             except KeyError:
 116                 pass
 117             cwg['raw'] = raw_T
 118         if temperature_config is not None:
 119             config_cwg = _h5_create_group(cwg, 'config')
 120             storage = _HDF5_Storage()
 121             storage.save(config=temperature_config, group=config_cwg)
 122         if processed_T is not None:
 123             try:
 124                 del cwg['processed']
 125             except KeyError:
 126                 pass
 127             cwg['processed'] = processed_T
 128
 129 def T_load(filename, group='/'):
 130     assert group.endswith('/')
 131     raw_T = processed_T = None
 132     with _h5py.File(filename, 'a') as f:
 133         try:
 134             raw_T = f[group+'raw'][...]
 135         except KeyError:
 136             pass
 137         temperature_config = _TemperatureConfig(storage=_HDF5_Storage(
 138                 filename=filename, group=group+'config/'))
 139         try:
 140             processed_T = f[group+'processed'][...]
 141         except KeyError:
 142             pass
 143     temperature_config.load()
 144     return (raw_T, temperature_config, processed_T)
 145
 146 def T_plot(raw_T=None, processed_T=None):
 147     if not _matplotlib:
 148         raise _matplotlib_import_error
 149     figure = _matplotlib_pyplot.figure()
 150     timestamp = _time.strftime('%H%M%S')
 151     if raw_T is None:
 152         if processed_T is None:
 153             return  # nothing to plot
 154         axes1 = None
 155         axes2 = figure.add_subplot(1, 1, 1)
 156     elif processed_T is None:
 157         axes1 = figure.add_subplot(1, 1, 1)
 158         axes2 = None
 159     else:
 160         axes1 = figure.add_subplot(2, 1, 1)
 161         axes2 = figure.add_subplot(2, 1, 2)
 162     if axes1:
 163         axes1.set_title('Raw Temperatures %s' % timestamp)
 164         axes1.plot(raw_T, label='raw')
 165     if axes2:
 166         axes2.set_title('Processed Temperatures %s' % timestamp)
 167         axes2.plot(processed_T, label='processed')
 168     if hasattr(figure, 'show'):
 169         figure.show()