--- /dev/null
+.PHONY : all check dist clean
+
+all : dummy_py
+
+dummy_py : setup.py FFT_tools.py
+ python setup.py install --home=~
+ echo "dummy for Makefile dependencies" > $@
+
+check : all
+ python FFT_tools.py
+
+dist :
+ python setup.py sdist
+ scp dist/data_logger*tar.gz einstein:public_html/code/python/
+
+clean :
+ python setup.py clean
+ rm -rf build dist data_logger.egg-info
+ rm -f dummy_py *.pyc
--- /dev/null
+== Dependencies ==
+
+linfit (depends in turn on scipy.stats.linregress)
+
+*** Command line interface ***
+
+** bumps **
+
+Create a two-column ascii data file for fitting bumps with (for example)
+ $ int16s_to_ascii_array ~/rsrch/data/z_piezo_calib/20080116/20080116090005_bump_surface_Z_piezo_output \
+ ~/rsrch/data/z_piezo_calib/20080116/20080116090005_bump_surface_Deflection_input > d
+
+Then fit slope with
+ $ python calibcant_bump_analyze.py -c d > stdout 2> stderr
+ $ cat stdout
+ 0.00958650972452
+ $ cat stderr
+ not cutting
+
+** vibs **
+
+Create a two-column ascii data file for fitting bumps with (for example)
+ $ cd ~/src/comedi_simult_AIO
+ $ cmd -c0 -n1 -F200000 -N2097152 > vib.raw # 2097152 = 2^21 points ~= 10 seconds of data
+ $ unitary_avg_fft_from_raw.py vib.raw
+Creates vib.raw.psd, with frequency in Hz in the first column, and amplitude in Volts^2/Hz in the second column.
+Then get the variance with
+ $ cd ~/src/calibrate_cantilever
+ $ python psd_filter_analyze.py -s9 -tmedian -g ~/src/comedi_simult_AIO/vib.raw 2> check.gp
+ <Lorentzian-fit-variance>
+You can check the validity of the fit with
+ $ gnuplot check.gp -
+And adjust the min and max of the fitted frequency range as neccessary
+ $ python psd_filter_analyze.py -s9 -tmedian -m1000 -M8500 -g ~/src/comedi_simult_AIO/vib.raw 2> check.gp
+ <new Lorentzian-fit-variance>
+Until the fit looks reasonable.
+
+
+** calibrating **
+
+Once you've assembled a few slopes, vibs, and Ts, put it all together to calculate k with
+ $ python calibcant_analyze.py -vC 0.02,0.03,0.025 22.5,22.1 6e-9,5.5e-9 > stdout 2> stderr
+ $ cat stdout
+ 456.069575483
+ $ cat stderr
+ Variable (units) : mean +/- std. dev. (relative error)
+ Cantilever k (pN/nm) : 456.07 +/- 146.671 (0.321599)
+ photoSensitivity**2 (V/nm)**2 : 0.000641667 +/- 0.000204464 (0.318645)
+ T (K) : 295.45 +/- 0.2 (0.000676933)
+ 1/Vphoto**2 (1/V)**2 : 1.74242e+08 +/- 7.57576e+06 (0.0434783)
--- /dev/null
+#!/usr/bin/python
+#
+# calibcant - tools for thermally calibrating AFM cantilevers
+#
+# Copyright (C) 2007,2008, William Trevor King
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as
+# published by the Free Software Foundation; either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+# See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
+# 02111-1307, USA.
+#
+# The author may be contacted at <wking@drexel.edu> on the Internet, or
+# write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
+# Philadelphia PA 19104, USA.
+
+"""
+Separate the more general T_analyze() from the other T_*()
+functions in calibcant. Also provide a command line interface
+for analyzing data acquired through other workflows.
+
+The relevant physical quantities are :
+ T Temperature at which thermal vibration measurements were aquired
+"""
+
+import numpy
+import common # common module for the calibcant package
+import config # config module for the calibcant package
+import data_logger
+import z_piezo_utils
+import linfit
+
+def C_to_K(celsius) :
+ "Convert Celsius -> Kelvin."
+ return celsius + 273.15
+
+def T_analyze(T, convert_to_K=C_to_K) :
+ """
+ Not much to do here, just convert to Kelvin.
+ Uses C_to_K (defined above) by default.
+ """
+ try : # if T is an array, convert element by element
+ for i in range(len(T)) :
+ T[i] = convert_to_K(T[i])
+ except TypeError : # otherwise, make an array from a single T
+ T = [convert_to_K(T)]
+ return T
+
+def T_save(T, log_dir=None) :
+ """
+ Save either a single T (if you are crazy :p),
+ or an array of them (more reasonable)
+ """
+ T = numpy.array(T, dtype=numpy.float)
+ if log_dir != None :
+ log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
+ log_name="T_float")
+ log.write_binary(T.tostring())
+ if LOG_DATA != None :
+ log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
+ log_name="T_float")
+ log.write_binary(T.tostring())
+
+def T_load(datafile=None) :
+ """
+ Load the saved T array (possibly of length 1), and return it. If
+ datafile == None, return an array of one config.DEFAULT_TEMP
+ instead.
+ """
+ if datafile == None :
+ return numpy.array([config.DEFAULT_TEMP], dtype=numpy.float)
+ else :
+ dl = data_logger.data_load()
+ return dl.read_binary(datafile)
+
+def T_plot(T, plotVerbose) :
+ """
+ Umm, just print the temperature?
+ """
+ if plotVerbose or PYLAB_VERBOSE or TEXT_VERBOSE :
+ print "Temperature ", T
+
+def T_load_analyze_tweaked(tweak_file, convert_to_K=C_to_K, textVerboseFile=None) :
+ "Load all the T array files from a tweak file and return a single array"
+ Ts = []
+ for line in file(tweak_file, 'r') :
+ parsed = line.split()
+ path = parsed[0].split('\n')[0]
+ if textVerboseFile != None :
+ print >> textVerboseFile, "Reading data from %s" % (path)
+ # read the data
+ data = T_load(path)
+ Ts.extend(data)
+ return numpy.array(Ts, dtype=numpy.float)
+
+# commandline interface functions
+import scipy.io, sys
+
+def read_data(ifile):
+ "ifile can be a filename string or open (seekable) file object"
+ if ifile == None : ifile = sys.stdin
+ unlabeled_data=scipy.io.read_array(ifile)
+ return unlabeled_data
+
+if __name__ == '__main__' :
+ # command line interface
+ from optparse import OptionParser
+
+ usage_string = ('%prog <input-file>\n'
+ '2008, W. Trevor King.\n'
+ '\n'
+ 'There are two operation modes, one to analyze a single T (temperature) file,\n'
+ 'and one to analyze tweak files.\n'
+ '\n'
+ 'Single file mode (the default) :\n'
+ 'Reads in single column ASCII file of temperatures and... prints them back out.\n'
+ 'No need to do this, but the option is available for consistency with the other\n'
+ 'calibcant modules.\n'
+ '\n'
+ 'Tweak file mode:\n'
+ 'Runs the same analysis as in single file mode for each T file in\n'
+ 'a tweak file. Each line in the tweak file specifies a single T file.\n'
+ 'A T file contains a sequence of 32 bit floats representing temperature in K.\n'
+ )
+ parser = OptionParser(usage=usage_string, version='%prog 0.1')
+ parser.add_option('-C', '--celsius', dest='celsius',
+ help='Use Celsius input temperatures instead of Kelvin (defaul %default)\n',
+ action='store_true', default=False)
+ parser.add_option('-o', '--output-file', dest='ofilename',
+ help='write output to FILE (default stdout)',
+ type='string', metavar='FILE')
+ parser.add_option('-c', '--comma-out', dest='comma_out', action='store_true',
+ help='Output comma-seperated values (default %default)',
+ default=False)
+ parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
+ help='Run in tweak-file mode',
+ default=False)
+ parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
+ help='Print lots of debugging information',
+ default=False)
+
+ options,args = parser.parse_args()
+ parser.destroy()
+ assert len(args) >= 1, "Need an input file"
+
+ ifilename = args[0]
+
+ if options.ofilename != None :
+ ofile = file(options.ofilename, 'w')
+ else :
+ ofile = sys.stdout
+ if options.verbose == True :
+ vfile = sys.stderr
+ else :
+ vfile = None
+ config.TEXT_VERBOSE = options.verbose
+ config.PYLAB_VERBOSE = False
+ config.GNUPLOT_VERBOSE = False
+ if options.celsius :
+ convert_to_K = C_to_K
+ else :
+ convert_to_K = lambda T : T # no-change function
+
+ if options.tweakmode == False :
+ data = read_data(ifilename)
+ Ts = T_analyze(T, convert_to_K)
+ else : # tweak file mode
+ Ts = T_load_analyze_tweaked(ifilename, convert_to_K, textVerboseFile=vfile)
+
+ if options.comma_out :
+ sep = ','
+ else :
+ sep = '\n'
+ common.write_array(ofile, Ts, sep)
+
+ if options.ofilename != None :
+ ofile.close()
--- /dev/null
+#!/usr/bin/python
+#
+# calibcant - tools for thermally calibrating AFM cantilevers
+#
+# Copyright (C) 2007,2008, William Trevor King
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as
+# published by the Free Software Foundation; either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+# See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
+# 02111-1307, USA.
+#
+# The author may be contacted at <wking@drexel.edu> on the Internet, or
+# write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
+# Philadelphia PA 19104, USA.
+
+"""
+Separate the more general calib_analyze() from the other calib_*()
+functions in calibcant. Also provide a command line interface
+for analyzing data acquired through other workflows.
+
+The relevent physical quantities are :
+ Vzp_out Output z-piezo voltage (what we generate)
+ Vzp Applied z-piezo voltage (after external ZPGAIN)
+ Zp The z-piezo position
+ Zcant The cantilever vertical deflection
+ Vphoto The photodiode vertical deflection voltage (what we measure)
+ Fcant The force on the cantilever
+ T The temperature of the cantilever and surrounding solution
+ (another thing we measure)
+ k_b Boltzmann's constant
+
+Which are related by the parameters :
+ zpGain Vzp_out / Vzp
+ zpSensitivity Zp / Vzp
+ photoSensitivity Vphoto / Zcant
+ k_cant Fcant / Zcant
+"""
+
+import numpy
+import common # common module for the calibcant package
+import config # config module for the calibcant package
+import T_analyze # T_analyze module for the calibcant package
+
+kb = 1.3806504e-23 # Boltzmann's constant in J/K
+
+def calib_analyze(bumps, Ts, vibs) :
+ """
+ Analyze data from get_calibration_data()
+ return (k, k_s, !!!a2_r, T_r, one_o_Vp2_r)
+ Inputs (all are arrays of recorded data) :
+ bumps measured (V_photodiode / nm_tip) proportionality constant
+ Ts measured temperature (K)
+ vibs measured V_photodiode variance in free solution
+ Outputs :
+ k cantilever spring constant (in N/m, or equivalently nN/nm)
+ k_s standard deviation in our estimate of k
+ !!!a2_r relative error in a**2
+ !!!T_r relative error in T
+ !!!one_o_Vp2_r relative error in 1/Vphotodiode_variance
+ Notes :
+ We're assuming vib is mostly from thermal cantilever vibrations
+ (and then only from vibrations in the single vertical degree of freedom),
+ and not from other noise sources.
+ The various relative errors are returned to help you gauge the
+ largest source of random error in your measurement of k.
+ If one of them is small, don't bother repeating that measurment too often.
+ If one is large, try repeating that measurement more.
+ Remember that you need enough samples to have a valid error estimate in
+ the first place, and that none of this addresses any systematic errors.
+ """
+ photoSensitivity2 = bumps**2
+ one_o_Vphoto2 = 1/vibs
+
+ photoSensitivity2_m = photoSensitivity2.mean()
+ T_m = Ts.mean()
+ one_o_Vphoto2_m = one_o_Vphoto2.mean()
+ # Vphoto / photoSensitivity = x
+ # k = kb T / <x**2> = kb T photoSensitiviy**2 * (1e9nm/m)**2 /
+ # <Vphoto_std**2>
+ #
+ # units, photoSensitivity = Vphoto/(Zcant in nm),
+ # so Vphoto/photoSensitivity = Zcant in nm
+ # so k = J/K * K / nm^2 * (1e9nm/m)**2 = N/m
+ k = kb * T_m * photoSensitivity2_m * one_o_Vphoto2_m * 1e18
+
+ # propogation of errors !!!
+ # first, get standard deviations
+ photoSensitivity2_s = photoSensitivity2.std()
+ T_s = Ts.std()
+ one_o_Vphoto2_s = one_o_Vphoto2.std()
+ # !!!!now, get relative errors
+ photoSensitivity2_r = photoSensitivity2_s / photoSensitivity2_m
+ T_r = T_s / T_m
+ one_o_Vphoto2_r = one_o_Vphoto2_s / one_o_Vphoto2_m
+
+ k_s = k*(photoSensitivity2_r**2 + T_r**2 + one_o_Vphoto2_r**2)**0.5
+
+ return (k, k_s,
+ photoSensitivity2_m, photoSensitivity2_s,
+ T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s)
+
+def string_errors(k_m, k_s,
+ photoSensitivity2_m, photoSensitivity2_s,
+ T_m, T_s,
+ one_o_Vphoto2_m, one_o_Vphoto2_s) :
+ k_r = k_s / k_m
+ photoSensitivity2_r = photoSensitivity2_s / photoSensitivity2_m
+ T_r = T_s / T_m
+ one_o_Vphoto2_r = one_o_Vphoto2_s / one_o_Vphoto2_m
+ string = "Variable (units) : mean +/- std. dev. (relative error)\n"
+ string += "Cantilever k (pN/nm) : %g +/- %g (%g)\n" \
+ % (k_m, k_s, k_r)
+ string += "photoSensitivity**2 (V/nm)**2 : %g +/- %g (%g)\n" \
+ % (photoSensitivity2_m, photoSensitivity2_s, photoSensitivity2_r)
+ string += "T (K) : %g +/- %g (%g)\n" \
+ % (T_m, T_s, T_r)
+ string += "1/Vphoto**2 (1/V)**2 : %g +/- %g (%g)" \
+ % (one_o_Vphoto2_m, one_o_Vphoto2_s, one_o_Vphoto2_r)
+ return string
+
+def calib_plot(bumps, Ts, vibs) :
+ from pylab import figure, subplot, plot, title, show
+ figure()
+ subplot(311)
+ plot(bumps, 'g.-')
+ title('Photodiode sensitivity (V/nm)')
+ subplot(312)
+ plot(Ts, 'r.-')
+ title('Temperature (K)')
+ subplot(313)
+ plot(vibs, 'b.-')
+ title('Thermal deflection variance (Volts**2)')
+ show()
+
+# commandline interface functions
+import scipy.io, sys
+
+def array_from_string(string):
+ ret = []
+ for num in string.split(',') :
+ ret.append(float(num))
+ assert len(ret) > 0
+ return numpy.array(ret)
+
+def read_data(ifile):
+ "ifile can be a filename string or open (seekable) file object"
+ unlabeled_data=scipy.io.read_array(file)
+ return unlabeled_data
+
+def get_array(string, filename, name) :
+ "get an array from supplied command line options"
+ if string != None :
+ array = array_from_string(string)
+ elif filename != None :
+ array = read_data(filename)
+ else :
+ raise Exception, "no %s information given" % (name)
+ return array
+
+if __name__ == '__main__' :
+ # command line interface
+ from optparse import OptionParser
+
+ usage_string = ('%prog <bumps> <temps> <vibs>\n'
+ '2008, W. Trevor King.\n'
+ '\n'
+ 'Takes arrays of Vphotodiode sensitivity (V/nm), Temperature (K), \n'
+ 'and Vibration variance (V**2) as comma seperated lists.\n'
+ 'Returns the cantilever spring constant (pN/nm).\n'
+ 'for example:\n'
+ ' $ %prog -b 0.02,0.03,0.025 -t 298.2,300.1 -v 6e-9,5.5e-9\n'
+ )
+ parser = OptionParser(usage=usage_string, version='%prog '+common.VERSION)
+ parser.add_option('-b', '--bump-string', dest='bump_string',
+ help='comma seperated photodiode sensitivities (V/nm)',
+ type='string', metavar='BUMPS')
+ parser.add_option('-t', '--temp-string', dest='temp_string',
+ help='comma seperated temperatures (K)',
+ type='string', metavar='TEMPS')
+ parser.add_option('-v', '--vib-string', dest='vib_string',
+ help='comma seperated vibration variances (V**2)',
+ type='string', metavar='VIBS')
+ parser.add_option('-B', '--bump-file', dest='bump_file',
+ help='comma seperated photodiode sensitivities (V/nm)',
+ type='string', metavar='BUMPFILE')
+ parser.add_option('-T', '--temp-file', dest='temp_file',
+ help='comma seperated temperatures (K)',
+ type='string', metavar='TEMPFILE')
+ parser.add_option('-V', '--vib-file', dest='vib_file',
+ help='comma seperated vibration variances (V**2)',
+ type='string', metavar='VIBFILE')
+ parser.add_option('-C', '--celsius', dest='celsius',
+ help='Use Celsius input temperatures instead of Kelvin (defaul %default)\n',
+ action='store_true', default=False)
+ parser.add_option('-o', '--output-file', dest='ofilename',
+ help='write output to FILE (default stdout)',
+ type='string', metavar='FILE')
+ parser.add_option('-p', '--plot-inputs', dest='plot',
+ help='plot the input arrays to check their distribution',
+ action='store_true', default=False)
+ parser.add_option('', '--verbose', dest='verbose', action='store_true',
+ help='print lots of debugging information',
+ default=False)
+
+ options,args = parser.parse_args()
+ parser.destroy()
+
+ bumps = get_array(options.bump_string, options.bump_file, "bump")
+ temps = get_array(options.temp_string, options.temp_file, "temp")
+ vibs = get_array(options.vib_string, options.vib_file, "vib")
+
+ if options.plot :
+ calib_plot(bumps, temps, vibs)
+
+ if options.celsius :
+ for i in range(len(temps)) :
+ temps[i] = T_analyze.C_to_K(temps[i])
+
+ km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s = calib_analyze(bumps, temps, vibs)
+
+ if options.verbose :
+ print >> sys.stderr, \
+ string_errors(km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s)
+
+ if options.ofilename != None :
+ print >> file(options.ofilename, 'w'), km
+ else :
+ print km
--- /dev/null
+#!/usr/bin/python
+#
+# calibcant - tools for thermally calibrating AFM cantilevers
+#
+# Copyright (C) 2007,2008, William Trevor King
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as
+# published by the Free Software Foundation; either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+# See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
+# 02111-1307, USA.
+#
+# The author may be contacted at <wking@drexel.edu> on the Internet, or
+# write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
+# Philadelphia PA 19104, USA.
+
+"""
+Aquire, save, and load cantilever calibration bump data.
+For measuring photodiode sensitivity.
+
+W. Trevor King Dec. 2007 - Oct. 2008
+
+The relevent physical quantities are :
+ Vzp_out Output z-piezo voltage (what we generate)
+ Vzp Applied z-piezo voltage (after external ZPGAIN)
+ Zp The z-piezo position
+ Zcant The cantilever vertical deflection
+ Vphoto The photodiode vertical deflection voltage (what we measure)
+
+Which are related by the parameters :
+ zpGain Vzp_out / Vzp
+ zpSensitivity Zp / Vzp
+ photoSensitivity Vphoto / Zcant
+
+Cantilever calibration assumes a pre-calibrated z-piezo (zpSensitivity) and
+amplifier (zpGain). In our lab, the z-piezo is calibrated by imaging a
+calibration sample, which has features with well defined sizes, and the gain
+is set with a knob on the Nanoscope.
+
+photoSensitivity is measured by bumping the cantilever against the surface,
+where Zp = Zcant (see the bump_*() family of functions)
+The measured slope Vphoto/Vout is converted to photoSensitivity via
+Vphoto/Vzp_out * Vzp_out/Vzp * Vzp/Zp * Zp/Zcant = Vphoto/Zcant
+ (measured) (1/zpGain) (1/zpSensitivity) (1) (photoSensitivity)
+
+We do all these measurements a few times to estimate statistical errors.
+
+The functions are layed out in the families:
+ bump_*()
+For each family, * can be any of :
+ aquire get real-world data
+ save store real-world data to disk
+ load get real-world data from disk
+ analyze interperate the real-world data.
+ plot show a nice graphic to convince people we're working :p
+ load_analyze_tweaked
+ read a file with a list of paths to previously saved real world data
+ load each file using *_load(), analyze using *_analyze(), and
+ optionally plot using *_plot().
+ Intended for re-processing old data.
+A family name without any _* extension (e.g. bump()),
+ runs *_aquire(), *_save(), *_analyze(), *_plot().
+"""
+
+import numpy
+import time
+import data_logger
+import z_piezo_utils
+import FFT_tools
+import linfit
+from calibcant_bump_analyze import bump_analyze
+
+LOG_DATA = True # quietly grab all real-world data and log to LOG_DIR
+LOG_DIR = '/home/wking/rsrch/data/calibrate_cantilever'
+
+TEXT_VERBOSE = True # for debugging
+
+
+# bump family
+
+def bump_aquire(zpiezo, push_depth, npoints, freq) :
+ """
+ Ramps closer push_depth and returns to the original position.
+ Inputs:
+ zpiezo an opened zpiezo.zpiezo instance
+ push_depth distance to approach, in nm
+ npoints number points during the approach and during the retreat
+ freq rate at which data is aquired
+ log_dir directory to log data to (see data_logger.py).
+ None to turn off logging (see also the global LOG_DATA).
+ Returns the aquired ramp data dictionary, with data in DAC/ADC bits.
+ """
+ # generate the bump output
+ start_pos = zpiezo.curPos()
+ pos_dist = zpiezo.pos_nm2out(push_depth) - zpiezo.pos_nm2out(0)
+ close_pos = start_pos + pos_dist
+ appr = linspace(start_pos, close_pos, npoints)
+ retr = linspace(close_pos, start_pos, npoints)
+ out = concatenate((appr, retr))
+ # run the bump, and measure deflection
+ if TEXT_VERBOSE :
+ print "Bump %g nm" % push_depth
+ data = zpiezo.ramp(out, freq)
+ # default saving, so we have a log in-case the operator is lazy ;)
+ if LOG_DATA == True :
+ log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
+ log_name="bump_surface")
+ log.write_dict_of_arrays(data)
+ return data
+
+def bump_save(data, log_dir) :
+ "Save the dictionary data, using data_logger.data_log()"
+ if log_dir != None :
+ log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
+ log_name="bump")
+ log.write_dict_of_arrays(data)
+
+def bump_load(datafile) :
+ "Load the dictionary data, using data_logger.date_load()"
+ dl = data_logger.data_load()
+ data = dl.read_dict_of_arrays(path)
+ return data
+
+def bump_plot(data, plotVerbose) :
+ "Plot the bump (Vphoto vs Vzp) if plotVerbose or PYLAB_VERBOSE == True"
+ if plotVerbose or PYLAB_VERBOSE :
+ _import_pylab()
+ _pylab.figure(BASE_FIGNUM)
+ _pylab.plot(data["Z piezo output"], data["Deflection input"],
+ '.', label='bump')
+ _pylab.title("bump surface")
+ _pylab.legend(loc='upper left')
+ _flush_plot()
+
+def bump(zpiezo, push_depth, npoints=1024, freq=100e3,
+ log_dir=None,
+ plotVerbose=False) :
+ """
+ Wrapper around bump_aquire(), bump_save(), bump_analyze(), bump_plot()
+ """
+ data = bump_aquire(zpiezo, push_depth, npoints, freq)
+ bump_save(data, log_dir)
+ photoSensitivity = bump_analyze(data, zpiezo.gain, zpiezo.sensitivity,
+ zpiezo.pos_out2V, zpiezo.def_in2V)
+ bump_plot(data, plotVerbose)
+ return photoSensitivity
+
+def bump_load_analyze_tweaked(tweak_file, zpGain=_usual_zpGain,
+ zpSensitivity=_usual_zpSensitivity,
+ Vzp_out2V=_usual_Vzp_out2V,
+ Vphoto_in2V=_usual_Vphoto_in2V,
+ plotVerbose=False) :
+ "Load the output file of tweak_calib_bump.sh, return an array of slopes"
+ photoSensitivity = []
+ for line in file(tweak_file, 'r') :
+ parsed = line.split()
+ path = parsed[0].split('\n')[0]
+ # read the data
+ full_data = bump_load(path)
+ if len(parsed) == 1 :
+ data = full_data # use whole bump
+ else :
+ # use the listed sections
+ zp = []
+ df = []
+ for rng in parsed[1:] :
+ p = rng.split(':')
+ starti = int(p[0])
+ stopi = int(p[1])
+ zp.extend(full_data['Z piezo output'][starti:stopi])
+ df.extend(full_data['Deflection input'][starti:stopi])
+ data = {'Z piezo output': array(zp),
+ 'Deflection input':array(df)}
+ pSi = bump_analyze(data, zpGain, zpSensitivity,
+ Vzp_out2V, Vphoto_in2V, plotVerbose)
+ photoSensitivity.append(pSi)
+ bump_plot(data, plotVervose)
+ return array(photoSensitivity, dtype=numpy.float)
+
--- /dev/null
+#!/usr/bin/python
+#
+# calibcant - tools for thermally calibrating AFM cantilevers
+#
+# Copyright (C) 2007,2008, William Trevor King
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as
+# published by the Free Software Foundation; either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+# See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
+# 02111-1307, USA.
+#
+# The author may be contacted at <wking@drexel.edu> on the Internet, or
+# write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
+# Philadelphia PA 19104, USA.
+
+"""
+Separate the more general bump_analyze() from the other bump_*()
+functions in calibcant. Also provide a command line interface
+for analyzing data acquired through other workflows.
+
+The relevant physical quantities are :
+ Vzp_out Output z-piezo voltage (what we generate)
+ Vzp Applied z-piezo voltage (after external ZPGAIN)
+ Zp The z-piezo position
+ Zcant The cantilever vertical deflection
+ Vphoto The photodiode vertical deflection voltage (what we measure)
+
+Which are related by the parameters :
+ zpGain Vzp_out / Vzp
+ zpSensitivity Zp / Vzp
+ photoSensitivity Vphoto / Zcant
+"""
+
+import numpy
+import common # common module for the calibcant package
+import config # config module for the calibcant package
+import data_logger
+import z_piezo_utils
+import linfit
+
+def bump_analyze(data, zpGain=config.zpGain,
+ zpSensitivity=config.zpSensitivity,
+ Vzp_out2V=config.Vzp_out2V,
+ Vphoto_in2V=config.Vphoto_in2V,
+ textVerboseFile=None, plotVerbose=False) :
+ """
+ Return the slope of the bump ;).
+ Inputs:
+ data dictionary of data in DAC/ADC bits
+ Vzp_out2V function that converts output DAC bits to Volts
+ Vphoto_in2V function that converts input ADC bits to Volts
+ zpGain zpiezo applied voltage per output Volt
+ zpSensitivity nm zpiezo response per applied Volt
+ Returns:
+ photoSensitivity (Vphoto/Zcant) in Volts/nm
+ Checks for strong correlation (r-value) and low randomness chance (p-value)
+
+ With the current implementation, the data is regressed in DAC/ADC bits
+ and THEN converted, so we're assuming that both conversions are LINEAR.
+ if they aren't, rewrite to convert before the regression.
+ """
+ scale_Vzp_bits2V = Vzp_out2V(1) - Vzp_out2V(0)
+ scale_Vphoto_bits2V = Vphoto_in2V(1) - Vphoto_in2V(0)
+ Vphoto2Vzp_out_bit, intercept = \
+ linfit.linregress(x=data["Z piezo output"],
+ y=data["Deflection input"],
+ plotVerbose=plotVerbose)
+ Vphoto2Vzp_out = Vphoto2Vzp_out_bit * scale_Vphoto_bits2V / scale_Vzp_bits2V
+
+ # 1 / (Vzp/Vzp_out * Zp/Vzp * Zcant/Zp )
+ Vzp_out2Zcant = 1.0/ (zpGain * zpSensitivity) # * 1
+ return Vphoto2Vzp_out * Vzp_out2Zcant
+
+def bump_save(data, log_dir) :
+ "Save the dictionary data, using data_logger.data_log()"
+ if log_dir != None :
+ log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
+ log_name="bump")
+ log.write_dict_of_arrays(data)
+
+def bump_load(datafile) :
+ "Load the dictionary data, using data_logger.date_load()"
+ dl = data_logger.data_load()
+ data = dl.read_dict_of_arrays(datafile)
+ return data
+
+def bump_plot(data, plotVerbose) :
+ "Plot the bump (Vphoto vs Vzp) if plotVerbose or PYLAB_VERBOSE == True"
+ if plotVerbose or config.PYLAB_VERBOSE :
+ common._import_pylab()
+ common._pylab.figure(config.BASE_FIGNUM)
+ common._pylab.plot(data["Z piezo output"], data["Deflection input"],
+ '.', label='bump')
+ common._pylab.title("bump surface")
+ common._pylab.legend(loc='upper left')
+ common._flush_plot()
+
+def bump_load_analyze_tweaked(tweak_file, zpGain=config.zpGain,
+ zpSensitivity=config.zpSensitivity,
+ Vzp_out2V=config.Vzp_out2V,
+ Vphoto_in2V=config.Vphoto_in2V,
+ textVerboseFile=None, plotVerbose=False) :
+ "Load the output file of tweak_calib_bump.sh, return an array of slopes"
+ photoSensitivity = []
+ for line in file(tweak_file, 'r') :
+ parsed = line.split()
+ path = parsed[0].split('\n')[0]
+ if textVerboseFile != None :
+ print >> textVerboseFile, "Reading data from %s with ranges %s" % (path, parsed[1:])
+ # read the data
+ full_data = bump_load(path)
+ if len(parsed) == 1 :
+ data = full_data # use whole bump
+ else :
+ # use the listed sections
+ zp = []
+ df = []
+ for rng in parsed[1:] :
+ p = rng.split(':')
+ starti = int(p[0])
+ stopi = int(p[1])
+ zp.extend(full_data['Z piezo output'][starti:stopi])
+ df.extend(full_data['Deflection input'][starti:stopi])
+ data = {'Z piezo output': numpy.array(zp),
+ 'Deflection input': numpy.array(df)}
+ pSi = bump_analyze(data, zpGain, zpSensitivity,
+ Vzp_out2V, Vphoto_in2V, plotVerbose)
+ photoSensitivity.append(pSi)
+ bump_plot(data, plotVerbose)
+ return numpy.array(photoSensitivity, dtype=numpy.float)
+
+# commandline interface functions
+import scipy.io, sys
+
+def read_data(ifile):
+ "ifile can be a filename string or open (seekable) file object"
+ if ifile == None : ifile = sys.stdin
+ unlabeled_data=scipy.io.read_array(ifile)
+ data = {}
+ data['Z piezo output'] = unlabeled_data[:,0]
+ data['Deflection input'] = unlabeled_data[:,1]
+ return data
+
+def remove_further_than(data, zp_crit) :
+ ndata = {}
+ ndata['Z piezo output'] = []
+ ndata['Deflection input'] = []
+ for zp,df in zip(data['Z piezo output'],data['Deflection input']) :
+ if zp > zp_crit :
+ ndata['Z piezo output'].append(zp)
+ ndata['Deflection input'].append(df)
+ return ndata
+
+if __name__ == '__main__' :
+ # command line interface
+ from optparse import OptionParser
+
+ usage_string = ('%prog <input-file>\n'
+ '2008, W. Trevor King.\n'
+ '\n'
+ 'There are two operation modes, one to analyze a single bump file,\n'
+ 'and one to analyze tweak files.\n'
+ '\n'
+ 'Single file mode (the default) :\n'
+ 'Scales raw DAC/ADC bit data and fits a straight line.\n'
+ 'Returns photodiode sensitivity Vphotodiode/Zcantilever in V/nm.\n'
+ '<input-file> should be whitespace-delimited, 2 column ASCII\n'
+ 'without a header line. e.g: "<zp_DAC>\\t<deflection_ADC>\\n"\n'
+ '\n'
+ 'Tweak file mode:\n'
+ 'Runs the same analysis as in single file mode for each bump in\n'
+ 'a tweak file. Each line in the tweak file specifies a single bump.\n'
+ 'The format of a line is a series of whitespace-separated fields--\n'
+ 'a base file path followed by optional point index ranges, e.g.:\n'
+ '20080919/20080919132500_bump_surface 10:651 1413:2047\n'
+ 'which only discards all points outside the index ranges [10,651)\n'
+ 'and [1413,2047) (indexing starts at 0).\n'
+ )
+ parser = OptionParser(usage=usage_string, version='%prog '+common.VERSION)
+ parser.add_option('-C', '--cut-contact', dest='cut',
+ help='bilinear fit to cut out contact region (currently only available in single-file mode)',
+ action='store_true', default=False)
+ parser.add_option('-o', '--output-file', dest='ofilename',
+ help='write output to FILE (default stdout)',
+ type='string', metavar='FILE')
+ parser.add_option('-c', '--comma-out', dest='comma_out', action='store_true',
+ help='Output comma-seperated values (default %default)',
+ default=False)
+ parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
+ help='Run in tweak-file mode',
+ default=False)
+ parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
+ help='Print lots of debugging information',
+ default=False)
+
+ options,args = parser.parse_args()
+ parser.destroy()
+ assert len(args) >= 1, "Need an input file"
+
+ ifilename = args[0]
+
+ if options.ofilename != None :
+ ofile = file(options.ofilename, 'w')
+ else :
+ ofile = sys.stdout
+ if options.verbose == True :
+ vfile = sys.stderr
+ else :
+ vfile = None
+ config.TEXT_VERBOSE = options.verbose
+ config.PYLAB_VERBOSE = False
+ config.GNUPLOT_VERBOSE = False
+
+ if options.tweakmode == False :
+ data = read_data(ifilename)
+ if options.cut :
+ ddict = {'approach':data} # although it may be any combination of approach and retract
+ try :
+ params = z_piezo_utils.analyzeSurfPosData(ddict, retAllParams=True)
+ a,b,c,d = tuple(params) # model : f(x) = x<c ? a + b*x : (a+b*c) + d*(x-c)
+ print >> sys.stderr, "fit with", params, ". using zp < %d" % c
+ data = remove_further_than(data, c)
+ except z_piezo_utils.poorFit, s :
+ # data not very bilinear, so don't cut anything off.
+ print >> sys.stderr, "use everything"
+
+ photoSensitivity = bump_analyze(data, textVerboseFile=vfile)
+
+ print >> ofile, photoSensitivity
+ else : # tweak file mode
+ slopes = bump_load_analyze_tweaked(ifilename, textVerboseFile=vfile)
+ if options.comma_out :
+ sep = ','
+ else :
+ sep = '\n'
+ common.write_array(ofile, slopes, sep)
+
+ if options.ofilename != None :
+ ofile.close()
--- /dev/null
+#!/usr/bin/python
+
+"""
+Aquire and analyze cantilever calibration data.
+
+W. Trevor King Dec. 2007-Jan. 2008
+
+GPL BOILERPLATE
+
+
+The relevent physical quantities are :
+ Vzp_out Output z-piezo voltage (what we generate)
+ Vzp Applied z-piezo voltage (after external ZPGAIN)
+ Zp The z-piezo position
+ Zcant The cantilever vertical deflection
+ Vphoto The photodiode vertical deflection voltage (what we measure)
+ Fcant The force on the cantilever
+ T The temperature of the cantilever and surrounding solution
+ (another thing we measure)
+ k_b Boltzmann's constant
+
+Which are related by the parameters :
+ zpGain Vzp_out / Vzp
+ zpSensitivity Zp / Vzp
+ photoSensitivity Vphoto / Zcant
+ k_cant Fcant / Zcant
+
+Cantilever calibration assumes a pre-calibrated z-piezo (zpSensitivity) and
+a amplifier (zpGain).
+In our lab, the z-piezo is calibrated by imaging a calibration sample,
+which has features with well defined sizes, and the gain is set with a knob
+on the Nanoscope.
+
+photoSensitivity is measured by bumping the cantilever against the surface,
+where Zp = Zcant (see the bump_*() family of functions)
+The measured slope Vphoto/Vout is converted to photoSensitivity via
+Vphoto/Vzp_out * Vzp_out/Vzp * Vzp/Zp * Zp/Zcant = Vphoto/Zcant
+ (measured) (1/zpGain) (1/zpSensitivity) (1) (photoSensitivity)
+
+k_cant is measured by watching the cantilever vibrate in free solution
+(see the vib_*() family of functions)
+The average energy of the cantilever in the vertical direction is given
+by the equipartition theorem.
+ 1/2 k_b T = 1/2 k_cant Zcant**2
+ so k_cant = k_b T / Zcant**2
+ but Zcant = Vphoto / photoSensitivity
+ so k_cant = k_b T * photoSensitivty**2 / Vphoto**2
+We measured photoSensitivity above.
+We can either measure T using an external function (see temperature.py),
+or just estimate it (see the T_*() family of functions).
+ (room temp ~22 deg C, accurate to +/- 5 deg, so 5/273.15 ~= 2%)
+Finally, a time series of Vphoto while we're far from the surface and not
+changing Vzp_out will give us Vphoto.
+
+We do all these measurements a few times to estimate statistical errors.
+
+The functions are layed out in the families:
+ bump_*(), vib_*(), T_*(), and calib_*()
+where calib_{save|load|analyze}() deal with derived data, not real-world data.
+For each family, * can be any of :
+ aquire get real-world data
+ save store real-world data to disk
+ load get real-world data from disk
+ analyze interperate the real-world data.
+ plot show a nice graphic to convince people we're working :p
+ load_analyze_tweaked
+ read a file with a list of paths to previously saved real world data
+ load each file using *_load(), analyze using *_analyze(), and
+ optionally plot using *_plot().
+ Intended for re-processing old data.
+A family name without any _* extension (e.g. bump()),
+ runs *_aquire(), *_save(), *_analyze(), *_plot().
+
+Also defines the two positioning functions:
+ move_just_onto_surface() and move_far_from_surface()
+"""
+
+import numpy
+import time
+import data_logger
+import z_piezo_utils
+import FFT_tools
+import linfit
+import GnuplotBiDir # used for fitting lorentzian
+
+kb = 1.3806504e-23 # Boltzmann's constant in J/K
+DEFAULT_TEMP = 22 # assume 22 deg C
+
+LOG_DATA = True # quietly grab all real-world data and log to LOG_DIR
+LOG_DIR = '/home/wking/rsrch/data/calibrate_cantilever'
+
+GNUFIT_DATA_BASE='./calibrate_cantilever_fitdata'
+
+TEXT_VERBOSE = True # for debugging
+GNUPLOT_VERBOSE = True # turn on fit check plotting
+PYLAB_VERBOSE = True # turn on plotting
+PYLAB_INTERACTIVE = True # select between draw() and show() for flushing plots
+BASE_FIGNUM = 20 # to avoid writing to already existing figures
+
+# handle extra verbose input modules, only imported if we need them
+_flush_plot = None
+_final_flush_plot = None
+_pylab = None
+def _import_pylab() :
+ "Import pylab plotting functions for when we need to plot"
+ global _pylab
+ global _flush_plot
+ global _final_flush_plot
+ if _pylab == None :
+ import pylab as _pylab
+ if _flush_plot == None :
+ if PYLAB_INTERACTIVE :
+ _flush_plot = _pylab.draw
+ else :
+ def _flush_plot () : pass
+ if _final_flush_plot == None :
+ if PYLAB_INTERACTIVE :
+ _final_flush_plot = _pylab.draw
+ else :
+ _final_flush_plot = _pylab.show
+
+
+# HACK
+# make sure you make a system note (add_system_note) if you change these
+# in case you don't have access to a z_piezo.z_piezo for conversion functions
+_usual_zpGain=20
+_usual_zpSensitivity=7.41
+_usual_zeroVphoto_bits = 2**15
+_usual_zeroVzp_bits = 2**15
+def _usual_Vphoto_in2V(inp) :
+ return (float(inp)-float(_usual_zeroVphoto_bits))*(10.0/2**15)
+def _usual_Vzp_out2V(inp) :
+ return (float(inp)-float(_usual_zeroVzp_bits))*(10.0/2**15)
+
+def C_to_K(celsius) :
+ "Convert Celsius -> Kelvin."
+ return celsius + 273.15
+
+# bump family
+
+def bump_aquire(zpiezo, push_depth, npoints, freq) :
+ """
+ Ramps closer push_depth and returns to the original position.
+ Inputs:
+ zpiezo an opened zpiezo.zpiezo instance
+ push_depth distance to approach, in nm
+ npoints number points during the approach and during the retreat
+ freq rate at which data is aquired
+ log_dir directory to log data to (see data_logger.py).
+ None to turn off logging (see also the global LOG_DATA).
+ Returns the aquired ramp data dictionary, with data in DAC/ADC bits.
+ """
+ # generate the bump output
+ start_pos = zpiezo.curPos()
+ pos_dist = zpiezo.pos_nm2out(push_depth) - zpiezo.pos_nm2out(0)
+ close_pos = start_pos + pos_dist
+ appr = linspace(start_pos, close_pos, npoints)
+ retr = linspace(close_pos, start_pos, npoints)
+ out = concatenate((appr, retr))
+ # run the bump, and measure deflection
+ if TEXT_VERBOSE :
+ print "Bump %g nm" % push_depth
+ data = zpiezo.ramp(out, freq)
+ # default saving, so we have a log in-case the operator is lazy ;)
+ if LOG_DATA == True :
+ log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
+ log_name="bump_surface")
+ log.write_dict_of_arrays(data)
+ return data
+
+def bump_save(data, log_dir) :
+ "Save the dictionary data, using data_logger.data_log()"
+ if log_dir != None :
+ log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
+ log_name="bump")
+ log.write_dict_of_arrays(data)
+
+def bump_load(datafile) :
+ "Load the dictionary data, using data_logger.date_load()"
+ dl = data_logger.data_load()
+ data = dl.read_dict_of_arrays(path)
+ return data
+
+def bump_analyze(data, zpGain=_usual_zpGain,
+ zpSensitivity=_usual_zpSensitivity,
+ Vzp_out2V=_usual_Vzp_out2V,
+ Vphoto_in2V=_usual_Vphoto_in2V,
+ plotVerbose=False) :
+ """
+ Return the slope of the bump ;).
+ Inputs:
+ data dictinary of data in DAC/ADC bits
+ pos_out2V function that converts output DAC bits to Volts
+ def_in2V funtion that converts input ADC bits to Volts
+ zpGain zpiezo applied voltage per output Volt
+ zpSensitivity nm zpiezo response per applied Volt
+ Returns:
+ photoSensitivity (Vphoto/Zcant) in Volts/nm
+ Checks for strong correlation (r-value) and low randomness chance (p-value)
+
+ With the current implementation, the data is regressed in DAC/ADC bits
+ and THEN converted, so we're assuming that both conversions are LINEAR.
+ if they aren't, rewrite to convert before the regression.
+ """
+ scale_Vzp_bits2V = Vzp_out2V(1) - Vzp_out2V(0)
+ scale_Vphoto_bits2V = Vzp_in2V(1) - Vphoto_in2V(0)
+ Vphoto2Vzp_out_bit, intercept = \
+ linfit.linregress(x=data["Z piezo output"],
+ y=data["Deflection input"],
+ plotVerbose=plotVerbose)
+ Vphoto2Vzp_out = Vphoto2Vzp_out * scale_Vphoto_bits2V / scale_Vzp_bits2V
+
+ # 1 / (Vzp/Vzp_out * Zp/Vzp * Zcant/Zp )
+ Vzp_out2Zcant = 1.0/ (zpiezo_gain * zpiezo_sensitivity) # * 1
+ return Vphoto2Vzp_out * Vzp_out2Zcant
+
+def bump_plot(data, plotVerbose) :
+ "Plot the bump (Vphoto vs Vzp) if plotVerbose or PYLAB_VERBOSE == True"
+ if plotVerbose or PYLAB_VERBOSE :
+ _import_pylab()
+ _pylab.figure(BASE_FIGNUM)
+ _pylab.plot(data["Z piezo output"], data["Deflection input"],
+ '.', label='bump')
+ _pylab.title("bump surface")
+ _pylab.legend(loc='upper left')
+ _flush_plot()
+
+def bump(zpiezo, push_depth, npoints=1024, freq=100e3,
+ log_dir=None,
+ plotVerbose=False) :
+ """
+ Wrapper around bump_aquire(), bump_save(), bump_analyze(), bump_plot()
+ """
+ data = bump_aquire(zpiezo, push_depth, npoints, freq)
+ bump_save(data, log_dir)
+ photoSensitivity = bump_analyze(data, zpiezo.gain, zpiezo.sensitivity,
+ zpiezo.pos_out2V, zpiezo.def_in2V)
+ bump_plot(data, plotVerbose)
+ return photoSensitivity
+
+def bump_load_analyze_tweaked(tweak_file, zpGain=_usual_zpGain,
+ zpSensitivity=_usual_zpSensitivity,
+ Vzp_out2V=_usual_Vzp_out2V,
+ Vphoto_in2V=_usual_Vphoto_in2V,
+ plotVerbose=False) :
+ "Load the output file of tweak_calib_bump.sh, return an array of slopes"
+ photoSensitivity = []
+ for line in file(tweak_file, 'r') :
+ parsed = line.split()
+ path = parsed[0].split('\n')[0]
+ # read the data
+ full_data = bump_load(path)
+ if len(parsed) == 1 :
+ data = full_data # use whole bump
+ else :
+ # use the listed sections
+ zp = []
+ df = []
+ for rng in parsed[1:] :
+ p = rng.split(':')
+ starti = int(p[0])
+ stopi = int(p[1])
+ zp.extend(full_data['Z piezo output'][starti:stopi])
+ df.extend(full_data['Deflection input'][starti:stopi])
+ data = {'Z piezo output': array(zp),
+ 'Deflection input':array(df)}
+ pSi = bump_analyze(data, zpGain, zpSensitivity,
+ Vzp_out2V, Vphoto_in2V, plotVerbose)
+ photoSensitivity.append(pSi)
+ bump_plot(data, plotVervose)
+ return array(photoSensitivity, dtype=numpy.float)
+
+# T family.
+# Fairly stubby, since a one shot Temp measurement is a common thing.
+# We just wrap that to provide a consistent interface.
+
+def T_aquire(get_T=None) :
+ """
+ Measure the current temperature of the sample,
+ or, if get_T == None, fake it by returning DEFAULT_TEMP
+ """
+ if get_T == None :
+ if TEXT_VERBOSE :
+ print "Fake temperature %g" % DEFAULT_TEMP
+ return DEFAULT_TEMP
+ else :
+ if TEXT_VERBOSE :
+ print "Measure temperature"
+ return get_T()
+
+def T_save(T, log_dir=None) :
+ """
+ Save either a single T (if you are crazy :p),
+ or an array of them (more reasonable)
+ """
+ T = numpy.array(T, dtype=numpy.float)
+ if log_dir != None :
+ log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
+ log_name="T_float")
+ log.write_binary(T.tostring())
+ if LOG_DATA != None :
+ log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
+ log_name="T_float")
+ log.write_binary(T.tostring())
+
+def T_load(datafile=None) :
+ """
+ Load the saved T array (possibly of length 1), and return it.
+ If datafile == None, return an array of one DEFAULT_TEMP instead.
+ """
+ if datafile == None :
+ return numpy.array([DEFAULT_TEMP], dtype=numpy.float)
+ else :
+ return fromfile(file=cleanpath, dtype=numpy.float)
+
+def T_analyze(T, convert_to_K=C_to_K) :
+ """
+ Not much to do here, just convert to Kelvin.
+ Uses C_to_K (defined above) by default.
+ """
+ try : # if T is an array, convert element by element
+ for i in range(len(T)) :
+ T[i] = convert_to_K(T[i])
+ except TypeError :
+ T = convert_to_K(T)
+ return T
+
+def T_plot(T, plotVerbose) :
+ """
+ Umm, just print the temperature?
+ """
+ if plotVerbose or PYLAB_VERBOSE or TEXT_VERBOSE :
+ print "Temperature ", T
+
+def T(get_T=None, convert_to_K=C_to_K, log_dir=None, plotVerbose=None) :
+ """
+ Wrapper around T_aquire(), T_save(), T_analyze(), T_plot()
+ """
+ T_raw = T_aquire(get_T)
+ T_save(T_raw, log_dir)
+ T_ret = T_analyze(T_raw, convert_to_K)
+ T_plot(T_raw, plotVerbose)
+ return T_ret
+
+def T_load_analyze_tweaked(tweak_file=None, convert_to_K=C_to_K,
+ plotVerbose=False) :
+ "Read the T files listed in tweak_file, return an array of Ts in Kelvin"
+ if tweak_file != None :
+ Tlist=[]
+ for filepath in file(datafile, 'r') :
+ cleanpath = filepath.split('\n')[0]
+ Tlist.extend(T_load(cleanpath))
+ Tlist = numpy.array(Tlist, dtype=numpy.float)
+ T_raw = array(Tlist, dtype=numpy.float)
+ for Ti in T_raw :
+ T_plot(T, plotVerbose)
+ else :
+ T_raw = T_load(None)
+ T = T_analyze(T_raw, convert_to_K)
+ return T
+
+# vib family
+
+def vib_aquire(zpiezo, time, freq) :
+ """
+ Record data for TIME seconds at FREQ Hz from ZPIEZO at it's current position.
+ """
+ # round up to the nearest power of two, for efficient FFT-ing
+ nsamps = ceil_pow_of_two(time*freq)
+ time = nsamps / freq
+ # take some data, keeping the position voltage at it's current value
+ out = ones((nsamps,), dtype=uint16) * zpiezo.curPos()
+ if TEXT_VERBOSE :
+ print "get %g seconds of data" % time
+ data = zpiezo.ramp(out, freq)
+ if LOG_DATA :
+ log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
+ log_name="vib_%gHz" % freq)
+ log.write_dict_of_arrays(data)
+ return data
+
+def vib_save(data, log_dir=None) :
+ "Save the dictionary data, using data_logger.data_log()"
+ if log_dir :
+ log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
+ log_name="vib_%gHz" % freq)
+ log.write_dict_of_arrays(data)
+
+def vib_load(datafile=None) :
+ "Load the dictionary data, using data_logger.date_load()"
+ dl = data_logger.data_load()
+ data = dl.read_dict_of_arrays(path)
+ return data
+
+def vib_plot(data, freq,
+ chunk_size=2048, overlap=True, window=FFT_tools.window_hann,
+ plotVerbose=True) :
+ """
+ If plotVerbose or PYLAB_VERBOSE == True, plot 3 subfigures:
+ Time series (Vphoto vs sample index) (show any major drift events),
+ A histogram of Vphoto, with a gaussian fit (demonstrate randomness), and
+ FFTed Vphoto data (Vphoto vs Frequency) (show major noise components).
+ """
+ if plotVerbose or PYLAB_VERBOSE :
+ _import_pylab()
+ _pylab.hold(False)
+ _pylab.figure(BASE_FIGNUM+2)
+
+ # plot time series
+ _pylab.subplot(311)
+ _pylab.plot(data["Deflection input"], 'r.')
+ _pylab.title("free oscillation")
+
+ # plot histogram distribution and gaussian fit
+ _pylab.subplot(312)
+ n, bins, patches = \
+ _pylab.hist(data["Deflection input"], bins=30,
+ normed=1, align='center')
+ gaus = numpy.zeros((len(bins),))
+ mean = data["Deflection input"].mean()
+ std = data["Deflection input"].std()
+ pi = numpy.pi
+ exp = numpy.exp
+ for i in range(len(bins)) :
+ gaus[i] = (2*pi)**(-.5)/std * exp(-0.5*((bins[i]-mean)/std)**2)
+ _pylab.hold(True)
+ _pylab.plot(bins, gaus, 'r-');
+ _pylab.hold(False)
+
+ # plot FFTed data
+ _pylab.subplot(313)
+ AC_data = data["Deflection input"] - mean
+ (freq_axis, power) = \
+ FFT_tools.unitary_avg_power_spectrum(AC_data, freq, chunk_size,
+ overlap, window)
+ _pylab.semilogy(freq_axis, power, 'r.-')
+ _pylab.hold(True)
+ filtered_power = FFT_tools.windowed_filter(power,
+ FFT_tools.windowed_median_point, s=5)
+ _pylab.semilogy(freq_axis, filtered_power, 'b.-')
+ _flush_plot()
+
+def vib_analyze_naive(data, Vphoto_in2V=_usual_Vphoto_in2V,
+ zeroVphoto_bits=_usual_zeroVphoto_bits,
+ plotVerbose=False) :
+ """
+ Calculate the variance of the raw data, and convert it to Volts**2.
+ This method is simple and easy to understand, but it highly succeptible
+ to noise, drift, etc.
+
+ Vphoto_in2V is a function converting Vphoto values from bits to Volts.
+ This function is assumed linear, see bump_analyze().
+ zeroVphoto_bits is the bit value corresponding to Vphoto = zero Volts.
+ """
+ Vphoto_std_bits = data["Deflection input"].std()
+ Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zpiezo_zeroV_in)
+ return Vphoto_std**2
+
+def vib_analyze(data, freq, minFreq=500, maxFreq=7000,
+ chunk_size=4096, overlap=True, window=FFT_tools.window_hann,
+ median_filter_width=5,
+ Vphoto_in2V=_usual_Vphoto_in2V, plotVerbose=False) :
+ """
+ Calculate the variance in the raw data due to the cantilever's
+ thermal oscillation and convert it to Volts**2.
+
+ Improves on vib_analyze_naive() by first converting Vphoto(t) to
+ frequency space, and fitting a Lorentzian in the relevant frequency
+ range (see cantilever_calib.pdf for derivation).
+
+ The conversion to frequency space generates an average power spectrum
+ by breaking the data into windowed chunks and averaging the power
+ spectrums for the chunks together.
+ See FFT_tools.avg_power_spectrum
+
+ Vphoto_in2V is a function converting Vphoto values from bits to Volts.
+ This function is assumed linear, see bump_analyze().
+ zeroVphoto_bits is the bit value corresponding to Vphoto = zero Volts.
+ """
+ Vphoto_t = numpy.zeros((len(data["Deflection input"]),),
+ dtype=numpy.float)
+ # convert the data from bits to volts
+ if TEXT_VERBOSE :
+ print "Converting %d bit values to voltages" % len(Vphoto_t)
+ for i in range(len(data["Deflection input"])) :
+ Vphoto_t[i] = Vphoto_in2V(data["Deflection input"][i])
+
+ # Compute the average power spectral density per unit time (in uV**2/Hz)
+ if TEXT_VERBOSE :
+ print "Compute the averaged power spectral density in uV**2/Hz"
+ (freq_axis, power) = \
+ FFT_tools.unitary_avg_power_spectrum((Vphoto_t - Vphoto_t.mean())*1e6,
+ freq, chunk_size,overlap, window)
+
+ # cut out the relevent frequency range
+ if TEXT_VERBOSE :
+ print "Cut the frequency range %g to %g Hz" % (minFreq, maxFreq)
+ imin = 0
+ while freq_axis[imin] < minFreq : imin += 1
+ imax = imin
+ while freq_axis[imax] < maxFreq : imax += 1
+ assert imax >= imin + 10 , \
+ "Less than 10 points in freq range (%g,%g)" % (minFreq, maxFreq)
+
+ # We don't need to filter, because taking data is so cheap...
+ # median filter the relevent range
+ #filtered_power = FFT_tools.windowed_filter(power[imin:imax],
+ # FFT_tools.windowed_median_point,
+ # s=median_filter_width)
+ filtered_power = power[imin:imax]
+
+ # check
+ if (plotVerbose or PYLAB_VERBOSE) and False :
+ if TEXT_VERBOSE :
+ print "Plot the FFT data for checking"
+ vib_plot(data, freq, chunk_size, overlap, window, plotVerbose)
+ _import_pylab()
+ _pylab.figure(5)
+ #_pylab.semilogy(freq_axis, power, 'r.-')
+ #print "len ", len(freq_axis), len(power), freq_axis.min(), filtered_power.min()
+ _pylab.semilogy(freq_axis, power, 'r.-',
+ freq_axis[imin:imax], power[imin:imax], 'b.-',
+ freq_axis[imin:imax], filtered_power, 'g.-')
+ _flush_plot()
+
+ # save about-to-be-fitted data to a temporary file
+ if TEXT_VERBOSE :
+ print "Save data in temp file for fitting"
+ datafilename = "%s_%d" % (GNUFIT_DATA_BASE, time.time())
+ fd = file(datafilename, 'w')
+ for i in range(imin, imax) :
+ fd.write("%g\t%g\n" % (freq_axis[i], filtered_power[i-imin]))
+ fd.close()
+
+ # generate guesses for Lorentzian parameters A,B,C
+ if TEXT_VERBOSE :
+ print "Fit lorentzian"
+ max_fp_index = numpy.argmin(filtered_power)
+ max_fp = filtered_power[max_fp_index]
+ half_max_index = 0
+ while filtered_power[half_max_index] < max_fp/2 :
+ half_max_index += 1
+ # Lorentzian L(x) = c / ((a**2-x**2)**2 + (b*x)**2)
+ # peak centered at a, so
+ A_guess = freq_axis[max_fp_index+imin]
+ # Half width w on lower side when L(a-w) = L(a)/2
+ # (a**2 - (a-w)**2)**2 + (b*(a-w))**2 = 2*(b*a)**2
+ # Let W=(a-w)**2, A=a**2, and B=b**2
+ # (A - W)**2 + BW = 2*AB
+ # W**2 - 2AW + A**2 + BW = 2AB
+ # W**2 + (B-2A)W + (A**2-2AB) = 0
+ # W = (2A-B)/2 * [1 +/- sqrt(1 - 4(A**2-2AB)/(B-2A)**2]
+ # = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
+ # (a-w)**2 = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
+ # so w is a disaster ;)
+ # For some values of A and B (non-underdamped), W is imaginary or negative.
+ B_guess = A_guess/2
+ # for underdamped, max value at x ~= a, where f(a) = c/(ab)**2, so
+ C_guess = max_fp * (A_guess*B_guess)**2
+
+ # fit Lorentzian using Gnuplot's 'fit' command
+ g = GnuplotBiDir.Gnuplot()
+ # The Lorentzian is the predicted one-sided power spectral density
+ # For an oscillator whose position x obeys
+ # m x'' + gamma x' + kx = F_thermal(t)
+ # A is the ?driving noise?
+ # B is gamma, the damping term
+ # C is omega_0, the resonant frequency
+ # Fitting the PSD of V = photoSensitivity*x just rescales the parameters
+ g("f(x) = C / ((A**2-x**2)**2 + (x*B)**2)")
+ ## The mass m is given by m = k_B T / (pi A)
+ ## The spring constant k is given by k = m (omega_0)**2
+ ## The quality factor Q is given by Q = omega_0 m / gamma
+ g("A=%g; B=%g; C=%g" % (A_guess, B_guess, C_guess))
+ g("fit f(x) '%s' via A,B,C" % datafilename)
+ A=float(g.getvar('A'))
+ B=float(g.getvar('B'))
+ C=float(g.getvar('C'))
+ if TEXT_VERBOSE :
+ print "Fitted f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with"
+ print "A = %g, \t B = %g, \t C = %g" % (A, B, C)
+ if plotVerbose or GNUPLOT_VERBOSE :
+ if TEXT_VERBOSE :
+ print "Fitted f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with"
+ print "A = %g, \t B = %g, \t C = %g" % (A, B, C)
+ # write all the ft data now
+ fd = file(datafilename, 'w')
+ for i in range(len(freq_axis)) :
+ fd.write("%g\t%g\n" % (freq_axis[i], power[i]))
+ fd.write("\n") # blank line to separate fit data.
+ # write the fit data
+ for i in range(imin, imax) :
+ x = freq_axis[i]
+ fd.write("%g\t%g\n" % (freq_axis[i],
+ C / ((A**2-x**2)**2 + (x*B)**2) ) )
+ fd.close()
+ print Vphoto_t.var(), A, B, C, ((numpy.pi * C)/(2 * B * A**2) * 1e-12)
+ g("set terminal epslatex color solid")
+ g("set output 'calibration.tex")
+ g("set size 2,2") # multipliers on default 5"x3.5"
+ #g("set title 'Thermal calibration'")
+ g("set logscale y")
+ g("set xlabel 'Frequency (Hz)'")
+ g("set ylabel 'Power ($\mu$V$^2$/Hz)'")
+ g("set xrange [0:20000]")
+ g("g(x) = x < %g || x > %g ? 1/0 : f(x)" % (minFreq, maxFreq)) # only plot on fitted region
+ g("plot '%s' title 'Vibration data', g(x) title 'Lorentzian fit'" % datafilename)
+
+ g("set mouse")
+ g("pause mouse 'click with mouse'")
+ g.getvar("MOUSE_BUTTON")
+ del(g)
+
+ # Integrating the the power spectral density per unit time (power) over the
+ # frequency axis [0, fN] returns the total signal power per unit time
+ # int_0^fN power(f)df = 1/T int_0^T |x(t)**2| dt
+ # = <V(t)**2>, the variance for our AC signal.
+ # The variance from our fitted Lorentzian is the area under the Lorentzian
+ # <V(t)**2> = (pi*C) / (2*B*A**2)
+ # Our A is in uV**2, so convert back to Volts
+ return (numpy.pi * C) / (2 * B * A**2) * 1e-12
+
+def vib(zpiezo, time=1, freq=50e3, log_dir=None, plotVerbose=False) :
+ """
+ Wrapper around vib_aquire(), vib_save(), vib_analyze(), vib_plot()
+ """
+ data = vib_aquire(zpiezo, time, freq)
+ vib_save(data, log_dir)
+ Vphoto_var = vib_analyze(data, freq)
+ vib_plot(data, plotVerbose)
+ return Vphoto_var
+
+def vib_load_analyze_tweaked(tweak_file, minFreq=1000, maxFreq=7000,
+ chunk_size=2048, overlap=True,
+ window=FFT_tools.window_hann,
+ median_filter_width=5,
+ Vphoto_in2V=_usual_Vphoto_in2V,
+ plotVerbose=False) :
+ """
+ Read the vib files listed in tweak_file.
+ Return an array of Vphoto variances (due to the cantilever) in Volts**2
+ """
+ dl = data_logger.data_load()
+ Vphoto_var = []
+ for path in file(tweak_file, 'r') :
+ if path[-1] == '\n':
+ path = path.split('\n')[0]
+ # read the data
+ data = dl.read_dict_of_arrays(path)
+ freq = float(path.split('_')[-1].split('Hz')[0])
+ if TEXT_VERBOSE and False :
+ print "Analyzing '%s' at %g Hz" % (path, freq)
+ # analyze
+ Vphoto_var.append(vib_analyze(data, freq, minFreq, maxFreq,
+ chunk_size, overlap, window,
+ median_filter_width,
+ Vphoto_in2V, plotVerbose))
+ return numpy.array(Vphoto_var, dtype=numpy.float)
+
+
+# A few positioning functions, so we can run bump_aquire() and vib_aquire()
+# with proper spacing relative to the surface.
+
+def move_just_onto_surface(stepper, zpiezo, Depth_nm=100, setpoint=2) :
+ """
+ Uses z_piezo_utils.getSurfPos() to pinpoint the position of the surface.
+ Adjusts the stepper position as required to get within stepper_tol nm
+ of the surface.
+ Then set Vzp to place the cantilever Depth_nm onto the surface.
+
+ If getSurfPos() fails to find the surface, backs off (for safety)
+ and steps in (without moving the zpiezo) until Vphoto > setpoint.
+ """
+ stepper_tol = 250 # nm, generous estimate of the fullstep stepsize
+
+ if TEXT_VERBOSE :
+ print "moving just onto surface"
+ # Zero the piezo
+ if TEXT_VERBOSE :
+ print "zero the z piezo output"
+ zpiezo.jumpToPos(zpiezo.pos_nm2out(0))
+ # See if we're near the surface already
+ if TEXT_VERBOSE :
+ print "See if we're starting near the surface"
+ try :
+ dist = zpiezo.pos_out2nm( \
+ z_piezo_utils.getSurfPos(zpiezo, zpiezo.def_V2in(setpoint))
+ )
+ except (z_piezo_utils.tooClose, z_piezo_utils.poorFit), string :
+ if TEXT_VERBOSE :
+ print "distance failed with: ", string
+ print "Back off 200 half steps"
+ # Back away 200 steps
+ stepper.step_rel(-400)
+ stepper.step_rel(200)
+ sp = zpiezo.def_V2in(setpoint) # sp = setpoint in bits
+ zpiezo.updateInputs()
+ cd = zpiezo.curDef() # cd = current deflection in bits
+ if TEXT_VERBOSE :
+ print "Single stepping approach"
+ while cd < sp :
+ if TEXT_VERBOSE :
+ print "deflection %g < setpoint %g. step closer" % (cd, sp)
+ stepper.step_rel(2) # Full step in
+ zpiezo.updateInputs()
+ cd = zpiezo.curDef()
+ # Back off two steps (protecting against backlash)
+ if TEXT_VERBOSE :
+ print "Step back 4 half steps to get off the setpoint"
+ stepper.step_rel(-200)
+ stepper.step_rel(196)
+ # get the distance to the surface
+ zpiezo.updateInputs()
+ if TEXT_VERBOSE :
+ print "get surf pos, with setpoint %g (%d)" % (setpoint, zpiezo.def_V2in(setpoint))
+ for i in range(20) : # HACK, keep stepping back until we get a distance
+ try :
+ dist = zpiezo.pos_out2nm(getSurfPos(zpiezo, zpiezo.def_V2in(setpoint)))
+ except (tooClose, poorFit), string :
+ stepper.step_rel(-200)
+ stepper.step_rel(198)
+ continue
+ break
+ if i >= 19 :
+ print "tried %d times, still too close! bailing" % i
+ print "probably an invalid setpoint."
+ raise Exception, "weirdness"
+ if TEXT_VERBOSE :
+ print 'distance to surface ', dist, ' nm'
+ # fine tune the stepper position
+ while dist < -stepper_tol : # step back if we need to
+ stepper.step_rel(-200)
+ stepper.step_rel(198)
+ dist = zpiezo.pos_out2nm(getSurfPos(zpiezo, zpiezo.def_V2in(setpoint)))
+ if TEXT_VERBOSE :
+ print 'distance to surface ', dist, ' nm, step back'
+ while dist > stepper_tol : # and step forward if we need to
+ stepper.step_rel(2)
+ dist = zpiezo.pos_out2nm(getSurfPos(zpiezo, zpiezo.def_V2in(setpoint)))
+ if TEXT_VERBOSE :
+ print 'distance to surface ', dist, ' nm, step closer'
+ # now adjust the zpiezo to place us just onto the surface
+ target = dist + Depth_nm
+ zpiezo.jumpToPos(zpiezo.pos_nm2out(target))
+ # and we're there :)
+ if TEXT_VERBOSE :
+ print "We're %g nm into the surface" % Depth_nm
+
+def move_far_from_surface(stepper, um_back=50) :
+ """
+ Step back a specified number of microns.
+ (uses very rough estimate of step distance at the moment)
+ """
+ step_nm = 100
+ steps = int(um_back*1000/step_nm)
+ print "step back %d steps" % steps
+ stepper.step_rel(-steps)
+
+
+# and finally, the calib family
+
+def calib_aquire(stepper, zpiezo, get_T=None,
+ approach_setpoint=2,
+ bump_start_depth=100, bump_push_depth=200,
+ bump_npoints=1024, bump_freq=100e3,
+ T_convert_to_K=C_to_K,
+ vib_time=1, vib_freq=50e3,
+ num_bumps=10, num_Ts=10, num_vibs=20,
+ log_dir=None, plotVerbose=False) :
+ """
+ Aquire data for calibrating a cantilever in one function.
+ return (bump, T, vib), each of which is an array.
+ Inputs :
+ stepper a stepper.stepper_obj for coarse Z positioning
+ zpiezo a z_piezo.z_piezo for fine positioning and deflection readin
+ setpoint maximum allowed deflection (in Volts) during approaches
+ num_bumps number of 'a's (see Outputs)
+ push_depth_nm depth of each push when generating a
+ num_temps number of 'T's (see Outputs)
+ num_vibs number of 'vib's (see Outputs)
+ log_dir directory to log data to. Default 'None' disables logging.
+ Outputs (all are arrays of recorded data) :
+ bumps measured (V_photodiode / nm_tip) proportionality constant
+ Ts measured temperature (K)
+ vibs measured V_photodiode variance in free solution
+ """
+ # get bumps
+ move_just_onto_surface(stepper, zpiezo,
+ bump_start_depth, approach_setpoint)
+ bumps=zeros((num_bumps,))
+ for i in range(num_bumps) :
+ bumps[i] = bump(zpiezo, bump_push_depth, bump_npoints, bump_freq,
+ log_dir, plotVerbose)
+ if TEXT_VERBOSE :
+ print bumps
+
+ # get Ts
+ if get_T == None :
+ get_T = lambda : DEFAULT_TEMP
+ assert T_convert_to_K == C_to_K, "You didn't define get_T()!"
+ Ts=zeros((num_Ts,), dtype=float)
+ for i in range(num_Ts) :
+ Ts[i] = T(get_T, T_convert_to_K, log_dir, plotVerbose)
+ time.sleep(1) # wait a bit to get an independent temperature measure
+ print Ts
+
+ # get vibs
+ move_far_from_surface(stepper)
+ vibs=zeros((num_vibs,))
+ for i in range(num_vibs) :
+ vibs[i] = vib(zpiezo, vib_time, vib_freq, log_dir=log_dir)
+ print vibs
+
+ if LOG_DATA != None :
+ data = {'bump':bumps, 'T':Ts, 'vib':vibs}
+ log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
+ log_name="calib")
+ log.write_dict_of_arrays(data)
+
+ return (bumps, Ts, vibs)
+
+def calib_save(bumps, Ts, vibs, log_dir) :
+ """
+ Save a dictonary with the bump, T, and vib data.
+ """
+ if log_dir != None :
+ data = {'bump':bumps, 'T':Ts, 'vib':vibs}
+ log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
+ log_name="calib")
+ log.write_dict_of_arrays(data)
+
+def calib_load(datafile) :
+ "Load the dictionary data, using data_logger.date_load()"
+ dl = data_logger.data_load()
+ data = dl.read_dict_of_arrays(path)
+ return (data['bump'], data['T'], data['vib'])
+
+def calib_save_analysis(k, k_s,
+ photoSensitivity2_m, photoSensitivity2_s,
+ T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s,
+ log_dir=None) :
+ if log_dir != None :
+ log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
+ log_name="calib_analysis_text")
+ log.write_binary((
+ "k (N/m) : %g +/- %g\n" % (k, k_s)
+ + "photoSensitivity**2 (V/nm)**2 : %g +/- %g\n" % \
+ (photoSensitivity2_m, photoSensitivity2_s)
+ + "T (K) : %g +/- %g\n" % (T_m, T_s)
+ + "1/Vp**2 (1/V)**2 : %g +/- %g\n" % \
+ (one_o_Vphoto2_m, one_o_Vphoto2_s)
+ ))
+
+def calib_plot(bumps, Ts, vibs, plotVerbose) :
+ if plotVerbose or PYLAB_VERBOSE :
+ _import_pylab()
+ _pylab.figure(BASE_FIGNUM+4)
+ _pylab.subplot(311)
+ _pylab.plot(bumps, 'g.')
+ _pylab.title('Photodiode sensitivity (V/nm)')
+ _pylab.subplot(312)
+ _pylab.plot(Ts, 'r.')
+ _pylab.title('Temperature (K)')
+ _pylab.subplot(313)
+ _pylab.plot(vibs, 'b.')
+ _pylab.title('Thermal deflection variance (Volts**2)')
+ _flush_plot()
+
+def calib(stepper, zpiezo, tempController=None,
+ setpoint=2.0,
+ num_bumps=10, push_depth_nm=300,
+ num_temps=10,
+ num_vibs=10,
+ log_dir=None) :
+ """
+ Calibrate a cantilever in one function.
+ The I-don't-care-about-the-details black box version :p.
+ return (k, k_s)
+ Inputs:
+ (see calib_aquire())
+ Outputs :
+ k cantilever spring constant (in N/m, or equivalently nN/nm)
+ k_s standard deviation in our estimate of k
+ Notes :
+ See get_calibration_data() for the data aquisition code
+ See analyze_calibration_data() for the analysis code
+ """
+ a, T, vib = calib_aquire(stepper, zpiezo, tempController,
+ setpoint,
+ num_bumps=num_bumps,
+ push_depth_nm=push_depth_nm,
+ num_temps=num_temps,
+ num_vibs=num_vibs,
+ log_dir=log_dir)
+ calib_save(a, T, vib, log_dir)
+ k,k_s,ps2_m, ps2_s,T_m,T_s,one_o_Vp2_m,one_o_Vp2_s = \
+ calib_analyze(a, T, vib, log_dir=log_dir)
+ calib_save_analysis(k,k_s,ps2_m, ps2_s,T_m,T_s,one_o_Vp2_m,one_o_Vp2_s,
+ log_dir)
+ calib_plot(a, T, vib, plotVerbose)
+ return (k, k_s)
+
+def calib_load_analyze_tweaks(bump_tweaks, vib_tweaks,
+ T_tweaks=None,
+ log_dir=None,
+ plotVerbose=True) :
+ a = read_tweaked_bumps(bump_tweaks)
+ vib = V_photo_variance_from_file(vib_tweaks)
+ return analyze_calibration_data(a, T, vib, log_dir=log_dir)
+
+
--- /dev/null
+import config
+
+VERSION="0.2"
+
+# handle extra verbose input modules, only imported if we need them
+_flush_plot = None
+_final_flush_plot = None
+_pylab = None
+
+def _import_pylab() :
+ """Import pylab plotting functions for when we need to plot. This
+ function can be called multiple times, and ensures that the pylab
+ setup is only imported once. It defines the functions
+ _flush_plot() to be called after incremental changes,
+ _final_flush_plot(), to be called at the end of any graphical processing, and
+ _pylab, a reference to the pylab module."""
+ global _pylab
+ global _flush_plot
+ global _final_flush_plot
+ if _pylab == None :
+ import pylab as _pylab
+ if _flush_plot == None :
+ if config.PYLAB_INTERACTIVE :
+ _flush_plot = _pylab.draw
+ else :
+ def _flush_plot () : pass
+ if _final_flush_plot == None :
+ if config.PYLAB_INTERACTIVE :
+ _final_flush_plot = _pylab.draw
+ else :
+ _final_flush_plot = _pylab.show
+
+def write_array(ofile, array, seperator):
+ """Erite an iterable array seperated by seperator.
+ Terminate the output with an endline."""
+ strings = [repr(x) for x in array]
+ ofile.write(seperator.join(strings))
+ ofile.write('\n')
--- /dev/null
+"""Define some variables to configure the package for a particular lab
+and workflow."""
+
+import z_piezo_utils
+
+DEFAULT_TEMP = 22 # assume 22 deg C
+LOG_DATA = True # quietly grab all real-world data and log to LOG_DIR
+LOG_DIR = '/home/wking/rsrch/data/calibrate_cantilever'
+GNUFIT_DATA_BASE='./calibrate_cantilever_fitdata'
+TEXT_VERBOSE = True # for debugging
+GNUPLOT_VERBOSE = True # turn on fit check plotting
+PYLAB_VERBOSE = True # turn on plotting
+PYLAB_INTERACTIVE = True # select between draw() and show() for flushing plots
+BASE_FIGNUM = 20 # to avoid writing to already existing figures
+
+# zpGain zpiezo applied voltage per output Volt
+zpGain = z_piezo_utils._usual_zpGain
+# zpSensitivity nm zpiezo response per applied Volt
+zpSensitivity = z_piezo_utils._usual_zpSensitivity
+# Vzp_out2V function that converts output DAC bits to Volts
+Vzp_out2V = z_piezo_utils._usual_Vzp_out2V
+# Vphoto_in2V function that converts input ADC bits to Volts
+Vphoto_in2V = z_piezo_utils._usual_Vphoto_in2V
+# zeroVphoto_bits ADC bit output for a 0V input
+zeroVphoto_bits = z_piezo_utils._usual_zeroVphoto_bits
--- /dev/null
+#!/usr/bin/python
+
+"""
+windowed_filter(array, winfunc, s)
+with winfuncs :
+ windowed_mean_point
+ windowed_median_point
+
+Command line bindings for filtering 2 column ACSII data.
+"""
+
+import numpy
+
+def windowed_mean_point(array, i, width=50) :
+ """
+ Filter data with a windowed average (mean).
+ The ith point is replaced by the average of the points with indices in
+ the range i +/- width (inclusive).
+
+ This function returns the ith point in the filtered array,
+ given the raw array and window half-width.
+
+ The weights on the average are uniform at the moment.
+ """
+ x=0
+ n=0
+ for j in range(i-width, i+width) :
+ if j >= 0 and j < len(array) :
+ x+=array[j]
+ n+=1
+ x /= float(n)
+ return x
+
+def windowed_median_point(array, i, width=9) :
+ """
+ Filter data with a windowed median.
+ The ith point is replaced by the median of the points with indices in
+ the range i +/- width (inclusive).
+
+ This function returns the ith point in the filtered array,
+ given the raw array and window half-width.
+ """
+ imin = i-width
+ if imin < 0 : imin = 0
+ imax = i+width
+ if imax >= len(array) : imax = len(array-1)
+ slice = array[imin:imax+1].copy()
+ slice.sort()
+ imid = numpy.floor((imax-imin)/2)
+ return slice[imid]
+
+def windowed_filter(array, winfunc, width=None) :
+ """
+ Filter data with a windowing function winfunc.
+ The ith point is replaced by the winfunc(array, i, s).
+
+ See the windowed_* functions for possible winfunc options.
+ """
+ out = array.copy()
+ if width == None :
+ win = lambda i : winfunc(array, i) # user winfunc's default s
+ else :
+ win = lambda i : winfunc(array, i, width)
+ for i in range(len(out)) :
+ out[i] = win(i)
+ return out
+
+
+# commandline interface functions
+import scipy.io, sys
+
+def read_data(ifile):
+ """
+ ifile can be a filename string or open (seekable) file object.
+ returns (column 1 array, column 2 array)
+ """
+ if ifile == None : ifile = sys.stdin
+ data=scipy.io.read_array(ifile)
+ return (data[:,0], data[:,1])
+
+def write_data(ofile, x, y):
+ """
+ ofile can be a filename string or open (seekable) file object.
+ """
+ data = numpy.zeros((len(x),2))
+ data[:,0] = x
+ data[:,1] = y
+ if ofile == None : ofile = sys.stdout
+ scipy.io.write_array(ofile, data, separator='\t')
+
+if __name__ == '__main__' :
+ # command line interface
+ from optparse import OptionParser
+
+ usage_string = ('%prog <input-file>\n'
+ '2008, W. Trevor King.\n'
+ '\n'
+ 'Apply various filters to data y(x).'
+ '<input-file> should be whitespace-delimited, 2 column ASCII\n'
+ 'without a header line.\n'
+ 'e.g: "<x>\\t<y>\\n"')
+ parser = OptionParser(usage=usage_string, version='%prog 0.1')
+ parser.add_option('-o', '--output-file', dest='ofilename',
+ help='write output to FILE (default stdout)',
+ type='string', metavar='FILE')
+ parser.add_option('-w', '--width', dest='width', default=None,
+ help='window width (i +/- width, inclusive)',
+ type='int', metavar='WIDTH')
+ parser.add_option('-t', '--type', dest='type', default='mean',
+ help='filter type (default %default)',
+ type='string', metavar='TYPE')
+ parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
+ help='Print lots of debugging information',
+ default=False)
+
+ options,args = parser.parse_args()
+ parser.destroy()
+ assert len(args) >= 1, "Need an input file"
+
+ ifilename = args[0]
+
+ x,y = read_data(ifilename)
+
+ if options.type == 'mean' :
+ winfunc = windowed_mean_point
+ elif options.type == 'median' :
+ winfunc = windowed_median_point
+ else :
+ raise Exception, "unrecognized window type '%s'" % options.type
+
+ y = windowed_filter(y, winfunc, width=options.width)
+
+ write_data(options.ofilename, x, y)
--- /dev/null
+#!/usr/bin/python
+
+if __name__ == "__main__" :
+
+ import sys
+ import filter, calibcant_vib_analyze
+ from optparse import OptionParser
+
+ usage_string = ('%prog <input-file>\n'
+ '2008, W. Trevor King.\n'
+ '\n'
+ 'Apply various filters to data y(x).'
+ '<input-file> should be whitespace-delimited, 2 column ASCII\n'
+ 'without a header line.\n'
+ 'e.g: "<x>\\t<y>\\n"')
+ parser = OptionParser(usage=usage_string, version='%prog 0.1')
+ parser.add_option('-o', '--output-file', dest='ofilename',
+ help='write output to FILE (default stdout)',
+ type='string', metavar='FILE')
+ parser.add_option('-w', '--width', dest='width', default=None,
+ help='window width (i +/- width, inclusive)',
+ type='int', metavar='WIDTH')
+ parser.add_option('-t', '--type', dest='type', default='mean',
+ help='filter type (default %default)',
+ type='string', metavar='TYPE')
+ parser.add_option('-m', '--min-freq', dest='min_freq', default=500,
+ help='minimum frequency in Hz (default %default)',
+ type='float', metavar='FREQ')
+ parser.add_option('-M', '--max-freq', dest='max_freq', default=7000,
+ help='maximum frequency in Hz (default %default)',
+ type='float', metavar='FREQ')
+ parser.add_option('-g', '--gnuplot', dest='gnuplot', action='store_true',
+ help='Print gnuplot fit check script to stderr',
+ default=False)
+ parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
+ help='Print lots of debugging information',
+ default=False)
+
+ options,args = parser.parse_args()
+ parser.destroy()
+ assert len(args) >= 1, "Need an input file"
+
+ ifilename = args[0]
+
+ # mirror filter behavior:
+ freq,psd = filter.read_data(ifilename)
+
+ if options.type == 'mean' :
+ winfunc = filter.windowed_mean_point
+ elif options.type == 'median' :
+ winfunc = filter.windowed_median_point
+ else :
+ raise Exception, "unrecognized window type '%s'" % options.type
+
+ psd = filter.windowed_filter(psd, winfunc, width=options.width)
+
+ # mirror calibcant_vib_analyze behavior:
+ (A,B,C) = calibcant_vib_analyze.fit_psd(freq, psd,
+ minFreq=options.min_freq,
+ maxFreq=options.max_freq)
+ if options.verbose :
+ print >> sys.stderr, "Fit f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with"
+ print >> sys.stderr, "A = %g, \t B = %g, \t C = %g" % (A, B, C)
+ if options.gnuplot :
+ print >> sys.stderr, \
+ calibcant_vib_analyze.gnuplot_check_fit(ifilename, A, B, C)
+
+ area = calibcant_vib_analyze.lorentzian_area(A,B,C)
+
+ if options.ofilename != None :
+ print >> file(options.ofilename, 'w'), area
+ else :
+ print area
--- /dev/null
+#!/usr/bin/python
+#
+# calibcant - tools for thermally calibrating AFM cantilevers
+#
+# Copyright (C) 2007,2008, William Trevor King
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as
+# published by the Free Software Foundation; either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+# See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
+# 02111-1307, USA.
+#
+# The author may be contacted at <wking@drexel.edu> on the Internet, or
+# write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
+# Philadelphia PA 19104, USA.
+
+"""
+Separate the more general vib_analyze() from the other vib_*()
+functions in calibcant. Also provide a command line interface
+for analyzing data acquired through other workflows.
+
+The relevent physical quantities are :
+ Vphoto The photodiode vertical deflection voltage (what we measure)
+"""
+
+import os, time, numpy
+import GnuplotBiDir # used for fitting lorentzian
+import common # common module for the calibcant package
+import config # config module for the calibcant package
+import data_logger
+import FFT_tools
+
+def vib_analyze_naive(deflection_bits, Vphoto_in2V=config.Vphoto_in2V,
+ zeroVphoto_bits=config.zeroVphoto_bits) :
+ """
+ Calculate the variance of the raw data, and convert it to Volts**2.
+ This method is simple and easy to understand, but it highly succeptible
+ to noise, drift, etc.
+
+ Vphoto_in2V is a function converting Vphoto values from bits to Volts.
+ This function is assumed linear, see bump_analyze().
+ zeroVphoto_bits is the bit value corresponding to Vphoto = zero Volts.
+ """
+ Vphoto_std_bits = deflection_bits.std()
+ Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zpiezo_zeroV_in)
+ return Vphoto_std**2
+
+def vib_analyze(deflection_bits, freq, minFreq=500, maxFreq=7000,
+ chunk_size=4096, overlap=True, window=FFT_tools.window_hann,
+ Vphoto_in2V=config.Vphoto_in2V, plotVerbose=False) :
+ """
+ Calculate the variance in the raw data due to the cantilever's
+ thermal oscillation and convert it to Volts**2.
+
+ Improves on vib_analyze_naive() by first converting Vphoto(t) to
+ frequency space, and fitting a Lorentzian in the relevant frequency
+ range (see cantilever_calib.pdf for derivation).
+
+ The conversion to frequency space generates an average power spectrum
+ by breaking the data into windowed chunks and averaging the power
+ spectrums for the chunks together.
+ See FFT_tools.avg_power_spectrum
+
+ Vphoto_in2V is a function converting Vphoto values from bits to Volts.
+ This function is assumed linear, see bump_analyze().
+ zeroVphoto_bits is the bit value corresponding to Vphoto = zero Volts.
+ """
+ Vphoto_t = numpy.zeros((len(deflection_bits),),
+ dtype=numpy.float)
+ # convert the data from bits to volts
+ if config.TEXT_VERBOSE :
+ print "Converting %d bit values to voltages" % len(Vphoto_t)
+ for i in range(len(deflection_bits)) :
+ Vphoto_t[i] = Vphoto_in2V(deflection_bits[i])
+
+ # Compute the average power spectral density per unit time (in uV**2/Hz)
+ if config.TEXT_VERBOSE :
+ print "Compute the averaged power spectral density in uV**2/Hz"
+ (freq_axis, power) = \
+ FFT_tools.unitary_avg_power_spectrum((Vphoto_t - Vphoto_t.mean())*1e6,
+ freq, chunk_size,overlap, window)
+
+ A,B,C = fit_psd(freq_axis, power, minFreq, maxFreq)
+
+ if config.TEXT_VERBOSE :
+ print "Fitted f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with"
+ print "A = %g, \t B = %g, \t C = %g" % (A, B, C)
+
+ vib_plot(deflection_bits, freq_axis, power, A, B, C, plotVerbose=plotVerbose)
+
+ # Our A is in uV**2, so convert back to Volts
+ return lorentzian_area(A,B,C) * 1e-12
+
+def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) :
+ """
+ freq_axis : array of frequencies in Hz
+ psd_data : array of PSD amplitudes for the frequencies in freq_axis
+
+ Calculate the variance in the raw data due to the cantilever's
+ thermal oscillation and convert it to Volts**2.
+
+ Improves on vib_analyze_naive() by working on frequency space data
+ and fitting a Lorentzian in the relevant frequency range (see
+ cantilever_calib.pdf for derivation).
+
+ The conversion to frequency space generates an average power spectrum
+ by breaking the data into windowed chunks and averaging the power
+ spectrums for the chunks together.
+ See FFT_tools.unitary_avg_power_spectrum().
+ """
+ # cut out the relevent frequency range
+ if config.TEXT_VERBOSE :
+ print "Cut the frequency range %g to %g Hz" % (minFreq, maxFreq)
+ imin = 0
+ while freq_axis[imin] < minFreq : imin += 1
+ imax = imin
+ while freq_axis[imax] < maxFreq : imax += 1
+ assert imax >= imin + 10 , \
+ "Less than 10 points in freq range (%g,%g)" % (minFreq, maxFreq)
+
+ # save about-to-be-fitted data to a temporary file
+ datafilename = "%s_%d" % (config.GNUFIT_DATA_BASE, time.time())
+ fd = file(datafilename, 'w')
+ for i in range(imin, imax) :
+ fd.write("%g\t%g\n" % (freq_axis[i], psd_data[i]))
+ fd.close()
+
+ # generate guesses for Lorentzian parameters A,B,C
+ max_psd_index = numpy.argmax(psd_data[imin:imax]) + imin
+ max_psd = psd_data[max_psd_index]
+ half_max_index = imin
+ while psd_data[half_max_index] < max_psd/2 :
+ half_max_index += 1
+ res_freq = freq_axis[max_psd_index]
+ half_width = freq_axis[max_psd_index] - freq_axis[half_max_index]
+ # Lorentzian L(x) = C / ((A**2-x**2)**2 + (B*x)**2)
+ # is expected power spectrum for
+ # x'' + B x' + A^2 x'' = F_external(t)/m
+ # (A = omega_0)
+ # C = (2 k_B T B) / (pi m)
+ #
+ # A = resonant frequency
+ # peak at x_res = sqrt(A^2 - B^2/2)
+ # which could be complex if there isn't a peak (overdamped)
+ # peak height = C / (3 x_res^4 + A^4)
+ # Q = A/B
+ #
+ # guess Q = 5, and adjust from there
+ Q_guess = 5
+ # so x_res^2 = B^2 Q^2 - B^2/2 = (Q^2-1/2)B^2
+ # B = x_res / sqrt(Q^2-1/2)
+ if config.TEXT_VERBOSE :
+ print "params : %g, %g" % (res_freq, max_psd)
+ B_guess = res_freq / numpy.sqrt(Q_guess**2-0.5)
+ A_guess = Q_guess*B_guess
+ C_guess = max_psd * (-res_freq**4 + A_guess**4)
+ #
+ # Half width w on lower side when L(a-w) = L(a)/2
+ # (a**2 - (a-w)**2)**2 + (b*(a-w))**2 = 2*(b*a)**2
+ # Let W=(a-w)**2, A=a**2, and B=b**2
+ # (A - W)**2 + BW = 2*AB
+ # W**2 - 2AW + A**2 + BW = 2AB
+ # W**2 + (B-2A)W + (A**2-2AB) = 0
+ # W = (2A-B)/2 * [1 +/- sqrt(1 - 4(A**2-2AB)/(B-2A)**2]
+ # = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
+ # (a-w)**2 = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
+ # so w is a disaster ;)
+ # For some values of A and B (non-underdamped), W is imaginary or negative.
+ if config.TEXT_VERBOSE :
+ print "guesses : %g, %g, %g" % (A_guess, B_guess, C_guess)
+
+ # fit Lorentzian using Gnuplot's 'fit' command
+ g = GnuplotBiDir.Gnuplot()
+ # The Lorentzian is the predicted one-sided power spectral density
+ # For an oscillator whose position x obeys
+ # m x'' + gamma x' + kx = F_thermal(t)
+ # A is the ?driving noise?
+ # B is gamma, the damping term
+ # C is omega_0, the resonant frequency
+ # Fitting the PSD of V = photoSensitivity*x just rescales the parameters
+ g("f(x) = C / ((A**2-x**2)**2 + (B*x)**2)")
+ ## The mass m is given by m = k_B T / (pi A)
+ ## The spring constant k is given by k = m (omega_0)**2
+ ## The quality factor Q is given by Q = omega_0 m / gamma
+ g("A=%g; B=%g; C=%g" % (A_guess, B_guess, C_guess))
+ g("fit f(x) '%s' via A,B,C" % datafilename)
+ A=abs(float(g.getvar('A'))) # A and B only show up as squares in f(x)
+ B=abs(float(g.getvar('B'))) # so ensure we get positive values
+ C=float(g.getvar('C'))
+ os.remove(datafilename)
+ return (A, B, C)
+
+def lorentzian_area(A, B, C) :
+ # Integrating the the power spectral density per unit time (power) over the
+ # frequency axis [0, fN] returns the total signal power per unit time
+ # int_0^fN power(f)df = 1/T int_0^T |x(t)**2| dt
+ # = <V(t)**2>, the variance for our AC signal.
+ # The variance from our fitted Lorentzian is the area under the Lorentzian
+ # <V(t)**2> = (pi*C) / (2*B*A**2)
+ return (numpy.pi * C) / (2 * B * A**2)
+
+def gnuplot_check_fit(datafilename, A, B, C) :
+ """
+ return a string containing a gnuplot script displaying the fit.
+ """
+ string = 'f(x) = C / ((A**2-x**2)**2 + (x*B)**2)\n'
+ string += 'A= %g ; B= %g ; C= %g\n' % (A, B, C)
+ string += 'set logscale y\n'
+ string += "plot '%s', f(x)\n" % datafilename
+ return string
+
+def vib_save(data, log_dir=None) :
+ """Save the dictionary data, using data_logger.data_log()
+ data should be a dictionary of arrays with the fields
+ 'Z piezo output'
+ 'Z piezo input'
+ 'Deflection input'
+ """
+ if log_dir :
+ log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
+ log_name="vib_%gHz" % freq)
+ log.write_dict_of_arrays(data)
+
+def vib_load(datafile=None) :
+ """Load the dictionary data, using data_logger.date_load()"
+ data should be a dictionary of arrays with the fields
+ 'Z piezo output'
+ 'Z piezo input'
+ 'Deflection input'
+ """
+ dl = data_logger.data_load()
+ data = dl.read_dict_of_arrays(datafile)
+ return data
+
+def vib_plot(deflection_bits, freq_axis, power, A, B, C,
+ plotVerbose=True) :
+ """
+ If plotVerbose or config.PYLAB_VERBOSE == True, plot 3 subfigures:
+ Time series (Vphoto vs sample index) (show any major drift events),
+ A histogram of Vphoto, with a gaussian fit (demonstrate randomness), and
+ FFTed Vphoto data (Vphoto vs Frequency) (show major noise components).
+ """
+ if plotVerbose or config.PYLAB_VERBOSE :
+ common._import_pylab()
+ common._pylab.hold(False)
+ common._pylab.figure(BASE_FIGNUM+2)
+
+ # plot time series
+ common._pylab.subplot(311)
+ common._pylab.plot(data["Deflection input"], 'r.')
+ common._pylab.title("free oscillation")
+
+ # plot histogram distribution and gaussian fit
+ common._pylab.subplot(312)
+ n, bins, patches = \
+ common._pylab.hist(data["Deflection input"], bins=30,
+ normed=1, align='center')
+ gaus = numpy.zeros((len(bins),))
+ mean = data["Deflection input"].mean()
+ std = data["Deflection input"].std()
+ pi = numpy.pi
+ exp = numpy.exp
+ for i in range(len(bins)) :
+ gaus[i] = (2*pi)**(-.5)/std * exp(-0.5*((bins[i]-mean)/std)**2)
+ common._pylab.hold(True)
+ common._pylab.plot(bins, gaus, 'r-');
+ common._pylab.hold(False)
+
+ # plot FFTed data
+ common._pylab.subplot(313)
+ common._pylab.semilogy(freq_axis, power, 'r.-')
+ common._flush_plot()
+ if (plotVerbose or config.GNUPLOT_VERBOSE) and False : # TODO, clean up and double check...
+ # write all the ft data now
+ fd = file(datafilename, 'w')
+ for i in range(len(freq_axis)) :
+ fd.write("%g\t%g\n" % (freq_axis[i], power[i]))
+ fd.write("\n") # blank line to separate fit data.
+ # write the fit data
+ for i in range(imin, imax) :
+ x = freq_axis[i]
+ fd.write("%g\t%g\n" % (freq_axis[i],
+ C / ((A**2-x**2)**2 + (x*B)**2) ) )
+ fd.close()
+ print Vphoto_t.var(), A, B, C, ((numpy.pi * C)/(2 * B * A**2) * 1e-12)
+ g("set terminal epslatex color solid")
+ g("set output 'calibration.tex'")
+ g("set size 2,2") # multipliers on default 5"x3.5"
+ #g("set title 'Thermal calibration'")
+ g("set logscale y")
+ g("set xlabel 'Frequency (Hz)'")
+ g("set ylabel 'Power ($\mu$V$^2$/Hz)'")
+ g("set xrange [0:20000]")
+ g("g(x) = x < %g || x > %g ? 1/0 : f(x)" % (minFreq, maxFreq)) # only plot on fitted region
+ g("plot '%s' title 'Vibration data', g(x) title 'Lorentzian fit'" % datafilename)
+
+ g("set mouse")
+ g("pause mouse 'click with mouse'")
+ g.getvar("MOUSE_BUTTON")
+ del(g)
+
+def vib_load_analyze_tweaked(tweak_file, minFreq=1000, maxFreq=7000,
+ chunk_size=2048, overlap=True,
+ window=FFT_tools.window_hann,
+ Vphoto_in2V=config.Vphoto_in2V,
+ plotVerbose=False) :
+ """
+ Read the vib files listed in tweak_file.
+ Return an array of Vphoto variances (due to the cantilever) in Volts**2
+ """
+ dl = data_logger.data_load()
+ Vphoto_var = []
+ for path in file(tweak_file, 'r') :
+ if path[-1] == '\n':
+ path = path.split('\n')[0]
+ # read the data
+ data = vib_load(path)
+ deflection_bits = data['Deflection input']
+ freq = float(path.split('_')[-1].split('Hz')[0])
+ if config.TEXT_VERBOSE :
+ print "Analyzing '%s' at %g Hz" % (path, freq)
+ # analyze
+ Vphoto_var.append(vib_analyze(deflection_bits, freq, minFreq, maxFreq,
+ chunk_size, overlap, window,
+ Vphoto_in2V, plotVerbose))
+ return numpy.array(Vphoto_var, dtype=numpy.float)
+
+# commandline interface functions
+import scipy.io, sys
+
+def read_data(ifile):
+ """
+ ifile can be a filename string or open (seekable) file object.
+ returns an array of data values (1 column)
+ """
+ #returns (column 1 array, column 2 array)
+ #"""
+ if ifile == None : ifile = sys.stdin
+ unlabeled_data=scipy.io.read_array(ifile)
+ return unlabeled_data #(unlabeled_data[:,0], unlabeled_data[:,1])
+
+if __name__ == '__main__' :
+ # command line interface
+ from optparse import OptionParser
+
+ usage_string = ('%prog <input-file>\n'
+ '2008, W. Trevor King.\n'
+ '\n'
+ 'Deflection power spectral density (PSD) data (X^2/Hz)\n'
+ 'and returns the variance in X units (X^2)'
+ '<input-file> should be whitespace-delimited, 2 column ASCII\n'
+ 'without a header line.\n'
+ 'e.g: "<frequency_Hz>\\t<deflection_psd X^2/Hz>\\n"\n')
+ parser = OptionParser(usage=usage_string, version='%prog 0.1')
+ parser.add_option('-f', '--sample-freq', dest='freq', default=100e3,
+ help='sample frequency in Hz (default %default)',
+ type='float', metavar='FREQ')
+ parser.add_option('-m', '--min-freq', dest='min_freq', default=500,
+ help='minimum frequency in Hz (default %default)',
+ type='float', metavar='FREQ')
+ parser.add_option('-M', '--max-freq', dest='max_freq', default=7000,
+ help='maximum frequency in Hz (default %default)',
+ type='float', metavar='FREQ')
+ parser.add_option('-o', '--output-file', dest='ofilename',
+ help='write output to FILE (default stdout)',
+ type='string', metavar='FILE')
+ parser.add_option('-c', '--comma-out', dest='comma_out', action='store_true',
+ help='Output comma-seperated values (default %default)',
+ default=False)
+ parser.add_option('-g', '--gnuplot', dest='gnuplot', action='store_true',
+ help='Print gnuplot fit check script to stderr',
+ default=False)
+ parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
+ help='Run in tweak-file mode',
+ default=False)
+ parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
+ help='Print lots of debugging information',
+ default=False)
+
+ options,args = parser.parse_args()
+ parser.destroy()
+ assert len(args) >= 1, "Need an input file"
+
+ ifilename = args[0]
+
+ if options.ofilename != None :
+ ofile = file(options.ofilename, 'w')
+ else :
+ ofile = sys.stdout
+ if options.verbose == True :
+ vfile = sys.stderr
+ else :
+ vfile = None
+ config.TEXT_VERBOSE = options.verbose
+ config.PYLAB_VERBOSE = False
+ config.GNUPLOT_VERBOSE = options.gnuplot
+
+ if options.tweakmode == False :
+ data = read_data(ifilename)
+ deflection_bits = data['Deflection input']
+ Vphoto_var = vib_analyze(deflection_bits, freq=options.freq,
+ minFreq=options.min_freq,
+ maxFreq=options.max_freq)
+ print >> ofile, Vphoto_var
+ else : # tweak mode
+ Vphoto_vars = vib_load_analyze_tweaked(ifilename,
+ minFreq=options.min_freq,
+ maxFreq=options.max_freq)
+ if options.comma_out :
+ sep = ','
+ else :
+ sep = '\n'
+ common.write_array(ofile, Vphoto_vars, sep)
+
+ if options.ofilename != None :
+ ofile.close()
projects / calibcant.git / commitdiff