3 # calibcant - tools for thermally calibrating AFM cantilevers
5 # Copyright (C) 2007,2008, William Trevor King
7 # This program is free software; you can redistribute it and/or
8 # modify it under the terms of the GNU General Public License as
9 # published by the Free Software Foundation; either version 3 of the
10 # License, or (at your option) any later version.
12 # This program is distributed in the hope that it will be useful, but
13 # WITHOUT ANY WARRANTY; without even the implied warranty of
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
15 # See the GNU General Public License for more details.
17 # You should have received a copy of the GNU General Public License
18 # along with this program; if not, write to the Free Software
19 # Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
22 # The author may be contacted at <wking@drexel.edu> on the Internet, or
23 # write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
24 # Philadelphia PA 19104, USA.
27 Separate the more general calib_analyze() from the other calib_*()
28 functions in calibcant. Also provide a command line interface
29 for analyzing data acquired through other workflows.
31 The relevent physical quantities are :
32 Vzp_out Output z-piezo voltage (what we generate)
33 Vzp Applied z-piezo voltage (after external ZPGAIN)
34 Zp The z-piezo position
35 Zcant The cantilever vertical deflection
36 Vphoto The photodiode vertical deflection voltage (what we measure)
37 Fcant The force on the cantilever
38 T The temperature of the cantilever and surrounding solution
39 (another thing we measure)
40 k_b Boltzmann's constant
42 Which are related by the parameters :
44 zpSensitivity Zp / Vzp
45 photoSensitivity Vphoto / Zcant
50 from splittable_kwargs import splittableKwargsFunction, \
51 make_splittable_kwargs_function
54 import common # common module for the calibcant package
55 import config # config module for the calibcant package
56 import T_analyze # T_analyze module for the calibcant package
58 kb = 1.3806504e-23 # Boltzmann's constant in J/K
60 #@splittableKwargsFunction((calib_plot, 'bumps', 'Ts', 'vibs'))
61 # Some of the child functions aren't yet defined, so postpone
62 # make-splittable until later in the module.
63 def calib_analyze(bumps, Ts, vibs, **kwargs) :
65 Analyze data from get_calibration_data()
66 return (k, k_s, !!!a2_r, T_r, one_o_Vp2_r)
67 Inputs (all are arrays of recorded data) :
68 bumps measured (V_photodiode / nm_tip) proportionality constant
69 Ts measured temperature (K)
70 vibs measured V_photodiode variance in free solution (V**2)
72 k cantilever spring constant (in N/m, or equivalently nN/nm)
73 k_s standard deviation in our estimate of k
74 !!!a2_r relative error in a**2
75 !!!T_r relative error in T
76 !!!one_o_Vp2_r relative error in 1/Vphotodiode_variance
78 We're assuming vib is mostly from thermal cantilever vibrations
79 (and then only from vibrations in the single vertical degree of freedom),
80 and not from other noise sources.
81 The various relative errors are returned to help you gauge the
82 largest source of random error in your measurement of k.
83 If one of them is small, don't bother repeating that measurment too often.
84 If one is large, try repeating that measurement more.
85 Remember that you need enough samples to have a valid error estimate in
86 the first place, and that none of this addresses any systematic errors.
88 calib_plot_kwargs, = calib_analyze._splitargs(calib_analyze, kwargs)
89 photoSensitivity2 = bumps**2
90 one_o_Vphoto2 = 1/vibs
92 photoSensitivity2_m = photoSensitivity2.mean()
94 one_o_Vphoto2_m = one_o_Vphoto2.mean()
95 # Vphoto / photoSensitivity = x
96 # k = kb T / <x**2> = kb T photoSensitiviy**2 * (1e9nm/m)**2 /
99 # units, photoSensitivity = Vphoto/(Zcant in nm),
100 # so Vphoto/photoSensitivity = Zcant in nm
101 # so k = J/K * K / nm^2 * (1e9nm/m)**2 = N/m
102 k = kb * T_m * photoSensitivity2_m * one_o_Vphoto2_m * 1e18
104 # propogation of errors !!!
105 # first, get standard deviations
106 photoSensitivity2_s = photoSensitivity2.std()
108 one_o_Vphoto2_s = one_o_Vphoto2.std()
109 # !!!!now, get relative errors
110 photoSensitivity2_r = photoSensitivity2_s / photoSensitivity2_m
112 one_o_Vphoto2_r = one_o_Vphoto2_s / one_o_Vphoto2_m
114 k_s = k*(photoSensitivity2_r**2 + T_r**2 + one_o_Vphoto2_r**2)**0.5
116 calib_plot(bumps, Ts, vibs, **calib_plot_kwargs)
119 photoSensitivity2_m, photoSensitivity2_s,
120 T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s)
122 @splittableKwargsFunction()
123 def string_errors(k_m, k_s,
124 photoSensitivity2_m, photoSensitivity2_s,
126 one_o_Vphoto2_m, one_o_Vphoto2_s) :
128 photoSensitivity2_r = photoSensitivity2_s / photoSensitivity2_m
130 one_o_Vphoto2_r = one_o_Vphoto2_s / one_o_Vphoto2_m
131 string = "Variable (units) : mean +/- std. dev. (relative error)\n"
132 string += "Cantilever k (N/m) : %g +/- %g (%g)\n" \
134 string += "photoSensitivity**2 (V/nm)**2 : %g +/- %g (%g)\n" \
135 % (photoSensitivity2_m, photoSensitivity2_s, photoSensitivity2_r)
136 string += "T (K) : %g +/- %g (%g)\n" \
138 string += "1/Vphoto**2 (1/V)**2 : %g +/- %g (%g)" \
139 % (one_o_Vphoto2_m, one_o_Vphoto2_s, one_o_Vphoto2_r)
142 @splittableKwargsFunction()
143 def calib_save(bumps, Ts, vibs, log_dir=None) :
145 Save a dictonary with the bump, T, and vib data.
148 data = {'bump':bumps, 'T':Ts, 'vib':vibs}
149 log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
151 log.write_dict_of_arrays(data)
153 def calib_load(datafile) :
154 "Load the dictionary data, using data_logger.date_load()"
155 dl = data_logger.data_load()
156 data = dl.read_dict_of_arrays(path)
157 return (data['bump'], data['T'], data['vib'])
159 def calib_save_analysis(k, k_s,
160 photoSensitivity2_m, photoSensitivity2_s,
161 T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s,
164 log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
165 log_name="calib_analysis_text")
167 "k (N/m) : %g +/- %g\n" % (k, k_s)
168 + "photoSensitivity**2 (V/nm)**2 : %g +/- %g\n" % \
169 (photoSensitivity2_m, photoSensitivity2_s)
170 + "T (K) : %g +/- %g\n" % (T_m, T_s)
171 + "1/Vp**2 (1/V)**2 : %g +/- %g\n" % \
172 (one_o_Vphoto2_m, one_o_Vphoto2_s)
175 @splittableKwargsFunction()
176 def calib_plot(bumps, Ts, vibs, plotVerbose=False) :
177 if plotVerbose or config.PYLAB_VERBOSE :
178 common._import_pylab()
179 common._pylab.figure(config.BASE_FIGNUM+4)
180 common._pylab.subplot(311)
181 common._pylab.plot(bumps, 'g.')
182 common._pylab.title('Photodiode sensitivity (V/nm)')
183 common._pylab.subplot(312)
184 common._pylab.plot(Ts, 'r.')
185 common._pylab.title('Temperature (K)')
186 common._pylab.subplot(313)
187 common._pylab.plot(vibs, 'b.')
188 common._pylab.title('Thermal deflection variance (Volts**2)')
191 make_splittable_kwargs_function(calib_analyze,
192 (calib_plot, 'bumps', 'Ts', 'vibs'))
194 @splittableKwargsFunction((calib_analyze, 'bumps', 'Ts', 'vibs'))
195 def calib_load_analyze_tweaks(bump_tweaks, vib_tweaks, T_tweaks=None) :
196 raise NotImplementedError
197 a = read_tweaked_bumps(bump_tweaks)
198 vib = V_photo_variance_from_file(vib_tweaks)
201 return analyze_calibration_data(a, T, vib, log_dir=log_dir)
203 # commandline interface functions
206 def array_from_string(string):
208 for num in string.split(',') :
209 ret.append(float(num))
211 return numpy.array(ret)
213 def read_data(ifile):
214 "ifile can be a filename string or open (seekable) file object"
215 unlabeled_data=scipy.io.read_array(file)
216 return unlabeled_data
218 def get_array(string, filename, name) :
219 "get an array from supplied command line options"
221 array = array_from_string(string)
222 elif filename != None :
223 array = read_data(filename)
225 raise Exception, "no %s information given" % (name)
228 if __name__ == '__main__' :
229 # command line interface
230 from optparse import OptionParser
232 usage_string = ('%prog <bumps> <temps> <vibs>\n'
233 '2008, W. Trevor King.\n'
235 'Takes arrays of Vphotodiode sensitivity (V/nm), Temperature (K), \n'
236 'and Vibration variance (V**2) as comma seperated lists.\n'
237 'Returns the cantilever spring constant (pN/nm).\n'
239 ' $ %prog -b 0.02,0.03,0.025 -t 298.2,300.1 -v 6e-9,5.5e-9\n'
241 parser = OptionParser(usage=usage_string, version='%prog '+common.VERSION)
242 parser.add_option('-b', '--bump-string', dest='bump_string',
243 help='comma seperated photodiode sensitivities (V/nm)',
244 type='string', metavar='BUMPS')
245 parser.add_option('-t', '--temp-string', dest='temp_string',
246 help='comma seperated temperatures (K)',
247 type='string', metavar='TEMPS')
248 parser.add_option('-v', '--vib-string', dest='vib_string',
249 help='comma seperated vibration variances (V**2)',
250 type='string', metavar='VIBS')
251 parser.add_option('-B', '--bump-file', dest='bump_file',
252 help='comma seperated photodiode sensitivities (V/nm)',
253 type='string', metavar='BUMPFILE')
254 parser.add_option('-T', '--temp-file', dest='temp_file',
255 help='comma seperated temperatures (K)',
256 type='string', metavar='TEMPFILE')
257 parser.add_option('-V', '--vib-file', dest='vib_file',
258 help='comma seperated vibration variances (V**2)',
259 type='string', metavar='VIBFILE')
260 parser.add_option('-C', '--celsius', dest='celsius',
261 help='Use Celsius input temperatures instead of Kelvin (defaul %default)\n',
262 action='store_true', default=False)
263 parser.add_option('-o', '--output-file', dest='ofilename',
264 help='write output to FILE (default stdout)',
265 type='string', metavar='FILE')
266 parser.add_option('-p', '--plot-inputs', dest='plot',
267 help='plot the input arrays to check their distribution',
268 action='store_true', default=False)
269 parser.add_option('', '--verbose', dest='verbose', action='store_true',
270 help='print lots of debugging information',
273 options,args = parser.parse_args()
276 bumps = get_array(options.bump_string, options.bump_file, "bump")
277 temps = get_array(options.temp_string, options.temp_file, "temp")
278 vibs = get_array(options.vib_string, options.vib_file, "vib")
281 calib_plot(bumps, temps, vibs)
284 for i in range(len(temps)) :
285 temps[i] = T_analyze.C_to_K(temps[i])
287 km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s = calib_analyze(bumps, temps, vibs)
290 print >> sys.stderr, \
291 string_errors(km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s)
293 if options.ofilename != None :
294 print >> file(options.ofilename, 'w'), km