A bunch of changes in one commit, sorry.
Moved to fledgling splittable_kwargs system to make default argument
maintenance easier. I expect the splittable_kwargs system still has
some growing to do, but it's already better than the old system.
Merged BE database from the calibcant subdir into the main BE database.
It was a mistake to create the database there in the first place.
--- /dev/null
+z_piezo_calib output confirmed with
+
+$ cd /home/wking/research/drexel/gyang/data/z_piezo_calib
+$ ipython --pylab
+In [1]: import z_piezo_calib
+In [2]: z_piezo_calib.analyze_calibration_data_from_tweaks(bump_tweaks="20081215_tweak.bump", vib_tweaks="20081215_tweak.vib", T_tweaks="20081215_tweak.temp", plotVerbose=False)
+a**2 : 0.000190122 +/- 1.86869e-06 (0.00982893)
+T : 295.15 +/- 5.68434e-14 (1.92592e-16)
+1/Vp**2 : 61416.8 +/- 17792.9 (0.289708)
+Out[2]:
+(0.0475823000538,
+ 0.0137928827301,
+ 0.00982893194938,
+ 1.92591627514e-16,
+ 0.289707547407)
+
+In [3]: z_piezo_calib.analyze_calibration_data_from_tweaks(bump_tweaks="20081215_tweak_second.bump", vib_tweaks="20081215_tweak_second.vib", T_tweaks="20081215_tweak_second.temp", plotVerbose=False)
+a**2 : 0.000217743 +/- 1.08297e-06 (0.00497362)
+T : 295.15 +/- 5.68434e-14 (1.92592e-16)
+1/Vp**2 : 38507.9 +/- 13142.7 (0.341299)
+Out[3]:
+(0.034168120631,
+ 0.0116627940459,
+ 0.00497362408195,
+ 1.92591627514e-16,
+ 0.341299306648)
+
+However, plotVerbose is broken (missing figure definition), so I'm
+gonna go fix that now.
+
--- /dev/null
+Content-type: text/plain
+
+
+Date: Fri, 19 Dec 2008 06:17:07 +0000
+
+
+From: W. Trevor King <wking@drexel.edu>
+
+
+In-reply-to: daf08023-7758-482e-a481-800872e84312
+
--- /dev/null
+I just reran Marissa's calibration. My original output (from the old
+gnuplot-fitted ~/.python/calibrate_cantilever.py) had produced
+
+ $ cat 20081215/20081215101516_analysis_text
+ k (N/m) : 0.0471212 +/- 0.0132663
+ a**2 (V/nm)**2 : 0.000191658 +/- 1.36722e-06
+ T (K) : 295.15 +/- 5.68434e-14
+ 1/Vp**2 (1/V)**2 : 60334.2 +/- 16980.8
+ $ cat 20081215/20081215124856_analysis_text
+ k (N/m) : 0.0346675 +/- 0.0118672
+ a**2 (V/nm)**2 : 0.000217743 +/- 1.08297e-06
+ T (K) : 295.15 +/- 5.68434e-14
+ 1/Vp**2 (1/V)**2 : 39070.7 +/- 13373.1
+
+Using all the defaults (and only the second take for the first
+calibration). However, once I'd set up the tweakfiles and ran
+20081215_calibrate.sh, I got
+
+ $ cat 20081215/20081215101516_analysis_text
+ Variable (units) : mean +/- std. dev. (relative error)
+ Cantilever k (N/m) : 0.0656644 +/- 0.00107743 (0.016408)
+ photoSensitivity**2 (V/nm)**2 : 0.000190122 +/- 1.86869e-06 (0.00982893)
+ T (K) : 295.15 +/- 5.68434e-14 (1.92592e-16)
+ 1/Vphoto**2 (1/V)**2 : 84756.3 +/- 1113.56 (0.0131383)
+ $ cat 20081215/20081215124856_analysis_text
+ Variable (units) : mean +/- std. dev. (relative error)
+ Cantilever k (N/m) : 0.0494809 +/- 0.000855265 (0.0172848)
+ photoSensitivity**2 (V/nm)**2 : 0.000217743 +/- 1.08297e-06 (0.00497362)
+ T (K) : 295.15 +/- 5.68434e-14 (1.92592e-16)
+ 1/Vphoto**2 (1/V)**2 : 55765.6 +/- 923.129 (0.0165537)
+
+Comparing the two, we see:
+ 101516:
+ k 0.0471212 -> 0.0656644 (with a 92% drop in error)
+ photoSensitivity 0.000191658 (little change)
+ T 295.15 (no change)
+ 1/Vphoto**2 60334.2 -> 84756.3 (with a 94% drop in error)
+ 124856:
+ k 0.0346675 -> 0.0494809 (with a 93% drop in error)
+ photoSensitivity 0.000217743 (no change)
+ T 295.15 (no change)
+ 1/Vphoto**2 39070.7 -> 55765.6 (with a 93% drop in error)
+
+I don't have hard evidence yet, but I am confident this discrepancy is
+due to poor fitting with the older calibration file.
+
--- /dev/null
+Content-type: text/plain
+
+
+Date: Fri, 19 Dec 2008 05:25:39 +0000
+
+
+From: W. Trevor King <wking@drexel.edu>
+
--- /dev/null
+Yikes, z_piezo_calib uses the actual variance (no fitting)! No wonder
+it was so ugly.
--- /dev/null
+Content-type: text/plain
+
+
+Date: Fri, 19 Dec 2008 06:22:45 +0000
+
+
+From: W. Trevor King <wking@drexel.edu>
+
+
+In-reply-to: 651ebc04-cc59-4996-b4f3-33bfd9c88c23
+
--- /dev/null
+Oops, I was running ~/.python/z_piezo_calib.calibrate_cantilever()
+during the experiment, not ~/.python/calibrate_cantilever.something().
+I was wondering why it was so difficult to get
+~./python/calibrate_cantilever.py working. I ended up getting
+
+ Cantilever k : 0.589927 +/- 0.790152 (1.33941)
+ photoSensitivity**2 : 0.000190122 +/- 1.86869e-06 (0.00982893)
+ T : 295.15 +/- 0 (0)
+ 1/Vphoto**2 : 761447 +/- 1.01986e+06 (1.33937)
+
+out of calibrate_cantilever.py for the first cantilever, but that is
+pretty rediculous... Maybe I'll come back and figure out what's going on
+after I check out the z_piezo_calib version...
--- /dev/null
+Content-type: text/plain
+
+
+Date: Fri, 19 Dec 2008 06:09:39 +0000
+
+
+From: W. Trevor King <wking@drexel.edu>
+
+
+In-reply-to: 651ebc04-cc59-4996-b4f3-33bfd9c88c23
+
--- /dev/null
+creator: W. Trevor King <wking@drexel.edu>
+
+
+reporter: W. Trevor King <wking@drexel.edu>
+
+
+severity: minor
+
+
+status: fixed
+
+
+summary: Suprising difference between old calibration and new calibration
+
+
+time: Fri, 19 Dec 2008 05:08:51 +0000
+
--- /dev/null
+Some files use config.TEXT_VERBOSE boolean, others use textVerboseFile option.
+Pick one method and go with it.
--- /dev/null
+Content-type: text/plain
+
+
+Date: Sun, 21 Dec 2008 03:07:58 +0000
+
+
+From: W. Trevor King <wking@drexel.edu>
+
--- /dev/null
+creator: W. Trevor King <wking@drexel.edu>
+
+
+reporter: W. Trevor King <wking@drexel.edu>
+
+
+severity: minor
+
+
+status: open
+
+
+summary: unify TEXT_VERBOSE
+
+
+time: Sun, 21 Dec 2008 03:07:08 +0000
+
-
-
-
-rcs_name=git
-
-
+rcs_name: git
+++ /dev/null
-
-
-
-rcs_name=git
-
-
-
+++ /dev/null
-Bugs Everywhere Tree 1 0
import config # config module for the calibcant package
import data_logger
import linfit
+from splittable_kwargs import splittableKwargsFunction, \
+ make_splittable_kwargs_function
def C_to_K(celsius) :
"Convert Celsius -> Kelvin."
return celsius + 273.15
+def K_to_K(kelvin) :
+ "Convert Kelvin -> Kelvin."
+ return kelvin
+
+@splittableKwargsFunction()
def T_analyze(T, convert_to_K=C_to_K) :
"""
Not much to do here, just convert to Kelvin.
T = [convert_to_K(T)]
return T
+@splittableKwargsFunction()
def T_save(T, log_dir=None) :
"""
Save either a single T (if you are crazy :p),
dl = data_logger.data_load()
return dl.read_binary(datafile)
-def T_plot(T, plotVerbose) :
+@splittableKwargsFunction()
+def T_plot(T, plotVerbose=False) :
"""
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) :
+@splittableKwargsFunction((T_analyze, 'T', 'convert_to_K'),
+ (T_plot, 'T'))
+def T_load_analyze_tweaked(tweak_file, convert_to_K=C_to_K, textVerboseFile=None, **kwargs) :
"Load all the T array files from a tweak file and return a single array"
+ T_analyze_kwargs,T_plot_kwargs = \
+ T_load_analyze_tweaked._splitargs(T_load_analyze_tweaked, kwargs)
Ts = []
for line in file(tweak_file, 'r') :
parsed = line.split()
# read the data
data = T_load(path)
Ts.extend(data)
+ T_analyze(Ts, convert_to_K=K_to_K)
return numpy.array(Ts, dtype=numpy.float)
# commandline interface functions
"""
import numpy
+from splittable_kwargs import splittableKwargsFunction, \
+ make_splittable_kwargs_function
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) :
+#@splittableKwargsFunction((calib_plot, 'bumps', 'Ts', 'vibs'))
+# Some of the child functions aren't yet defined, so postpone
+# make-splittable until later in the module.
+def calib_analyze(bumps, Ts, vibs, **kwargs) :
"""
Analyze data from get_calibration_data()
return (k, k_s, !!!a2_r, T_r, one_o_Vp2_r)
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.
"""
+ calib_plot_kwargs, = calib_analyze._splitargs(calib_analyze, kwargs)
photoSensitivity2 = bumps**2
one_o_Vphoto2 = 1/vibs
k_s = k*(photoSensitivity2_r**2 + T_r**2 + one_o_Vphoto2_r**2)**0.5
+ calib_plot(bumps, Ts, vibs, **calib_plot_kwargs)
+
return (k, k_s,
photoSensitivity2_m, photoSensitivity2_s,
T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s)
+@splittableKwargsFunction()
def string_errors(k_m, k_s,
photoSensitivity2_m, photoSensitivity2_s,
T_m, T_s,
% (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()
+@splittableKwargsFunction()
+def calib_save(bumps, Ts, vibs, log_dir=None) :
+ """
+ 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)
+ ))
+
+@splittableKwargsFunction()
+def calib_plot(bumps, Ts, vibs, plotVerbose=False) :
+ if plotVerbose or config.PYLAB_VERBOSE :
+ common.import_pylab()
+ common._pylab.figure(config.BASE_FIGNUM+4)
+ common._pylab.subplot(311)
+ common._pylab.plot(bumps, 'g.')
+ common._pylab.title('Photodiode sensitivity (V/nm)')
+ common._pylab.subplot(312)
+ common._pylab.plot(Ts, 'r.')
+ common._pylab.title('Temperature (K)')
+ common._pylab.subplot(313)
+ common._pylab.plot(vibs, 'b.')
+ common._pylab.title('Thermal deflection variance (Volts**2)')
+ common._flush_plot()
+
+make_splittable_kwargs_function(calib_analyze,
+ (calib_plot, 'bumps', 'Ts', 'vibs'))
+
+@splittableKwargsFunction((calib_analyze, 'bumps', 'Ts', 'vibs'))
+def calib_load_analyze_tweaks(bump_tweaks, vib_tweaks, T_tweaks=None) :
+ raise NotImplementedError
+ a = read_tweaked_bumps(bump_tweaks)
+ vib = V_photo_variance_from_file(vib_tweaks)
+ if T_tweaks == None:
+ pass
+ return analyze_calibration_data(a, T, vib, log_dir=log_dir)
# commandline interface functions
import scipy.io, sys
zpGain Vzp_out / Vzp
zpSensitivity Zp / Vzp
photoSensitivity Vphoto / Zcant
+
+photoSensitivity is measured by bumping the cantilever against the
+surface, where Zp = Zcant (see calibrate.bump_aquire()). The measured
+slope Vphoto/Vout is converted to photoSensitivity with bump_analyze().
"""
import numpy
import data_logger
import z_piezo_utils
import linfit
+from splittable_kwargs import splittableKwargsFunction, \
+ make_splittable_kwargs_function
+#@splittableKwargsFunction((bump_plot, 'data'))
+# Some of the child functions aren't yet defined, so postpone
+# make-splittable until later in the module.
def bump_analyze(data, zpGain=config.zpGain,
zpSensitivity=config.zpSensitivity,
Vzp_out2V=config.Vzp_out2V,
Vphoto_in2V=config.Vphoto_in2V,
- textVerboseFile=None, plotVerbose=False) :
+ **kwargs) :
"""
Return the slope of the bump ;).
Inputs:
and THEN converted, so we're assuming that both conversions are LINEAR.
if they aren't, rewrite to convert before the regression.
"""
+ bump_plot_kwargs, = bump_analyze._splitargs(bump_analyze, kwargs)
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
+ plotVerbose=False)
+ 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
+ bump_plot(data, **bump_plot_kwargs)
return Vphoto2Vzp_out * Vzp_out2Zcant
-def bump_save(data, log_dir) :
+@splittableKwargsFunction()
+def bump_save(data, log_dir=None) :
"Save the dictionary data, using data_logger.data_log()"
if log_dir != None :
log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
data = dl.read_dict_of_arrays(datafile)
return data
-def bump_plot(data, plotVerbose) :
+@splittableKwargsFunction()
+def bump_plot(data, plotVerbose=False) :
"Plot the bump (Vphoto vs Vzp) if plotVerbose or PYLAB_VERBOSE == True"
if plotVerbose or config.PYLAB_VERBOSE :
common._import_pylab()
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) :
+make_splittable_kwargs_function(bump_analyze, (bump_plot, 'data'))
+
+@splittableKwargsFunction((bump_analyze, 'data'))
+def bump_load_analyze_tweaked(tweak_file, **kwargs):
"Load the output file of tweak_calib_bump.sh, return an array of slopes"
+ bump_analyze_kwargs, = \
+ bump_load_analyze_tweaked._splitargs(bump_load_analyze_tweaked, kwargs)
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:])
+ if config.TEXT_VERBOSE :
+ print "Reading data from %s with ranges %s" % (path, parsed[1:])
# read the data
full_data = bump_load(path)
if len(parsed) == 1 :
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)
+ pSi = bump_analyze(data, **bump_analyze_kwargs)
photoSensitivity.append(pSi)
- bump_plot(data, plotVerbose)
return numpy.array(photoSensitivity, dtype=numpy.float)
# commandline interface functions
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_INTERACTIVE = False
config.PYLAB_VERBOSE = options.pylab
# data not very bilinear, so don't cut anything off.
print >> sys.stderr, "use everything"
- photoSensitivity = bump_analyze(data, textVerboseFile=vfile)
+ photoSensitivity = bump_analyze(data)
print >> ofile, photoSensitivity
else : # tweak file mode
- slopes = bump_load_analyze_tweaked(ifilename, textVerboseFile=vfile)
+ slopes = bump_load_analyze_tweaked(ifilename)
if options.comma_out :
sep = ','
else :
--- /dev/null
+#!/usr/bin/python
+
+"""
+Aquire and analyze cantilever calibration data.
+
+W. Trevor King Dec. 2007-Jan. 2008
+
+GPL BOILERPLATE
+
+
+The relevent physical quantities are :
+ Vzp_out Output z-piezo voltage (what we generate)
+ Vzp Applied z-piezo voltage (after external ZPGAIN)
+ Zp The z-piezo position
+ Zcant The cantilever vertical deflection
+ Vphoto The photodiode vertical deflection voltage (what we measure)
+ Fcant The force on the cantilever
+ T The temperature of the cantilever and surrounding solution
+ (another thing we measure or guess)
+ 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 bump_aquire() and the bump_analyze
+submodule).
+
+k_cant is measured by watching the cantilever vibrate in free solution
+(see the vib_aquire() and the vib_analyze submodule). 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 with the surface bumps. We can either
+measure T using an external function (see temperature.py), or just
+estimate it (see T_aquire() and the T_analyze submodule). Guessing
+room temp ~22 deg C is actually fairly reasonable. Assuming the
+actual fluid temperature is within +/- 5 deg, the error in the spring
+constant k_cant is within 5/273.15 ~= 2%. A time series of Vphoto
+while we're far from the surface and not changing Vzp_out will give us
+the average variance <Vphoto**2>.
+
+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().
+
+We also define the two positioning functions:
+ move_just_onto_surface() and move_far_from_surface()
+which make automating the calibration procedure more straightforward.
+"""
+
+import numpy
+import time
+import z_piezo_utils
+from splittable_kwargs import splittableKwargsFunction, \
+ make_splittable_kwargs_function
+
+import common
+import config
+import bump_analyze
+import T_analyze
+import vib_analyze
+import analyze
+
+# bump family
+
+@splittableKwargsFunction()
+def bump_aquire(zpiezo, push_depth=200, npoints=1024, freq=100e3) :
+ """
+ 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
+ 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 config.TEXT_VERBOSE :
+ print "Bump %g nm" % push_depth
+ data = zpiezo.ramp(out, freq)
+ return data
+
+@splittableKwargsFunction(bump_aquire,
+ (bump_analyze.bump_save, 'data'),
+ (bump_analyze.bump_analyze, 'data'))
+def bump(**kwargs):
+ """
+ Wrapper around bump_aquire(), bump_save(), bump_analyze()
+ """
+ bump_aquire_kwargs,bump_save_kwargs,bump_analyze_kwargs = \
+ bump._splitargs(bump, kwargs)
+ data = bump_aquire(**bump_aquire_kwargs)
+ bump_analyze.bump_save(data, **bump_save_kwargs)
+ photoSensitivity = bump_analyze.bump_analyze(data, **bump_analyze_kwargs)
+ return photoSensitivity
+
+# T family.
+# Fairly stubby, since a one shot Temp measurement is a common thing.
+# We just wrap that to provide a consistent interface.
+
+@splittableKwargsFunction()
+def T_aquire(get_T=None) :
+ """
+ Measure the current temperature of the sample,
+ or, if get_T == None, fake it by returning config.DEFAULT_TEMP
+ """
+ if get_T == None :
+ if config.TEXT_VERBOSE :
+ print "Fake temperature %g" % config.DEFAULT_TEMP
+ return config.DEFAULT_TEMP
+ else :
+ if config.TEXT_VERBOSE :
+ print "Measure temperature"
+ return get_T()
+
+@splittableKwargsFunction(T_aquire,
+ (T_analyze.T_save, 'T'),
+ (T_analyze.T_analyze, 'T'))
+def T(**kwargs):
+ """
+ Wrapper around T_aquire(), T_save(), T_analyze(), T_plot()
+ """
+ T_aquire_kwargs,T_save_kwargs,T_analyze_kwargs = \
+ T._splitargs(T, kwargs)
+ T_raw = T_aquire(**T_aquire_kwargs)
+ T_analyze.T_save(T_raw, **T_save_kwargs)
+ T_ret = T_analyze.T_analyze(T_raw, **T_analyze_kwargs)
+ return T_ret
+
+# vib family
+
+@splittableKwargsFunction()
+def vib_aquire(zpiezo, time=1, freq=50e3) :
+ """
+ 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 config.TEXT_VERBOSE :
+ print "get %g seconds of data" % time
+ data = zpiezo.ramp(out, freq)
+ data['sample frequency Hz'] = freq
+ return data
+
+@splittableKwargsFunction(vib_aquire,
+ (vib_analyze.vib_save, 'data'),
+ (vib_analyze.vib_analyze, 'deflection_bits', 'freq'))
+def vib(**kwargs) :
+ """
+ Wrapper around vib_aquire(), vib_save(), vib_analyze()
+ """
+ vib_aquire_kwargs,vib_save_kwargs,vib_analyze_kwargs = \
+ vib._splitargs(vib, kwargs)
+ data = vib_aquire(freq=freq, **vib_aquire_kwargs)
+ vib_analyze.vib_save(data, **vib_save_kwargs)
+ freq = data['sample frequency Hz']
+ Vphoto_var = vib_analyze.vib_analyze(deflection_bits=data, freq=freq,
+ **vib_analyze_kwargs)
+ return Vphoto_var
+
+# A few positioning functions, so we can run bump_aquire() and vib_aquire()
+# with proper spacing relative to the surface.
+
+@splittableKwargsFunction()
+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 config.TEXT_VERBOSE :
+ print "moving just onto surface"
+ # Zero the piezo
+ if config.TEXT_VERBOSE :
+ print "zero the z piezo output"
+ zpiezo.jumpToPos(zpiezo.pos_nm2out(0))
+ # See if we're near the surface already
+ if config.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 config.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 config.TEXT_VERBOSE :
+ print "Single stepping approach"
+ while cd < sp :
+ if config.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 config.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 config.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 config.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 config.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 config.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 config.TEXT_VERBOSE :
+ print "We're %g nm into the surface" % Depth_nm
+
+@splittableKwargsFunction()
+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
+
+@splittableKwargsFunction((move_just_onto_surface, 'stepper', 'zpiezo'),
+ (bump, 'zpiezo', 'freq', 'log_dir', 'Vphoto_in2V'),
+ (move_far_from_surface, 'stepper'),
+ (T, 'log_dir'),
+ (vib, 'zpiezo', 'log_dir', 'Vphoto_in2V'),
+ (analyze.calib_save, 'bumps','Ts','vibs','log_dir'))
+def calib_aquire(stepper, zpiezo, num_bumps=10, num_Ts=10, num_vibs=20,
+ bump_freq=100e3,
+ log_dir=config.LOG_DATA, Vphoto_in2V=config.Vphoto_in2V,
+ **kwargs):
+ """
+ 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
+ """
+ move_just_onto_surface_kwargs,bump_kwargs,move_far_from_surface_kwargs, \
+ T_kwargs,vib_kwargs,calib_save_kwargs = \
+ calib_aquire._splitargs(calib_aquire, kwargs)
+ # get bumps
+ move_just_onto_surface(stepper, zpiezo, **move_just_onto_surface_kwargs)
+ bumps=zeros((num_bumps,))
+ for i in range(num_bumps) :
+ bumps[i] = bump(zpiezo, freq=bump_freq, log_dir=log_dir,
+ Vphot_in2V=Vphoto_in2V, **bump_kwargs)
+ if config.TEXT_VERBOSE :
+ print bumps
+
+ move_far_from_surface(stepper, **move_far_from_surface_kwargs)
+
+ # get Ts
+ Ts=zeros((num_Ts,), dtype=float)
+ for i in range(num_Ts) :
+ Ts[i] = T(**T_kwargs)
+ time.sleep(1) # wait a bit to get an independent temperature measure
+ print Ts
+
+ # get vibs
+ vibs=zeros((num_vibs,))
+ for i in range(num_vibs) :
+ vibs[i] = vib(zpiezo, log_dir=log_dir, Vphoto_in2V=Vphoto_in2V,
+ **vib_kwargs)
+ print vibs
+
+ analyze.calib_save(bumps, Ts, vibs, log_dir, **calib_save_kwargs)
+
+ return (bumps, Ts, vibs)
+
+
+@splittableKwargsFunction( \
+ (calib_aquire, 'log_dir'),
+ (analyze.calib_analyze, 'bumps','Ts','vibs'))
+def calib(log_dir=None, **kwargs) :
+ """
+ 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
+ """
+ calib_aquire_kwargs,calib_analyze_kwargs = \
+ calib._splitargs(calib, kwargs)
+ a, T, vib = calib_aquire(**calib_aquire_kwargs)
+ k,k_s,ps2_m, ps2_s,T_m,T_s,one_o_Vp2_m,one_o_Vp2_s = \
+ analyze.calib_analyze(a, T, vib, log_dir=log_dir,
+ **calib_analyze_kwargs)
+ analyze.calib_save_analysis(k, k_s, ps2_m, ps2_s, T_m, T_s,
+ one_o_Vp2_m, one_o_Vp2_s, log_dir)
+ return (k, k_s)
+
+
+
+++ /dev/null
-#!/usr/bin/python
-
-"""
-Aquire and analyze cantilever calibration data.
-
-W. Trevor King Dec. 2007-Jan. 2008
-
-GPL BOILERPLATE
-
-
-The relevent physical quantities are :
- Vzp_out Output z-piezo voltage (what we generate)
- Vzp Applied z-piezo voltage (after external ZPGAIN)
- Zp The z-piezo position
- Zcant The cantilever vertical deflection
- Vphoto The photodiode vertical deflection voltage (what we measure)
- Fcant The force on the cantilever
- T The temperature of the cantilever and surrounding solution
- (another thing we measure)
- k_b Boltzmann's constant
-
-Which are related by the parameters :
- zpGain Vzp_out / Vzp
- zpSensitivity Zp / Vzp
- photoSensitivity Vphoto / Zcant
- k_cant Fcant / Zcant
-
-Cantilever calibration assumes a pre-calibrated z-piezo (zpSensitivity) and
-a amplifier (zpGain).
-In our lab, the z-piezo is calibrated by imaging a calibration sample,
-which has features with well defined sizes, and the gain is set with a knob
-on the Nanoscope.
-
-photoSensitivity is measured by bumping the cantilever against the surface,
-where Zp = Zcant (see the bump_*() family of functions)
-The measured slope Vphoto/Vout is converted to photoSensitivity via
-Vphoto/Vzp_out * Vzp_out/Vzp * Vzp/Zp * Zp/Zcant = Vphoto/Zcant
- (measured) (1/zpGain) (1/zpSensitivity) (1) (photoSensitivity)
-
-k_cant is measured by watching the cantilever vibrate in free solution
-(see the vib_*() family of functions)
-The average energy of the cantilever in the vertical direction is given
-by the equipartition theorem.
- 1/2 k_b T = 1/2 k_cant Zcant**2
- so k_cant = k_b T / Zcant**2
- but Zcant = Vphoto / photoSensitivity
- so k_cant = k_b T * photoSensitivty**2 / Vphoto**2
-We measured photoSensitivity above.
-We can either measure T using an external function (see temperature.py),
-or just estimate it (see the T_*() family of functions).
- (room temp ~22 deg C, accurate to +/- 5 deg, so 5/273.15 ~= 2%)
-Finally, a time series of Vphoto while we're far from the surface and not
-changing Vzp_out will give us Vphoto.
-
-We do all these measurements a few times to estimate statistical errors.
-
-The functions are layed out in the families:
- bump_*(), vib_*(), T_*(), and calib_*()
-where calib_{save|load|analyze}() deal with derived data, not real-world data.
-For each family, * can be any of :
- aquire get real-world data
- save store real-world data to disk
- load get real-world data from disk
- analyze interperate the real-world data.
- plot show a nice graphic to convince people we're working :p
- load_analyze_tweaked
- read a file with a list of paths to previously saved real world data
- load each file using *_load(), analyze using *_analyze(), and
- optionally plot using *_plot().
- Intended for re-processing old data.
-A family name without any _* extension (e.g. bump()),
- runs *_aquire(), *_save(), *_analyze(), *_plot().
-
-Also defines the two positioning functions:
- move_just_onto_surface() and move_far_from_surface()
-"""
-
-import numpy
-import time
-import data_logger
-import z_piezo_utils
-import FFT_tools
-import linfit
-import GnuplotBiDir # used for fitting lorentzian
-
-kb = 1.3806504e-23 # Boltzmann's constant in J/K
-DEFAULT_TEMP = 22 # assume 22 deg C
-
-LOG_DATA = True # quietly grab all real-world data and log to LOG_DIR
-LOG_DIR = '/home/wking/rsrch/data/calibrate_cantilever'
-
-GNUFIT_DATA_BASE='./calibrate_cantilever_fitdata'
-
-TEXT_VERBOSE = True # for debugging
-GNUPLOT_VERBOSE = True # turn on fit check plotting
-PYLAB_VERBOSE = True # turn on plotting
-PYLAB_INTERACTIVE = True # select between draw() and show() for flushing plots
-BASE_FIGNUM = 20 # to avoid writing to already existing figures
-
-# handle extra verbose input modules, only imported if we need them
-_flush_plot = None
-_final_flush_plot = None
-_pylab = None
-def _import_pylab() :
- "Import pylab plotting functions for when we need to plot"
- global _pylab
- global _flush_plot
- global _final_flush_plot
- if _pylab == None :
- import pylab as _pylab
- if _flush_plot == None :
- if PYLAB_INTERACTIVE :
- _flush_plot = _pylab.draw
- else :
- def _flush_plot () : pass
- if _final_flush_plot == None :
- if PYLAB_INTERACTIVE :
- _final_flush_plot = _pylab.draw
- else :
- _final_flush_plot = _pylab.show
-
-
-# HACK
-# make sure you make a system note (add_system_note) if you change these
-# in case you don't have access to a z_piezo.z_piezo for conversion functions
-_usual_zpGain=20
-_usual_zpSensitivity=7.41
-_usual_zeroVphoto_bits = 2**15
-_usual_zeroVzp_bits = 2**15
-def _usual_Vphoto_in2V(inp) :
- return (float(inp)-float(_usual_zeroVphoto_bits))*(10.0/2**15)
-def _usual_Vzp_out2V(inp) :
- return (float(inp)-float(_usual_zeroVzp_bits))*(10.0/2**15)
-
-def C_to_K(celsius) :
- "Convert Celsius -> Kelvin."
- return celsius + 273.15
-
-# bump family
-
-def bump_aquire(zpiezo, push_depth, npoints, freq) :
- """
- Ramps closer push_depth and returns to the original position.
- Inputs:
- zpiezo an opened zpiezo.zpiezo instance
- push_depth distance to approach, in nm
- npoints number points during the approach and during the retreat
- freq rate at which data is aquired
- log_dir directory to log data to (see data_logger.py).
- None to turn off logging (see also the global LOG_DATA).
- Returns the aquired ramp data dictionary, with data in DAC/ADC bits.
- """
- # generate the bump output
- start_pos = zpiezo.curPos()
- pos_dist = zpiezo.pos_nm2out(push_depth) - zpiezo.pos_nm2out(0)
- close_pos = start_pos + pos_dist
- appr = linspace(start_pos, close_pos, npoints)
- retr = linspace(close_pos, start_pos, npoints)
- out = concatenate((appr, retr))
- # run the bump, and measure deflection
- if TEXT_VERBOSE :
- print "Bump %g nm" % push_depth
- data = zpiezo.ramp(out, freq)
- # default saving, so we have a log in-case the operator is lazy ;)
- if LOG_DATA == True :
- log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
- log_name="bump_surface")
- log.write_dict_of_arrays(data)
- return data
-
-def bump_save(data, log_dir) :
- "Save the dictionary data, using data_logger.data_log()"
- if log_dir != None :
- log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
- log_name="bump")
- log.write_dict_of_arrays(data)
-
-def bump_load(datafile) :
- "Load the dictionary data, using data_logger.date_load()"
- dl = data_logger.data_load()
- data = dl.read_dict_of_arrays(path)
- return data
-
-def bump_analyze(data, zpGain=_usual_zpGain,
- zpSensitivity=_usual_zpSensitivity,
- Vzp_out2V=_usual_Vzp_out2V,
- Vphoto_in2V=_usual_Vphoto_in2V,
- plotVerbose=False) :
- """
- Return the slope of the bump ;).
- Inputs:
- data dictinary of data in DAC/ADC bits
- pos_out2V function that converts output DAC bits to Volts
- def_in2V funtion that converts input ADC bits to Volts
- zpGain zpiezo applied voltage per output Volt
- zpSensitivity nm zpiezo response per applied Volt
- Returns:
- photoSensitivity (Vphoto/Zcant) in Volts/nm
- Checks for strong correlation (r-value) and low randomness chance (p-value)
-
- With the current implementation, the data is regressed in DAC/ADC bits
- and THEN converted, so we're assuming that both conversions are LINEAR.
- if they aren't, rewrite to convert before the regression.
- """
- scale_Vzp_bits2V = Vzp_out2V(1) - Vzp_out2V(0)
- scale_Vphoto_bits2V = Vzp_in2V(1) - Vphoto_in2V(0)
- Vphoto2Vzp_out_bit, intercept = \
- linfit.linregress(x=data["Z piezo output"],
- y=data["Deflection input"],
- plotVerbose=plotVerbose)
- Vphoto2Vzp_out = Vphoto2Vzp_out * scale_Vphoto_bits2V / scale_Vzp_bits2V
-
- # 1 / (Vzp/Vzp_out * Zp/Vzp * Zcant/Zp )
- Vzp_out2Zcant = 1.0/ (zpiezo_gain * zpiezo_sensitivity) # * 1
- return Vphoto2Vzp_out * Vzp_out2Zcant
-
-def bump_plot(data, plotVerbose) :
- "Plot the bump (Vphoto vs Vzp) if plotVerbose or PYLAB_VERBOSE == True"
- if plotVerbose or PYLAB_VERBOSE :
- _import_pylab()
- _pylab.figure(BASE_FIGNUM)
- _pylab.plot(data["Z piezo output"], data["Deflection input"],
- '.', label='bump')
- _pylab.title("bump surface")
- _pylab.legend(loc='upper left')
- _flush_plot()
-
-def bump(zpiezo, push_depth, npoints=1024, freq=100e3,
- log_dir=None,
- plotVerbose=False) :
- """
- Wrapper around bump_aquire(), bump_save(), bump_analyze(), bump_plot()
- """
- data = bump_aquire(zpiezo, push_depth, npoints, freq)
- bump_save(data, log_dir)
- photoSensitivity = bump_analyze(data, zpiezo.gain, zpiezo.sensitivity,
- zpiezo.pos_out2V, zpiezo.def_in2V)
- bump_plot(data, plotVerbose)
- return photoSensitivity
-
-def bump_load_analyze_tweaked(tweak_file, zpGain=_usual_zpGain,
- zpSensitivity=_usual_zpSensitivity,
- Vzp_out2V=_usual_Vzp_out2V,
- Vphoto_in2V=_usual_Vphoto_in2V,
- plotVerbose=False) :
- "Load the output file of tweak_calib_bump.sh, return an array of slopes"
- photoSensitivity = []
- for line in file(tweak_file, 'r') :
- parsed = line.split()
- path = parsed[0].split('\n')[0]
- # read the data
- full_data = bump_load(path)
- if len(parsed) == 1 :
- data = full_data # use whole bump
- else :
- # use the listed sections
- zp = []
- df = []
- for rng in parsed[1:] :
- p = rng.split(':')
- starti = int(p[0])
- stopi = int(p[1])
- zp.extend(full_data['Z piezo output'][starti:stopi])
- df.extend(full_data['Deflection input'][starti:stopi])
- data = {'Z piezo output': array(zp),
- 'Deflection input':array(df)}
- pSi = bump_analyze(data, zpGain, zpSensitivity,
- Vzp_out2V, Vphoto_in2V, plotVerbose)
- photoSensitivity.append(pSi)
- bump_plot(data, plotVervose)
- return array(photoSensitivity, dtype=numpy.float)
-
-# T family.
-# Fairly stubby, since a one shot Temp measurement is a common thing.
-# We just wrap that to provide a consistent interface.
-
-def T_aquire(get_T=None) :
- """
- Measure the current temperature of the sample,
- or, if get_T == None, fake it by returning DEFAULT_TEMP
- """
- if get_T == None :
- if TEXT_VERBOSE :
- print "Fake temperature %g" % DEFAULT_TEMP
- return DEFAULT_TEMP
- else :
- if TEXT_VERBOSE :
- print "Measure temperature"
- return get_T()
-
-def T_save(T, log_dir=None) :
- """
- Save either a single T (if you are crazy :p),
- or an array of them (more reasonable)
- """
- T = numpy.array(T, dtype=numpy.float)
- if log_dir != None :
- log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
- log_name="T_float")
- log.write_binary(T.tostring())
- if LOG_DATA != None :
- log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
- log_name="T_float")
- log.write_binary(T.tostring())
-
-def T_load(datafile=None) :
- """
- Load the saved T array (possibly of length 1), and return it.
- If datafile == None, return an array of one DEFAULT_TEMP instead.
- """
- if datafile == None :
- return numpy.array([DEFAULT_TEMP], dtype=numpy.float)
- else :
- return fromfile(file=cleanpath, dtype=numpy.float)
-
-def T_analyze(T, convert_to_K=C_to_K) :
- """
- Not much to do here, just convert to Kelvin.
- Uses C_to_K (defined above) by default.
- """
- try : # if T is an array, convert element by element
- for i in range(len(T)) :
- T[i] = convert_to_K(T[i])
- except TypeError :
- T = convert_to_K(T)
- return T
-
-def T_plot(T, plotVerbose) :
- """
- Umm, just print the temperature?
- """
- if plotVerbose or PYLAB_VERBOSE or TEXT_VERBOSE :
- print "Temperature ", T
-
-def T(get_T=None, convert_to_K=C_to_K, log_dir=None, plotVerbose=None) :
- """
- Wrapper around T_aquire(), T_save(), T_analyze(), T_plot()
- """
- T_raw = T_aquire(get_T)
- T_save(T_raw, log_dir)
- T_ret = T_analyze(T_raw, convert_to_K)
- T_plot(T_raw, plotVerbose)
- return T_ret
-
-def T_load_analyze_tweaked(tweak_file=None, convert_to_K=C_to_K,
- plotVerbose=False) :
- "Read the T files listed in tweak_file, return an array of Ts in Kelvin"
- if tweak_file != None :
- Tlist=[]
- for filepath in file(datafile, 'r') :
- cleanpath = filepath.split('\n')[0]
- Tlist.extend(T_load(cleanpath))
- Tlist = numpy.array(Tlist, dtype=numpy.float)
- T_raw = array(Tlist, dtype=numpy.float)
- for Ti in T_raw :
- T_plot(T, plotVerbose)
- else :
- T_raw = T_load(None)
- T = T_analyze(T_raw, convert_to_K)
- return T
-
-# vib family
-
-def vib_aquire(zpiezo, time, freq) :
- """
- Record data for TIME seconds at FREQ Hz from ZPIEZO at it's current position.
- """
- # round up to the nearest power of two, for efficient FFT-ing
- nsamps = ceil_pow_of_two(time*freq)
- time = nsamps / freq
- # take some data, keeping the position voltage at it's current value
- out = ones((nsamps,), dtype=uint16) * zpiezo.curPos()
- if TEXT_VERBOSE :
- print "get %g seconds of data" % time
- data = zpiezo.ramp(out, freq)
- if LOG_DATA :
- log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
- log_name="vib_%gHz" % freq)
- log.write_dict_of_arrays(data)
- return data
-
-def vib_save(data, log_dir=None) :
- "Save the dictionary data, using data_logger.data_log()"
- if log_dir :
- log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
- log_name="vib_%gHz" % freq)
- log.write_dict_of_arrays(data)
-
-def vib_load(datafile=None) :
- "Load the dictionary data, using data_logger.date_load()"
- dl = data_logger.data_load()
- data = dl.read_dict_of_arrays(path)
- return data
-
-def vib_plot(data, freq,
- chunk_size=2048, overlap=True, window=FFT_tools.window_hann,
- plotVerbose=True) :
- """
- If plotVerbose or PYLAB_VERBOSE == True, plot 3 subfigures:
- Time series (Vphoto vs sample index) (show any major drift events),
- A histogram of Vphoto, with a gaussian fit (demonstrate randomness), and
- FFTed Vphoto data (Vphoto vs Frequency) (show major noise components).
- """
- if plotVerbose or PYLAB_VERBOSE :
- _import_pylab()
- _pylab.hold(False)
- _pylab.figure(BASE_FIGNUM+2)
-
- # plot time series
- _pylab.subplot(311)
- _pylab.plot(data["Deflection input"], 'r.')
- _pylab.title("free oscillation")
-
- # plot histogram distribution and gaussian fit
- _pylab.subplot(312)
- n, bins, patches = \
- _pylab.hist(data["Deflection input"], bins=30,
- normed=1, align='center')
- gaus = numpy.zeros((len(bins),))
- mean = data["Deflection input"].mean()
- std = data["Deflection input"].std()
- pi = numpy.pi
- exp = numpy.exp
- for i in range(len(bins)) :
- gaus[i] = (2*pi)**(-.5)/std * exp(-0.5*((bins[i]-mean)/std)**2)
- _pylab.hold(True)
- _pylab.plot(bins, gaus, 'r-');
- _pylab.hold(False)
-
- # plot FFTed data
- _pylab.subplot(313)
- AC_data = data["Deflection input"] - mean
- (freq_axis, power) = \
- FFT_tools.unitary_avg_power_spectrum(AC_data, freq, chunk_size,
- overlap, window)
- _pylab.semilogy(freq_axis, power, 'r.-')
- _pylab.hold(True)
- filtered_power = FFT_tools.windowed_filter(power,
- FFT_tools.windowed_median_point, s=5)
- _pylab.semilogy(freq_axis, filtered_power, 'b.-')
- _flush_plot()
-
-def vib_analyze_naive(data, Vphoto_in2V=_usual_Vphoto_in2V,
- zeroVphoto_bits=_usual_zeroVphoto_bits,
- plotVerbose=False) :
- """
- Calculate the variance of the raw data, and convert it to Volts**2.
- This method is simple and easy to understand, but it highly succeptible
- to noise, drift, etc.
-
- Vphoto_in2V is a function converting Vphoto values from bits to Volts.
- This function is assumed linear, see bump_analyze().
- zeroVphoto_bits is the bit value corresponding to Vphoto = zero Volts.
- """
- Vphoto_std_bits = data["Deflection input"].std()
- Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zpiezo_zeroV_in)
- return Vphoto_std**2
-
-def vib_analyze(data, freq, minFreq=500, maxFreq=7000,
- chunk_size=4096, overlap=True, window=FFT_tools.window_hann,
- median_filter_width=5,
- Vphoto_in2V=_usual_Vphoto_in2V, plotVerbose=False) :
- """
- Calculate the variance in the raw data due to the cantilever's
- thermal oscillation and convert it to Volts**2.
-
- Improves on vib_analyze_naive() by first converting Vphoto(t) to
- frequency space, and fitting a Lorentzian in the relevant frequency
- range (see cantilever_calib.pdf for derivation).
-
- The conversion to frequency space generates an average power spectrum
- by breaking the data into windowed chunks and averaging the power
- spectrums for the chunks together.
- See FFT_tools.avg_power_spectrum
-
- Vphoto_in2V is a function converting Vphoto values from bits to Volts.
- This function is assumed linear, see bump_analyze().
- zeroVphoto_bits is the bit value corresponding to Vphoto = zero Volts.
- """
- Vphoto_t = numpy.zeros((len(data["Deflection input"]),),
- dtype=numpy.float)
- # convert the data from bits to volts
- if TEXT_VERBOSE :
- print "Converting %d bit values to voltages" % len(Vphoto_t)
- for i in range(len(data["Deflection input"])) :
- Vphoto_t[i] = Vphoto_in2V(data["Deflection input"][i])
-
- # Compute the average power spectral density per unit time (in uV**2/Hz)
- if TEXT_VERBOSE :
- print "Compute the averaged power spectral density in uV**2/Hz"
- (freq_axis, power) = \
- FFT_tools.unitary_avg_power_spectrum((Vphoto_t - Vphoto_t.mean())*1e6,
- freq, chunk_size,overlap, window)
-
- # cut out the relevent frequency range
- if TEXT_VERBOSE :
- print "Cut the frequency range %g to %g Hz" % (minFreq, maxFreq)
- imin = 0
- while freq_axis[imin] < minFreq : imin += 1
- imax = imin
- while freq_axis[imax] < maxFreq : imax += 1
- assert imax >= imin + 10 , \
- "Less than 10 points in freq range (%g,%g)" % (minFreq, maxFreq)
-
- # We don't need to filter, because taking data is so cheap...
- # median filter the relevent range
- #filtered_power = FFT_tools.windowed_filter(power[imin:imax],
- # FFT_tools.windowed_median_point,
- # s=median_filter_width)
- filtered_power = power[imin:imax]
-
- # check
- if (plotVerbose or PYLAB_VERBOSE) and False :
- if TEXT_VERBOSE :
- print "Plot the FFT data for checking"
- vib_plot(data, freq, chunk_size, overlap, window, plotVerbose)
- _import_pylab()
- _pylab.figure(5)
- #_pylab.semilogy(freq_axis, power, 'r.-')
- #print "len ", len(freq_axis), len(power), freq_axis.min(), filtered_power.min()
- _pylab.semilogy(freq_axis, power, 'r.-',
- freq_axis[imin:imax], power[imin:imax], 'b.-',
- freq_axis[imin:imax], filtered_power, 'g.-')
- _flush_plot()
-
- # save about-to-be-fitted data to a temporary file
- if TEXT_VERBOSE :
- print "Save data in temp file for fitting"
- datafilename = "%s_%d" % (GNUFIT_DATA_BASE, time.time())
- fd = file(datafilename, 'w')
- for i in range(imin, imax) :
- fd.write("%g\t%g\n" % (freq_axis[i], filtered_power[i-imin]))
- fd.close()
-
- # generate guesses for Lorentzian parameters A,B,C
- if TEXT_VERBOSE :
- print "Fit lorentzian"
- max_fp_index = numpy.argmin(filtered_power)
- max_fp = filtered_power[max_fp_index]
- half_max_index = 0
- while filtered_power[half_max_index] < max_fp/2 :
- half_max_index += 1
- # Lorentzian L(x) = c / ((a**2-x**2)**2 + (b*x)**2)
- # peak centered at a, so
- A_guess = freq_axis[max_fp_index+imin]
- # Half width w on lower side when L(a-w) = L(a)/2
- # (a**2 - (a-w)**2)**2 + (b*(a-w))**2 = 2*(b*a)**2
- # Let W=(a-w)**2, A=a**2, and B=b**2
- # (A - W)**2 + BW = 2*AB
- # W**2 - 2AW + A**2 + BW = 2AB
- # W**2 + (B-2A)W + (A**2-2AB) = 0
- # W = (2A-B)/2 * [1 +/- sqrt(1 - 4(A**2-2AB)/(B-2A)**2]
- # = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
- # (a-w)**2 = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
- # so w is a disaster ;)
- # For some values of A and B (non-underdamped), W is imaginary or negative.
- B_guess = A_guess/2
- # for underdamped, max value at x ~= a, where f(a) = c/(ab)**2, so
- C_guess = max_fp * (A_guess*B_guess)**2
-
- # fit Lorentzian using Gnuplot's 'fit' command
- g = GnuplotBiDir.Gnuplot()
- # The Lorentzian is the predicted one-sided power spectral density
- # For an oscillator whose position x obeys
- # m x'' + gamma x' + kx = F_thermal(t)
- # A is the ?driving noise?
- # B is gamma, the damping term
- # C is omega_0, the resonant frequency
- # Fitting the PSD of V = photoSensitivity*x just rescales the parameters
- g("f(x) = C / ((A**2-x**2)**2 + (x*B)**2)")
- ## The mass m is given by m = k_B T / (pi A)
- ## The spring constant k is given by k = m (omega_0)**2
- ## The quality factor Q is given by Q = omega_0 m / gamma
- g("A=%g; B=%g; C=%g" % (A_guess, B_guess, C_guess))
- g("fit f(x) '%s' via A,B,C" % datafilename)
- A=float(g.getvar('A'))
- B=float(g.getvar('B'))
- C=float(g.getvar('C'))
- if TEXT_VERBOSE :
- print "Fitted f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with"
- print "A = %g, \t B = %g, \t C = %g" % (A, B, C)
- if plotVerbose or GNUPLOT_VERBOSE :
- if TEXT_VERBOSE :
- print "Fitted f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with"
- print "A = %g, \t B = %g, \t C = %g" % (A, B, C)
- # write all the ft data now
- fd = file(datafilename, 'w')
- for i in range(len(freq_axis)) :
- fd.write("%g\t%g\n" % (freq_axis[i], power[i]))
- fd.write("\n") # blank line to separate fit data.
- # write the fit data
- for i in range(imin, imax) :
- x = freq_axis[i]
- fd.write("%g\t%g\n" % (freq_axis[i],
- C / ((A**2-x**2)**2 + (x*B)**2) ) )
- fd.close()
- print Vphoto_t.var(), A, B, C, ((numpy.pi * C)/(2 * B * A**2) * 1e-12)
- g("set terminal epslatex color solid")
- g("set output 'calibration.tex")
- g("set size 2,2") # multipliers on default 5"x3.5"
- #g("set title 'Thermal calibration'")
- g("set logscale y")
- g("set xlabel 'Frequency (Hz)'")
- g("set ylabel 'Power ($\mu$V$^2$/Hz)'")
- g("set xrange [0:20000]")
- g("g(x) = x < %g || x > %g ? 1/0 : f(x)" % (minFreq, maxFreq)) # only plot on fitted region
- g("plot '%s' title 'Vibration data', g(x) title 'Lorentzian fit'" % datafilename)
-
- g("set mouse")
- g("pause mouse 'click with mouse'")
- g.getvar("MOUSE_BUTTON")
- del(g)
-
- # Integrating the the power spectral density per unit time (power) over the
- # frequency axis [0, fN] returns the total signal power per unit time
- # int_0^fN power(f)df = 1/T int_0^T |x(t)**2| dt
- # = <V(t)**2>, the variance for our AC signal.
- # The variance from our fitted Lorentzian is the area under the Lorentzian
- # <V(t)**2> = (pi*C) / (2*B*A**2)
- # Our A is in uV**2, so convert back to Volts
- return (numpy.pi * C) / (2 * B * A**2) * 1e-12
-
-def vib(zpiezo, time=1, freq=50e3, log_dir=None, plotVerbose=False) :
- """
- Wrapper around vib_aquire(), vib_save(), vib_analyze(), vib_plot()
- """
- data = vib_aquire(zpiezo, time, freq)
- vib_save(data, log_dir)
- Vphoto_var = vib_analyze(data, freq)
- vib_plot(data, plotVerbose)
- return Vphoto_var
-
-def vib_load_analyze_tweaked(tweak_file, minFreq=1000, maxFreq=7000,
- chunk_size=2048, overlap=True,
- window=FFT_tools.window_hann,
- median_filter_width=5,
- Vphoto_in2V=_usual_Vphoto_in2V,
- plotVerbose=False) :
- """
- Read the vib files listed in tweak_file.
- Return an array of Vphoto variances (due to the cantilever) in Volts**2
- """
- dl = data_logger.data_load()
- Vphoto_var = []
- for path in file(tweak_file, 'r') :
- if path[-1] == '\n':
- path = path.split('\n')[0]
- # read the data
- data = dl.read_dict_of_arrays(path)
- freq = float(path.split('_')[-1].split('Hz')[0])
- if TEXT_VERBOSE and False :
- print "Analyzing '%s' at %g Hz" % (path, freq)
- # analyze
- Vphoto_var.append(vib_analyze(data, freq, minFreq, maxFreq,
- chunk_size, overlap, window,
- median_filter_width,
- Vphoto_in2V, plotVerbose))
- return numpy.array(Vphoto_var, dtype=numpy.float)
-
-
-# A few positioning functions, so we can run bump_aquire() and vib_aquire()
-# with proper spacing relative to the surface.
-
-def move_just_onto_surface(stepper, zpiezo, Depth_nm=100, setpoint=2) :
- """
- Uses z_piezo_utils.getSurfPos() to pinpoint the position of the surface.
- Adjusts the stepper position as required to get within stepper_tol nm
- of the surface.
- Then set Vzp to place the cantilever Depth_nm onto the surface.
-
- If getSurfPos() fails to find the surface, backs off (for safety)
- and steps in (without moving the zpiezo) until Vphoto > setpoint.
- """
- stepper_tol = 250 # nm, generous estimate of the fullstep stepsize
-
- if TEXT_VERBOSE :
- print "moving just onto surface"
- # Zero the piezo
- if TEXT_VERBOSE :
- print "zero the z piezo output"
- zpiezo.jumpToPos(zpiezo.pos_nm2out(0))
- # See if we're near the surface already
- if TEXT_VERBOSE :
- print "See if we're starting near the surface"
- try :
- dist = zpiezo.pos_out2nm( \
- z_piezo_utils.getSurfPos(zpiezo, zpiezo.def_V2in(setpoint))
- )
- except (z_piezo_utils.tooClose, z_piezo_utils.poorFit), string :
- if TEXT_VERBOSE :
- print "distance failed with: ", string
- print "Back off 200 half steps"
- # Back away 200 steps
- stepper.step_rel(-400)
- stepper.step_rel(200)
- sp = zpiezo.def_V2in(setpoint) # sp = setpoint in bits
- zpiezo.updateInputs()
- cd = zpiezo.curDef() # cd = current deflection in bits
- if TEXT_VERBOSE :
- print "Single stepping approach"
- while cd < sp :
- if TEXT_VERBOSE :
- print "deflection %g < setpoint %g. step closer" % (cd, sp)
- stepper.step_rel(2) # Full step in
- zpiezo.updateInputs()
- cd = zpiezo.curDef()
- # Back off two steps (protecting against backlash)
- if TEXT_VERBOSE :
- print "Step back 4 half steps to get off the setpoint"
- stepper.step_rel(-200)
- stepper.step_rel(196)
- # get the distance to the surface
- zpiezo.updateInputs()
- if TEXT_VERBOSE :
- print "get surf pos, with setpoint %g (%d)" % (setpoint, zpiezo.def_V2in(setpoint))
- for i in range(20) : # HACK, keep stepping back until we get a distance
- try :
- dist = zpiezo.pos_out2nm(getSurfPos(zpiezo, zpiezo.def_V2in(setpoint)))
- except (tooClose, poorFit), string :
- stepper.step_rel(-200)
- stepper.step_rel(198)
- continue
- break
- if i >= 19 :
- print "tried %d times, still too close! bailing" % i
- print "probably an invalid setpoint."
- raise Exception, "weirdness"
- if TEXT_VERBOSE :
- print 'distance to surface ', dist, ' nm'
- # fine tune the stepper position
- while dist < -stepper_tol : # step back if we need to
- stepper.step_rel(-200)
- stepper.step_rel(198)
- dist = zpiezo.pos_out2nm(getSurfPos(zpiezo, zpiezo.def_V2in(setpoint)))
- if TEXT_VERBOSE :
- print 'distance to surface ', dist, ' nm, step back'
- while dist > stepper_tol : # and step forward if we need to
- stepper.step_rel(2)
- dist = zpiezo.pos_out2nm(getSurfPos(zpiezo, zpiezo.def_V2in(setpoint)))
- if TEXT_VERBOSE :
- print 'distance to surface ', dist, ' nm, step closer'
- # now adjust the zpiezo to place us just onto the surface
- target = dist + Depth_nm
- zpiezo.jumpToPos(zpiezo.pos_nm2out(target))
- # and we're there :)
- if TEXT_VERBOSE :
- print "We're %g nm into the surface" % Depth_nm
-
-def move_far_from_surface(stepper, um_back=50) :
- """
- Step back a specified number of microns.
- (uses very rough estimate of step distance at the moment)
- """
- step_nm = 100
- steps = int(um_back*1000/step_nm)
- print "step back %d steps" % steps
- stepper.step_rel(-steps)
-
-
-# and finally, the calib family
-
-def calib_aquire(stepper, zpiezo, get_T=None,
- approach_setpoint=2,
- bump_start_depth=100, bump_push_depth=200,
- bump_npoints=1024, bump_freq=100e3,
- T_convert_to_K=C_to_K,
- vib_time=1, vib_freq=50e3,
- num_bumps=10, num_Ts=10, num_vibs=20,
- log_dir=None, plotVerbose=False) :
- """
- Aquire data for calibrating a cantilever in one function.
- return (bump, T, vib), each of which is an array.
- Inputs :
- stepper a stepper.stepper_obj for coarse Z positioning
- zpiezo a z_piezo.z_piezo for fine positioning and deflection readin
- setpoint maximum allowed deflection (in Volts) during approaches
- num_bumps number of 'a's (see Outputs)
- push_depth_nm depth of each push when generating a
- num_temps number of 'T's (see Outputs)
- num_vibs number of 'vib's (see Outputs)
- log_dir directory to log data to. Default 'None' disables logging.
- Outputs (all are arrays of recorded data) :
- bumps measured (V_photodiode / nm_tip) proportionality constant
- Ts measured temperature (K)
- vibs measured V_photodiode variance in free solution
- """
- # get bumps
- move_just_onto_surface(stepper, zpiezo,
- bump_start_depth, approach_setpoint)
- bumps=zeros((num_bumps,))
- for i in range(num_bumps) :
- bumps[i] = bump(zpiezo, bump_push_depth, bump_npoints, bump_freq,
- log_dir, plotVerbose)
- if TEXT_VERBOSE :
- print bumps
-
- # get Ts
- if get_T == None :
- get_T = lambda : DEFAULT_TEMP
- assert T_convert_to_K == C_to_K, "You didn't define get_T()!"
- Ts=zeros((num_Ts,), dtype=float)
- for i in range(num_Ts) :
- Ts[i] = T(get_T, T_convert_to_K, log_dir, plotVerbose)
- time.sleep(1) # wait a bit to get an independent temperature measure
- print Ts
-
- # get vibs
- move_far_from_surface(stepper)
- vibs=zeros((num_vibs,))
- for i in range(num_vibs) :
- vibs[i] = vib(zpiezo, vib_time, vib_freq, log_dir=log_dir)
- print vibs
-
- if LOG_DATA != None :
- data = {'bump':bumps, 'T':Ts, 'vib':vibs}
- log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
- log_name="calib")
- log.write_dict_of_arrays(data)
-
- return (bumps, Ts, vibs)
-
-def calib_save(bumps, Ts, vibs, log_dir) :
- """
- Save a dictonary with the bump, T, and vib data.
- """
- if log_dir != None :
- data = {'bump':bumps, 'T':Ts, 'vib':vibs}
- log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
- log_name="calib")
- log.write_dict_of_arrays(data)
-
-def calib_load(datafile) :
- "Load the dictionary data, using data_logger.date_load()"
- dl = data_logger.data_load()
- data = dl.read_dict_of_arrays(path)
- return (data['bump'], data['T'], data['vib'])
-
-def calib_save_analysis(k, k_s,
- photoSensitivity2_m, photoSensitivity2_s,
- T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s,
- log_dir=None) :
- if log_dir != None :
- log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
- log_name="calib_analysis_text")
- log.write_binary((
- "k (N/m) : %g +/- %g\n" % (k, k_s)
- + "photoSensitivity**2 (V/nm)**2 : %g +/- %g\n" % \
- (photoSensitivity2_m, photoSensitivity2_s)
- + "T (K) : %g +/- %g\n" % (T_m, T_s)
- + "1/Vp**2 (1/V)**2 : %g +/- %g\n" % \
- (one_o_Vphoto2_m, one_o_Vphoto2_s)
- ))
-
-def calib_plot(bumps, Ts, vibs, plotVerbose) :
- if plotVerbose or PYLAB_VERBOSE :
- _import_pylab()
- _pylab.figure(BASE_FIGNUM+4)
- _pylab.subplot(311)
- _pylab.plot(bumps, 'g.')
- _pylab.title('Photodiode sensitivity (V/nm)')
- _pylab.subplot(312)
- _pylab.plot(Ts, 'r.')
- _pylab.title('Temperature (K)')
- _pylab.subplot(313)
- _pylab.plot(vibs, 'b.')
- _pylab.title('Thermal deflection variance (Volts**2)')
- _flush_plot()
-
-def calib(stepper, zpiezo, tempController=None,
- setpoint=2.0,
- num_bumps=10, push_depth_nm=300,
- num_temps=10,
- num_vibs=10,
- log_dir=None) :
- """
- Calibrate a cantilever in one function.
- The I-don't-care-about-the-details black box version :p.
- return (k, k_s)
- Inputs:
- (see calib_aquire())
- Outputs :
- k cantilever spring constant (in N/m, or equivalently nN/nm)
- k_s standard deviation in our estimate of k
- Notes :
- See get_calibration_data() for the data aquisition code
- See analyze_calibration_data() for the analysis code
- """
- a, T, vib = calib_aquire(stepper, zpiezo, tempController,
- setpoint,
- num_bumps=num_bumps,
- push_depth_nm=push_depth_nm,
- num_temps=num_temps,
- num_vibs=num_vibs,
- log_dir=log_dir)
- calib_save(a, T, vib, log_dir)
- k,k_s,ps2_m, ps2_s,T_m,T_s,one_o_Vp2_m,one_o_Vp2_s = \
- calib_analyze(a, T, vib, log_dir=log_dir)
- calib_save_analysis(k,k_s,ps2_m, ps2_s,T_m,T_s,one_o_Vp2_m,one_o_Vp2_s,
- log_dir)
- calib_plot(a, T, vib, plotVerbose)
- return (k, k_s)
-
-def calib_load_analyze_tweaks(bump_tweaks, vib_tweaks,
- T_tweaks=None,
- log_dir=None,
- plotVerbose=True) :
- a = read_tweaked_bumps(bump_tweaks)
- vib = V_photo_variance_from_file(vib_tweaks)
- return analyze_calibration_data(a, T, vib, log_dir=log_dir)
-
-
PYLAB_INTERACTIVE = True # select between draw() and show() for flushing plots
BASE_FIGNUM = 20 # to avoid writing to already existing figures
+
+# 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 for conversion
+# functions
+
# zpGain zpiezo applied voltage per output Volt
zpGain = z_piezo.DEFAULT_GAIN
# zpSensitivity nm zpiezo response per applied Volt
import config # config module for the calibcant package
import data_logger
import FFT_tools
+from splittable_kwargs import splittableKwargsFunction, \
+ make_splittable_kwargs_function
PLOT_GUESSED_LORENTZIAN=False
+@splittableKwargsFunction()
def vib_analyze_naive(deflection_bits, Vphoto_in2V=config.Vphoto_in2V,
zeroVphoto_bits=config.zeroVphoto_bits) :
"""
Vphoto_std_bits = deflection_bits.std()
Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zeroVphoto_bits)
return Vphoto_std**2
-
-def vib_analyze(deflection_bits, freq, minFreq=500, maxFreq=7000,
+
+#@splittableKwargsFunction((fit_psd, 'freq_axis', 'psd_data'),
+# (vib_plot, 'deflection_bits', 'freq_axis', 'power',
+# 'A', 'B', 'C', 'minFreq', 'maxFreq'))
+# Some of the child functions aren't yet defined, so postpone
+# make-splittable until later in the module.
+def vib_analyze(deflection_bits, freq,
chunk_size=4096, overlap=True, window=FFT_tools.window_hann,
- Vphoto_in2V=config.Vphoto_in2V, plotVerbose=False) :
+ Vphoto_in2V=config.Vphoto_in2V, **kwargs) :
"""
Calculate the variance in the raw data due to the cantilever's
thermal oscillation and convert it to Volts**2.
Vphoto_in2V is a function converting Vphoto values from bits to Volts.
This function is assumed linear, see bump_analyze().
"""
+ fit_psd_kwargs,vib_plot_kwargs = vib_analyze._splitargs(vib_analyze,kwargs)
+ if 'minFreq' in fit_psd_kwargs:
+ vib_plot_kwargs['minFreq'] = fit_psd_kwargs['minFreq']
+ if 'maxFreq' in fit_psd_kwargs:
+ vib_plot_kwargs['maxFreq'] = fit_psd_kwargs['maxFreq']
+
Vphoto_t = numpy.zeros((len(deflection_bits),),
dtype=numpy.float)
# convert the data from bits to volts
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)
+ A,B,C = fit_psd(freq_axis, power, **kwargs)
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,
- minFreq, maxFreq, plotVerbose=plotVerbose)
+ vib_plot(deflection_bits, freq_axis, power, A, B, C, **kwargs)
# Our A is in uV**2, so convert back to Volts**2
fitted_variance = lorentzian_area(A,B,C) * 1e-12
return min(fitted_variance, naive_variance)
+@splittableKwargsFunction()
def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) :
"""
freq_axis : array of frequencies in Hz
string += "plot '%s', f(x)\n" % datafilename
return string
+@splittableKwargsFunction()
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
'Deflection input'
"""
if log_dir :
+ freq = data['sample frequency Hz']
log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
log_name="vib_%gHz" % freq)
log.write_dict_of_arrays(data)
data = dl.read_dict_of_arrays(datafile)
return data
+@splittableKwargsFunction()
def vib_plot(deflection_bits, freq_axis, power, A, B, C,
- minFreq=None, maxFreq=None, plotVerbose=True) :
+ minFreq=None, maxFreq=None, plotVerbose=False) :
"""
If plotVerbose or config.PYLAB_VERBOSE == True, plot 3 subfigures:
Time series (Vphoto vs sample index) (show any major drift events),
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) :
+# Make vib_analyze splittable (was postponed until the child functions
+# were defined).
+make_splittable_kwargs_function(vib_analyze,
+ (fit_psd, 'freq_axis', 'psd_data'),
+ (vib_plot, 'deflection_bits', 'freq_axis', 'power', 'A', 'B', 'C',
+ 'minFreq', 'maxFreq'))
+
+@splittableKwargsFunction((vib_analyze, 'deflection_bits', 'freq'))
+def vib_load_analyze_tweaked(tweak_file, **kwargs) :
"""
Read the vib files listed in tweak_file.
Return an array of Vphoto variances (due to the cantilever) in Volts**2
"""
+ vib_analyze_kwargs, = \
+ vib_load_analyze_tweaked._splitargs(vib_load_analyze_tweaked, kwargs)
dl = data_logger.data_load()
Vphoto_var = []
for path in file(tweak_file, 'r') :
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))
+ Vphoto_var.append(vib_analyze(deflection_bits, freq,
+ **vib_analyze_kwargs))
return numpy.array(Vphoto_var, dtype=numpy.float)
# commandline interface functions
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_INTERACTIVE = False
config.PYLAB_VERBOSE = options.pylab
projects / calibcant.git / commitdiff