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 bump_analyze() from the other bump_*()
28 functions in calibcant. Also provide a command line interface
29 for analyzing data acquired through other workflows.
31 The relevant 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)
38 Which are related by the parameters :
40 zpSensitivity Zp / Vzp
41 photoSensitivity Vphoto / Zcant
43 photoSensitivity is measured by bumping the cantilever against the
44 surface, where Zp = Zcant (see calibrate.bump_aquire()). The measured
45 slope Vphoto/Vout is converted to photoSensitivity with bump_analyze().
50 import common # common module for the calibcant package
51 import config # config module for the calibcant package
53 from splittable_kwargs import splittableKwargsFunction, \
54 make_splittable_kwargs_function
57 #@splittableKwargsFunction((bump_fit, 'zpiezo_output_bits',
58 # 'deflection_input_bits'))
59 # Some of the child functions aren't yet defined, so postpone
60 # make-splittable until later in the module.
61 def bump_analyze(data, zpGain=config.zpGain,
62 zpSensitivity=config.zpSensitivity,
63 Vzp_out2V=config.Vzp_out2V,
64 Vphoto_in2V=config.Vphoto_in2V,
67 Return the slope of the bump ;).
69 data dictionary of data in DAC/ADC bits
70 Vzp_out2V function that converts output DAC bits to Volts
71 Vphoto_in2V function that converts input ADC bits to Volts
72 zpGain zpiezo applied voltage per output Volt
73 zpSensitivity nm zpiezo response per applied Volt
75 photoSensitivity (Vphoto/Zcant) in Volts/nm
76 Checks for strong correlation (r-value) and low randomness chance (p-value)
78 With the current implementation, the data is regressed in DAC/ADC bits
79 and THEN converted, so we're assuming that both conversions are LINEAR.
80 if they aren't, rewrite to convert before the regression.
82 bump_fit_kwargs, = bump_analyze._splitargs(bump_analyze, kwargs)
83 Vphoto2Vzp_out_bit = bump_fit(data['Z piezo output'],
84 data['Deflection input'],
86 scale_Vzp_bits2V = Vzp_out2V(1) - Vzp_out2V(0)
87 scale_Vphoto_bits2V = Vphoto_in2V(1) - Vphoto_in2V(0)
88 Vphoto2Vzp_out = Vphoto2Vzp_out_bit *scale_Vphoto_bits2V / scale_Vzp_bits2V
89 # 1 / (Vzp/Vzp_out * Zp/Vzp * Zcant/Zp )
90 Vzp_out2Zcant = 1.0/ (zpGain * zpSensitivity) # * 1
91 return Vphoto2Vzp_out * Vzp_out2Zcant
93 def limited_quadratic(x, params):
96 flat region (off-surface)
97 quadratic region (in-contact)
98 flat region (high-voltage-rail)
100 x_contact (x value for the surface-contact kink)
101 y_contact (y value for the surface-contact kink)
102 slope (dy/dx at the surface-contact kink)
103 quad (d**2 y / dx**2, allow decreasing sensitivity with increased x)
105 high_voltage_rail = 2**16 - 1 # bits
106 x_contact,y_contact,slope,quad = params
107 y = slope*(x-x_contact) + quad*(x-x_contact)**2+ y_contact
108 y = numpy.clip(y, y_contact, high_voltage_rail)
111 def limited_quadratic_param_guess(x, y) :
113 Guess rough parameters for a limited_quadratic model. Assumes the
114 bump approaches (raising the deflection as it does so) first.
115 Retracting after the approach is optional. Approximates the contact
116 position and an on-surface (high) position by finding first crossings
117 of thresholds 0.3 and 0.7 of the y value's total range. Not the
118 most efficient algorithm, but it seems fairly robust.
120 y_contact = float(y.min())
121 y_max = float(y.max())
123 y_low = y_contact + 0.3 * (y_max-y_contact)
124 y_high = y_contact + 0.7 * (y_max-y_contact)
128 while y[i] < y_high :
131 x_contact = float(x[i_low])
132 x_high = float(x[i_high])
133 slope = (y_high - y_contact) / (x_high - x_contact)
135 return (x_contact, y_contact, slope, quad)
137 def limited_quadratic_sensitivity(params):
139 Return the estimated sensitivity to small deflections according to
140 limited_quadratic fit parameters.
145 @splittableKwargsFunction()
146 def bump_fit(zpiezo_output_bits, deflection_input_bits,
147 paramGuesser=limited_quadratic_param_guess,
148 model=limited_quadratic,
149 sensitivity_from_fit_params=limited_quadratic_sensitivity,
151 x = zpiezo_output_bits
152 y = deflection_input_bits
153 def residual(p, y, x) :
154 return model(x, p) - y
155 paramGuess = paramGuesser(x, y)
156 p,cov,info,mesg,ier = \
157 scipy.optimize.leastsq(residual, paramGuess, args=(y, x),
158 full_output=True, maxfev=int(10e3))
159 if config.TEXT_VERBOSE :
160 print "Fitted params:",p
161 print "Covariance mx:",cov
165 print "Solution converged"
167 print "Solution did not converge"
168 if plotVerbose or config.PYLAB_VERBOSE :
169 yguess = model(x, paramGuess)
170 #yguess = None # Don't print the guess, since I'm convinced it's ok ;).
172 bump_plot(data={"Z piezo output":x, "Deflection input":y},
173 yguess=yguess, yfit=yfit, plotVerbose=plotVerbose)
174 return sensitivity_from_fit_params(p)
176 @splittableKwargsFunction()
177 def bump_save(data, log_dir=None) :
178 "Save the dictionary data, using data_logger.data_log()"
180 log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
182 log.write_dict_of_arrays(data)
184 def bump_load(datafile) :
185 "Load the dictionary data, using data_logger.date_load()"
186 dl = data_logger.data_load()
187 data = dl.read_dict_of_arrays(datafile)
190 @splittableKwargsFunction()
191 def bump_plot(data, yguess=None, yfit=None, plotVerbose=False) :
192 "Plot the bump (Vphoto vs Vzp) if plotVerbose or PYLAB_VERBOSE == True"
193 if plotVerbose or config.PYLAB_VERBOSE :
194 common._import_pylab()
195 common._pylab.figure(config.BASE_FIGNUM)
196 if yfit != None: # two subplot figure
197 common._pylab.subplot(211)
198 common._pylab.hold(False)
199 common._pylab.plot(data["Z piezo output"], data["Deflection input"],
201 common._pylab.hold(True)
203 common._pylab.plot(data["Z piezo output"], yguess,
206 common._pylab.plot(data["Z piezo output"], yfit,
208 common._pylab.title("bump surface")
209 common._pylab.legend(loc='upper left')
210 common._pylab.xlabel("Z piezo output voltage (bits)")
211 common._pylab.ylabel("Photodiode input voltage (bits)")
213 # second subplot for residual
214 common._pylab.subplot(212)
215 common._pylab.plot(data["Z piezo output"],
216 data["Deflection input"] - yfit,
217 'r-', label='residual')
218 common._pylab.legend(loc='upper right')
219 common._pylab.xlabel("Z piezo output voltage (bits)")
220 common._pylab.ylabel("Photodiode input voltage (bits)")
223 make_splittable_kwargs_function(bump_analyze,
224 (bump_fit, 'zpiezo_output_bits',
225 'deflection_input_bits'))
227 @splittableKwargsFunction((bump_analyze, 'data'))
228 def bump_load_analyze_tweaked(tweak_file, **kwargs):
229 "Load the output file of tweak_calib_bump.sh, return an array of slopes"
230 bump_analyze_kwargs, = \
231 bump_load_analyze_tweaked._splitargs(bump_load_analyze_tweaked, kwargs)
232 photoSensitivity = []
233 for line in file(tweak_file, 'r') :
234 parsed = line.split()
235 path = parsed[0].strip()
236 if path[0] == '#' : # a comment
238 if config.TEXT_VERBOSE :
239 print "Reading data from %s with ranges %s" % (path, parsed[1:])
241 full_data = bump_load(path)
242 if len(parsed) == 1 :
243 data = full_data # use whole bump
245 # use the listed sections
248 for rng in parsed[1:] :
252 zp.extend(full_data['Z piezo output'][starti:stopi])
253 df.extend(full_data['Deflection input'][starti:stopi])
254 data = {'Z piezo output': numpy.array(zp),
255 'Deflection input': numpy.array(df)}
256 pSi = bump_analyze(data, **bump_analyze_kwargs)
257 photoSensitivity.append(pSi)
258 return numpy.array(photoSensitivity, dtype=numpy.float)
260 # commandline interface functions
263 def read_data(ifile):
264 "ifile can be a filename string or open (seekable) file object"
265 if ifile == None : ifile = sys.stdin
266 unlabeled_data=scipy.io.read_array(ifile)
268 data['Z piezo output'] = unlabeled_data[:,0]
269 data['Deflection input'] = unlabeled_data[:,1]
272 def remove_further_than(data, zp_crit) :
274 ndata['Z piezo output'] = []
275 ndata['Deflection input'] = []
276 for zp,df in zip(data['Z piezo output'],data['Deflection input']) :
278 ndata['Z piezo output'].append(zp)
279 ndata['Deflection input'].append(df)
282 if __name__ == '__main__' :
283 # command line interface
284 from optparse import OptionParser
286 usage_string = ('%prog <input-file>\n'
287 '2008, W. Trevor King.\n'
289 'There are two operation modes, one to analyze a single bump file,\n'
290 'and one to analyze tweak files.\n'
292 'Single file mode (the default) :\n'
293 'Scales raw DAC/ADC bit data and fits a straight line.\n'
294 'Returns photodiode sensitivity Vphotodiode/Zcantilever in V/nm.\n'
295 '<input-file> should be whitespace-delimited, 2 column ASCII\n'
296 'without a header line. e.g: "<zp_DAC>\\t<deflection_ADC>\\n"\n'
299 'Runs the same analysis as in single file mode for each bump in\n'
300 'a tweak file. Each line in the tweak file specifies a single bump.\n'
301 'Blank lines and those beginning with a pound sign (#) are ignored.\n'
302 'The format of a line is a series of whitespace-separated fields--\n'
303 'a base file path followed by optional point index ranges, e.g.:\n'
304 '20080919/20080919132500_bump_surface 10:651 1413:2047\n'
305 'which only discards all points outside the index ranges [10,651)\n'
306 'and [1413,2047) (indexing starts at 0).\n'
308 parser = OptionParser(usage=usage_string, version='%prog '+common.VERSION)
309 parser.add_option('-o', '--output-file', dest='ofilename',
310 help='write output to FILE (default stdout)',
311 type='string', metavar='FILE')
312 parser.add_option('-c', '--comma-out', dest='comma_out', action='store_true',
313 help='Output comma-seperated values (default %default)',
315 parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
316 help='Produce pylab fit checks during execution',
318 parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
319 help='Run in tweak-file mode',
321 parser.add_option('-d', '--datalogger-mode', dest='datalogger_mode', action='store_true',
322 help='Run input files with datalogger.read_dict_of_arrays(). This is useful, for example, to test a single line from a tweakfile.',
324 parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
325 help='Print lots of debugging information',
328 options,args = parser.parse_args()
330 assert len(args) >= 1, "Need an input file"
334 if options.ofilename != None :
335 ofile = file(options.ofilename, 'w')
338 config.TEXT_VERBOSE = options.verbose
339 config.PYLAB_INTERACTIVE = False
340 config.PYLAB_VERBOSE = options.pylab
341 config.GNUPLOT_VERBOSE = False
343 if options.tweakmode == False :
344 if options.datalogger_mode:
345 data = bump_load(ifilename)
347 data = read_data(ifilename)
348 photoSensitivity = bump_analyze(data)
350 print >> ofile, photoSensitivity
351 else : # tweak file mode
352 slopes = bump_load_analyze_tweaked(ifilename)
353 if options.comma_out :
357 common.write_array(ofile, slopes, sep)
359 if options.ofilename != None :