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