Only calculate relative errors if a previous value exists.
[calibcant.git] / calibcant / vib.py
1 # calibcant - tools for thermally calibrating AFM cantilevers
2 #
3 # Copyright (C) 2011-2012 W. Trevor King <wking@drexel.edu>
4 #
5 # This file is part of calibcant.
6 #
7 # calibcant is free software: you can redistribute it and/or modify it under
8 # the terms of the GNU General Public License as published by the Free Software
9 # Foundation, either version 3 of the License, or (at your option) any later
10 # version.
11 #
12 # calibcant is distributed in the hope that it will be useful, but WITHOUT ANY
13 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
14 # A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
15 #
16 # You should have received a copy of the GNU General Public License along with
17 # calibcant.  If not, see <http://www.gnu.org/licenses/>.
18
19 """Acquire, save, and load cantilever thermal vibration data.
20
21 For measuring the cantilever's spring constant.
22
23 The relevent physical quantity is:
24   Vphoto   The photodiode vertical deflection voltage (what we measure)
25 """
26
27 from numpy import zeros as _zeros
28 from FFT_tools import ceil_pow_of_two as _ceil_pow_of_two
29 from pycomedi.constant import TRIG_SRC as _TRIG_SRC
30 from pycomedi.utility import inttrig_insn as _inttrig_insn
31 from pycomedi.utility import Reader as _Reader
32
33 from . import LOG as _LOG
34 from .vib_analyze import vib_analyze as _vib_analyze
35 from .vib_analyze import vib_save as _vib_save
36
37
38 def vib_acquire(piezo, vibration_config):
39     """Record thermal vibration data for `piezo` at its current position.
40
41     Inputs:
42       piezo             a pypiezo.afm.AFMPiezo instance
43       vibration_config  a .config._VibrationConfig instance
44     """
45     _LOG.debug('prepare vibration aquisition command')
46
47     # round up to the nearest power of two, for efficient FFT-ing
48     n_samps = _ceil_pow_of_two(
49         vibration_config['sample-time']*vibration_config['frequency'])
50     time = n_samps / vibration_config['frequency']
51     scan_period_ns = int(1e9 / vibration_config['frequency'])
52
53     input_channel = piezo.input_channel_by_name('deflection')
54     channel = input_channel.channel
55
56     channels = [channel]
57
58     dtype = piezo.deflection_dtype()
59     data = _zeros((n_samps, len(channels)), dtype=dtype)
60
61     cmd = channel.subdevice.get_cmd_generic_timed(
62         len(channels), scan_period_ns)
63     cmd.start_src = _TRIG_SRC.int
64     cmd.start_arg = 0
65     cmd.stop_src = _TRIG_SRC.count
66     cmd.stop_arg = n_samps
67     cmd.chanlist = channels
68     channel.subdevice.cmd = cmd
69     for i in range(3):
70         rc = channel.subdevice.command_test()
71         if rc == None: break
72         _LOG.debug('command test %d: %s' % (i,rc))
73
74     _LOG.info('get %g seconds of vibration data at %g Hz'
75               % (vibration_config['sample-time'],
76                  vibration_config['frequency']))
77     channel.subdevice.command()
78     reader = _Reader(channel.subdevice, data)
79     reader.start()
80     channel.subdevice.device.do_insn(_inttrig_insn(channel.subdevice))
81     reader.join()
82     data = data.reshape((data.size,))
83     return data
84
85 def vib(piezo, vibration_config, filename, group='/'):
86     """Wrapper around vib_acquire(), vib_analyze(), vib_save().
87
88     >>> import os
89     >>> import tempfile
90     >>> from h5config.storage.hdf5 import pprint_HDF5
91     >>> from pycomedi.device import Device
92     >>> from pycomedi.subdevice import StreamingSubdevice
93     >>> from pycomedi.channel import AnalogChannel
94     >>> from pycomedi.constant import AREF, SUBDEVICE_TYPE, UNIT
95     >>> from pypiezo.afm import AFMPiezo
96     >>> from pypiezo.base import InputChannel
97     >>> from pypiezo.config import ChannelConfig
98     >>> from .config import VibrationConfig
99
100     Setup an `AFMPiezo` instance.
101
102     >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
103     >>> os.close(fd)
104
105     >>> d = Device('/dev/comedi0')
106     >>> d.open()
107
108     >>> s_in = d.find_subdevice_by_type(SUBDEVICE_TYPE.ai,
109     ...     factory=StreamingSubdevice)
110
111     >>> channel = s_in.channel(0, factory=AnalogChannel, aref=AREF.diff)
112     >>> channel.range = channel.find_range(
113     ...     unit=UNIT.volt, min=-10, max=10)
114     >>> channel_config = ChannelConfig()
115     >>> channel_config['name'] = 'deflection'
116
117     >>> c = InputChannel(config=channel_config, channel=channel)
118     >>> c.setup_config()
119
120     >>> piezo = AFMPiezo(axes=[], inputs=[c])
121
122     Test a vibration:
123
124     >>> vibration_config = VibrationConfig()
125     >>> vib(piezo, vibration_config, filename, group='/vibration')
126     TODO: replace skipped example data with real-world values
127     >>> pprint_HDF5(filename)  # doctest: +ELLIPSIS, +REPORT_UDIFF
128     /
129       /vibration
130         /vibration/config
131           /vibration/config/deflection
132             <HDF5 dataset "channel": shape (), type "<i4">
133               0
134             <HDF5 dataset "conversion-coefficients": shape (2,), type "<f8">
135               [ -1.00000000e+01   3.05180438e-04]
136             <HDF5 dataset "conversion-origin": shape (), type "<f8">
137               0.0
138             <HDF5 dataset "device": shape (), type "|S12">
139               /dev/comedi0
140             <HDF5 dataset "inverse-conversion-coefficients": shape (2,), type "<f8">
141               [    0.    3276.75]
142             <HDF5 dataset "inverse-conversion-origin": shape (), type "<f8">
143               -10.0
144             <HDF5 dataset "maxdata": shape (), type "<i8">
145               65535
146             <HDF5 dataset "name": shape (), type "|S10">
147               deflection
148             <HDF5 dataset "range": shape (), type "<i4">
149               0
150             <HDF5 dataset "subdevice": shape (), type "<i4">
151               0
152           /vibration/config/vibration
153             <HDF5 dataset "chunk-size": shape (), type "<i4">
154               2048
155             <HDF5 dataset "frequency": shape (), type "<f8">
156               50000.0
157             <HDF5 dataset "maximum-fit-frequency": shape (), type "<f8">
158               25000.0
159             <HDF5 dataset "minimum-fit-frequency": shape (), type "<f8">
160               500.0
161             <HDF5 dataset "model": shape (), type "|S12">
162               Breit-Wigner
163             <HDF5 dataset "overlap": shape (), type "|b1">
164               False
165             <HDF5 dataset "sample-time": shape (), type "<i4">
166               1
167             <HDF5 dataset "window": shape (), type "|S4">
168               Hann
169         <HDF5 dataset "processed": shape (), type "<f8">
170           ...
171         /vibration/raw
172           <HDF5 dataset "deflection": shape (65536,), type "<u2">
173             [...]
174
175     Close the Comedi device.
176
177     >>> d.close()
178
179     Cleanup our temporary config file.
180
181     >>> os.remove(filename)
182     """
183     deflection_input_channel = piezo.input_channel_by_name('deflection')
184     
185     deflection_channel_config = deflection_input_channel.config
186
187     deflection = vib_acquire(piezo, vibration_config)
188     variance = _vib_analyze(
189         deflection, vibration_config, deflection_channel_config)
190     _vib_save(
191         filename, group, raw_vibration=deflection,
192         vibration_config=vibration_config,
193         deflection_channel_config=deflection_channel_config,
194         processed_vibration=variance)
195     return variance