From 21b604e13776d3b10e2471416cdbfe9f2a3deeb9 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Tue, 11 Nov 2008 10:51:25 -0500 Subject: [PATCH 1/1] Began versioning. --- Makefile | 19 + README | 50 ++ calibcant/T_analyze.py | 186 ++++++ calibcant/T_analyze.pyc | Bin 0 -> 6127 bytes calibcant/__init__.py | 0 calibcant/analyze.py | 238 ++++++++ calibcant/bump.py | 188 ++++++ calibcant/bump_analyze.py | 249 ++++++++ calibcant/calibrate_cantilever.py | 911 ++++++++++++++++++++++++++++++ calibcant/common.py | 38 ++ calibcant/common.pyc | Bin 0 -> 1618 bytes calibcant/config.py | 25 + calibcant/config.pyc | Bin 0 -> 866 bytes calibcant/filter.py | 133 +++++ calibcant/psd_filter_analyze.py | 73 +++ calibcant/vib_analyze.py | 425 ++++++++++++++ 16 files changed, 2535 insertions(+) create mode 100644 Makefile create mode 100644 README create mode 100755 calibcant/T_analyze.py create mode 100644 calibcant/T_analyze.pyc create mode 100644 calibcant/__init__.py create mode 100755 calibcant/analyze.py create mode 100644 calibcant/bump.py create mode 100755 calibcant/bump_analyze.py create mode 100644 calibcant/calibrate_cantilever.py create mode 100644 calibcant/common.py create mode 100644 calibcant/common.pyc create mode 100644 calibcant/config.py create mode 100644 calibcant/config.pyc create mode 100644 calibcant/filter.py create mode 100755 calibcant/psd_filter_analyze.py create mode 100755 calibcant/vib_analyze.py diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..007b7d5 --- /dev/null +++ b/Makefile @@ -0,0 +1,19 @@ +.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 diff --git a/README b/README new file mode 100644 index 0000000..bf5945b --- /dev/null +++ b/README @@ -0,0 +1,50 @@ +== 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 + +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 + +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) diff --git a/calibcant/T_analyze.py b/calibcant/T_analyze.py new file mode 100755 index 0000000..4b94002 --- /dev/null +++ b/calibcant/T_analyze.py @@ -0,0 +1,186 @@ +#!/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 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 \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() diff --git a/calibcant/T_analyze.pyc b/calibcant/T_analyze.pyc new file mode 100644 index 0000000000000000000000000000000000000000..83a378148a6c25f45179dff96cecc732f20e0d7d GIT binary patch literal 6127 zcmcIoOLH5?5$;`*APJKa^`f4RMzUprG6_<&U6nYpY+9ydS_-8CAlpGX*;-(R;EIb~ zXlFngD(VvX59FBqg&h6Qa>+6Ix(6Uo_9X|I5(e|=>6z)r*WLK}-;3qr7ytRXO~b$Q z`2GSPdv%IP6E)EWq(eFgsELm6htv#5{V8frjr!BnoF4UOs5vw0m#A4HJw-o(;VhY1 z($n-47|tO0a5@;$T_$}+=&lOgHPY{r{vGM-q}3$d4bnG-?!6HwaDD&R{BDu{KTxB3gZ%Xkz00C`=Z=2dywHJ%B~yJW7?8wh2MUd&LhOo79ltMs=J zi+IU8U2a|{bA=Mc+ITTV`7%~+jHWj5G-}?Y3_HydCGYWo$=@e)ON74Y7e1ioZOUdL zanRo$kxWrkVL2LfRz{sWvic$A%yX8StHNF16@9MC;u_PkO0RNB*cs*ajq2uwQ5}<+B2HB!inBN!yfT#)wO8a_ zp~yj2VCha}rM%bA+Ahg6tCCE$<1}ft(<)6e15#&- zy|``4dwHS!2}#ybI(D&&+b{b`VKj~^@_uJu`QwiA;$Sb$k8HJEZs3}QNzEaSQoa3w zO>nY$*~c((iLoli6Fx1ghH98@4-avDOvFwd?I-PhenZj4hQp+VC0t&0O>9BZHJP*O z$S|?LL^=5{{&w4t67eU*=K;%Y4$ZM3%f)(yXYCbU%CkdLIQ77!HtAcnc3(X<>0y#p zJA3?h_T?89kWqzsxmh%APw8iQ$$pMzeLwG-^`ip_W!)C-^)zm+i*Q*>5lcA1>rLS& zdAxq2scLWFO5pJza(VRl4zuL@QP}3MDSY@`dBk7Q8GJ%NR>|Qy0Y3<4iT_{yJxRoE871bIj%*6fFRfV>~ zZ*wzfia6^S$HNpFoyT`$&@U=g>rAEd~VZ)kDLOpm1VEd}{zE=+oas+pu@N1NIV^sI(ZUaKXb0AIiEG9e;oh zRMkUL07U5_yMan`WIRyrBPjpT*C*Kl+ot8{A%c6}?GAyu2SazFN;i(Jp9_&?K%I=Z zyIu$^*Jpozvh`*3?RM?Ur}b?SYGeER#%Mz5ynnDC_R`#aYl>EG&2J(t=Kmb$ST=r| zz7Vwewe0aH0Q`3hIxy`i`mtd2WzRrp3=tY*dOc0Yz+H?%0rW+cfEbA};1@82->{V7SQvr)H9%UN>=@pJ1&**loCF+z<}-|7>}T-03P*1H z5RLd#oTj3PkvFlCS)ZoG${m^bK#8|ZcvC!rXODwy7)64HEO*QHJlYzJ_#_*Xr0c&YHMLuz6oO#;g%>5bLL}8!d zA-7O$bkw~h>-F6l3-5k;b7SN8tLnR|Y82)W5_=4@_T>X-`1rV^oW(e*0&;Eas^S2G zG<2wZ$bOvBOh6v$A)K*VX6O3Fjkf^hZ6I6xc;4pJtL7Bq-Cgs!coOYh} zyP4XmKiJtBDitZx7;I7AWV%|dqPl`ItxylTePm*H0{RB^H^*vLukS1PNu2Qs;~Hmv15Fd6G=AoXy$Z_=J+P(VF%QXGnZ>3r zdQjr1TO)%h3l_#wbc^F!hk`wFMc?0ZNat^sTaLYJYcPzFS4h^9v%DBi>@BQIKFh`f z8aA7UJ5RPp2Ri_c#?p$sO7d`qmSMijf zKJ_Ph5^KZ{VYr8J+J|Vif&!x(9o`}RHTeI^+WQUuF-(5BVzbl60a42Ij49}E!cAp_2=ciGhd(q5T-l7v%0 zJ>)7efqVR6^-ID5z^{*EmLQs>KI6H!@*BV$&byb!-Ik8;|BKJYhiV_9k;|?6h%dha zi3p-d=j|wxusuwGyWKpiUFB_(@3Jkuih7B8l}CO20%3&=Zn}pQfS*|7Gkr3^zlPLq zaKjc>yNgDE+;ZHhCVA~+p5El<52Lk6O1LNzLO_)S{MsLRnGs)nF&bJP8=92i+IA3MSd+1%aSO_&QfC))^0N? zM@MlY5rKD+8dp%_1`T#b?Hy|A!w zdSPLH_VlSs;e2T^TnbA|7lI2x2^FVrucAL2Otj(Ipd7p#mcylBIh+rcf<+`&3#d}( VF<%bnF~1O;3onE`DxtR&{0|*3#25eo literal 0 HcmV?d00001 diff --git a/calibcant/__init__.py b/calibcant/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/calibcant/analyze.py b/calibcant/analyze.py new file mode 100755 index 0000000..19a4d96 --- /dev/null +++ b/calibcant/analyze.py @@ -0,0 +1,238 @@ +#!/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 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 / = kb T photoSensitiviy**2 * (1e9nm/m)**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 \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 diff --git a/calibcant/bump.py b/calibcant/bump.py new file mode 100644 index 0000000..0473485 --- /dev/null +++ b/calibcant/bump.py @@ -0,0 +1,188 @@ +#!/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 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) + diff --git a/calibcant/bump_analyze.py b/calibcant/bump_analyze.py new file mode 100755 index 0000000..56b135d --- /dev/null +++ b/calibcant/bump_analyze.py @@ -0,0 +1,249 @@ +#!/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 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 \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' + ' should be whitespace-delimited, 2 column ASCII\n' + 'without a header line. e.g: "\\t\\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> 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() diff --git a/calibcant/calibrate_cantilever.py b/calibcant/calibrate_cantilever.py new file mode 100644 index 0000000..6721605 --- /dev/null +++ b/calibcant/calibrate_cantilever.py @@ -0,0 +1,911 @@ +#!/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 + # = , the variance for our AC signal. + # The variance from our fitted Lorentzian is the area under the Lorentzian + # = (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) + + diff --git a/calibcant/common.py b/calibcant/common.py new file mode 100644 index 0000000..a93adfd --- /dev/null +++ b/calibcant/common.py @@ -0,0 +1,38 @@ +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') diff --git a/calibcant/common.pyc b/calibcant/common.pyc new file mode 100644 index 0000000000000000000000000000000000000000..cfa016f1961f062b3065abdbd101c014ef9ec178 GIT binary patch literal 1618 zcmb_cQI8Wh5T4ylauA3Fr2>g3SE|&lgp!uW3RR`YNsUwmv^gO`g`&;In>Dv<@7Ye% zNI7_*|ET|=KZut;!Z+UIuBGZzapLjJ*faBuXS{#B9S#3ffBjU@`G~Rpfa$&iS>Yei z3DI*fF+Go1%#%QSf%a*Y(C;`fpj9t02dtl`eVD6v7&DiE{bqdsA_U8anC=aT51wDf zWNuIdON2RzDJRkqMINDI;>|B1^hou|jKIP5fN`A8U=!j)2nLA$-v$SC`nt2eLis^3 z!0CC~V=%8HqG{6_>qWb(^QCBOrQEqhMmC;U=#)@qvwb)9brCbSb!m+Us$CYwJ zf$`zAC~_^95=&jmeNWj!7eW{7L8eL zYB}!ur9d^pKtX7Z>9{^<#d*q_SqefQ5%&PPNo5*2*`6T!#MxqkB2EBGEp6`Q0_E_k zmcL5dDOsKrrfG~Gw>y6g4*r()U5E~o(yJhT0I~>KpbLS)Xz97bgDL~l9Vv`1)!N^L z{N&ls4}V-tXU9j`!^7k0)1zRxviUY_|LdB=wyVtclV>&pnnl+^T?#`e+~*)P97Q*x zQ8bK3@#uUl66iTyPVW%Y%|O;fy|aY&5j}cwi_UtqPx$5Td-Nuv&rYw=tAx(_w8z7Y zXdjb`Y2W8P{G(2bx6ML$a2&%u139wFOZ1SyV)G@QM{aGt6Hc~}Jvx4QAs9Oj{d6pC zgC2z)H;~!*cH_lX`HKCZD{H)qvEx=5&QgvPuGY%dhLUYLWKg{{N~fG^(3=w75x4-s zbKTEj;G4izU;Ogz3(tOk^9Q`?EjrCz*6Y(WVZgTt7{1^r43LtqqT5k=1zg5RLj0#k z*?jtVmhl%mv|z@lG6tA2`0N(QCAz@(l0k%DeD%(I=?jxb literal 0 HcmV?d00001 diff --git a/calibcant/config.py b/calibcant/config.py new file mode 100644 index 0000000..2374ba6 --- /dev/null +++ b/calibcant/config.py @@ -0,0 +1,25 @@ +"""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 diff --git a/calibcant/config.pyc b/calibcant/config.pyc new file mode 100644 index 0000000000000000000000000000000000000000..f8711bbc7523b982f8afe9f0d8fec10f83b76032 GIT binary patch literal 866 zcmZuw-EP`26gDZPh0^tB-7cqGB26OFuB&Qc1cpR`L5eotSSB%Gk;akjfZ%2J5PP$} zz;}LFc>EYJBR@m zQCvj$6vd|qmn?H;@!VqBVg)AYMiI(i(3NzOu$O!Z+5qMfj~X?ht-wYK>}V zriE_+O*Z8zbCq~z?vhGQCka!o7H%ZcDWBz%xq8mrj7AGOWA0Q)mx3oXkMe{{H=&aX zO=EW@?m)jVNKCe;Vf z*Ix~hDc9ybeg~48Nv_Qn2#G~Rh`t1SLo&`bg5;1*=_)YIpyfYxLlOp$16>3nlsVJq zV08O!()2^$Cb3RusIMW=?me|Sp^cG-KMHJ~F!&ylaWHK3fwxaLc)*W1lz z$BgNGCmnAxB69s{tjj= 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 \n' + '2008, W. Trevor King.\n' + '\n' + 'Apply various filters to data y(x).' + ' should be whitespace-delimited, 2 column ASCII\n' + 'without a header line.\n' + 'e.g: "\\t\\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) diff --git a/calibcant/psd_filter_analyze.py b/calibcant/psd_filter_analyze.py new file mode 100755 index 0000000..29c9641 --- /dev/null +++ b/calibcant/psd_filter_analyze.py @@ -0,0 +1,73 @@ +#!/usr/bin/python + +if __name__ == "__main__" : + + import sys + import filter, calibcant_vib_analyze + from optparse import OptionParser + + usage_string = ('%prog \n' + '2008, W. Trevor King.\n' + '\n' + 'Apply various filters to data y(x).' + ' should be whitespace-delimited, 2 column ASCII\n' + 'without a header line.\n' + 'e.g: "\\t\\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 diff --git a/calibcant/vib_analyze.py b/calibcant/vib_analyze.py new file mode 100755 index 0000000..837045c --- /dev/null +++ b/calibcant/vib_analyze.py @@ -0,0 +1,425 @@ +#!/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 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 + # = , the variance for our AC signal. + # The variance from our fitted Lorentzian is the area under the Lorentzian + # = (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 \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)' + ' should be whitespace-delimited, 2 column ASCII\n' + 'without a header line.\n' + 'e.g: "\\t\\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() -- 2.26.2