4 Aquire and analyze cantilever calibration data.
6 W. Trevor King Dec. 2007-Jan. 2008
11 The relevent physical quantities are :
12 Vzp_out Output z-piezo voltage (what we generate)
13 Vzp Applied z-piezo voltage (after external ZPGAIN)
14 Zp The z-piezo position
15 Zcant The cantilever vertical deflection
16 Vphoto The photodiode vertical deflection voltage (what we measure)
17 Fcant The force on the cantilever
18 T The temperature of the cantilever and surrounding solution
19 (another thing we measure)
20 k_b Boltzmann's constant
22 Which are related by the parameters :
24 zpSensitivity Zp / Vzp
25 photoSensitivity Vphoto / Zcant
28 Cantilever calibration assumes a pre-calibrated z-piezo (zpSensitivity) and
30 In our lab, the z-piezo is calibrated by imaging a calibration sample,
31 which has features with well defined sizes, and the gain is set with a knob
34 photoSensitivity is measured by bumping the cantilever against the surface,
35 where Zp = Zcant (see the bump_*() family of functions)
36 The measured slope Vphoto/Vout is converted to photoSensitivity via
37 Vphoto/Vzp_out * Vzp_out/Vzp * Vzp/Zp * Zp/Zcant = Vphoto/Zcant
38 (measured) (1/zpGain) (1/zpSensitivity) (1) (photoSensitivity)
40 k_cant is measured by watching the cantilever vibrate in free solution
41 (see the vib_*() family of functions)
42 The average energy of the cantilever in the vertical direction is given
43 by the equipartition theorem.
44 1/2 k_b T = 1/2 k_cant Zcant**2
45 so k_cant = k_b T / Zcant**2
46 but Zcant = Vphoto / photoSensitivity
47 so k_cant = k_b T * photoSensitivty**2 / Vphoto**2
48 We measured photoSensitivity above.
49 We can either measure T using an external function (see temperature.py),
50 or just estimate it (see the T_*() family of functions).
51 (room temp ~22 deg C, accurate to +/- 5 deg, so 5/273.15 ~= 2%)
52 Finally, a time series of Vphoto while we're far from the surface and not
53 changing Vzp_out will give us Vphoto.
55 We do all these measurements a few times to estimate statistical errors.
57 The functions are layed out in the families:
58 bump_*(), vib_*(), T_*(), and calib_*()
59 where calib_{save|load|analyze}() deal with derived data, not real-world data.
60 For each family, * can be any of :
61 aquire get real-world data
62 save store real-world data to disk
63 load get real-world data from disk
64 analyze interperate the real-world data.
65 plot show a nice graphic to convince people we're working :p
67 read a file with a list of paths to previously saved real world data
68 load each file using *_load(), analyze using *_analyze(), and
69 optionally plot using *_plot().
70 Intended for re-processing old data.
71 A family name without any _* extension (e.g. bump()),
72 runs *_aquire(), *_save(), *_analyze(), *_plot().
74 Also defines the two positioning functions:
75 move_just_onto_surface() and move_far_from_surface()
84 import GnuplotBiDir # used for fitting lorentzian
86 kb = 1.3806504e-23 # Boltzmann's constant in J/K
87 DEFAULT_TEMP = 22 # assume 22 deg C
89 LOG_DATA = True # quietly grab all real-world data and log to LOG_DIR
90 LOG_DIR = '/home/wking/rsrch/data/calibrate_cantilever'
92 GNUFIT_DATA_BASE='./calibrate_cantilever_fitdata'
94 TEXT_VERBOSE = True # for debugging
95 GNUPLOT_VERBOSE = True # turn on fit check plotting
96 PYLAB_VERBOSE = True # turn on plotting
97 PYLAB_INTERACTIVE = True # select between draw() and show() for flushing plots
98 BASE_FIGNUM = 20 # to avoid writing to already existing figures
100 # handle extra verbose input modules, only imported if we need them
102 _final_flush_plot = None
104 def _import_pylab() :
105 "Import pylab plotting functions for when we need to plot"
108 global _final_flush_plot
110 import pylab as _pylab
111 if _flush_plot == None :
112 if PYLAB_INTERACTIVE :
113 _flush_plot = _pylab.draw
115 def _flush_plot () : pass
116 if _final_flush_plot == None :
117 if PYLAB_INTERACTIVE :
118 _final_flush_plot = _pylab.draw
120 _final_flush_plot = _pylab.show
124 # make sure you make a system note (add_system_note) if you change these
125 # in case you don't have access to a z_piezo.z_piezo for conversion functions
127 _usual_zpSensitivity=7.41
128 _usual_zeroVphoto_bits = 2**15
129 _usual_zeroVzp_bits = 2**15
130 def _usual_Vphoto_in2V(inp) :
131 return (float(inp)-float(_usual_zeroVphoto_bits))*(10.0/2**15)
132 def _usual_Vzp_out2V(inp) :
133 return (float(inp)-float(_usual_zeroVzp_bits))*(10.0/2**15)
135 def C_to_K(celsius) :
136 "Convert Celsius -> Kelvin."
137 return celsius + 273.15
141 def bump_aquire(zpiezo, push_depth, npoints, freq) :
143 Ramps closer push_depth and returns to the original position.
145 zpiezo an opened zpiezo.zpiezo instance
146 push_depth distance to approach, in nm
147 npoints number points during the approach and during the retreat
148 freq rate at which data is aquired
149 log_dir directory to log data to (see data_logger.py).
150 None to turn off logging (see also the global LOG_DATA).
151 Returns the aquired ramp data dictionary, with data in DAC/ADC bits.
153 # generate the bump output
154 start_pos = zpiezo.curPos()
155 pos_dist = zpiezo.pos_nm2out(push_depth) - zpiezo.pos_nm2out(0)
156 close_pos = start_pos + pos_dist
157 appr = linspace(start_pos, close_pos, npoints)
158 retr = linspace(close_pos, start_pos, npoints)
159 out = concatenate((appr, retr))
160 # run the bump, and measure deflection
162 print "Bump %g nm" % push_depth
163 data = zpiezo.ramp(out, freq)
164 # default saving, so we have a log in-case the operator is lazy ;)
165 if LOG_DATA == True :
166 log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
167 log_name="bump_surface")
168 log.write_dict_of_arrays(data)
171 def bump_save(data, log_dir) :
172 "Save the dictionary data, using data_logger.data_log()"
174 log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
176 log.write_dict_of_arrays(data)
178 def bump_load(datafile) :
179 "Load the dictionary data, using data_logger.date_load()"
180 dl = data_logger.data_load()
181 data = dl.read_dict_of_arrays(path)
184 def bump_analyze(data, zpGain=_usual_zpGain,
185 zpSensitivity=_usual_zpSensitivity,
186 Vzp_out2V=_usual_Vzp_out2V,
187 Vphoto_in2V=_usual_Vphoto_in2V,
190 Return the slope of the bump ;).
192 data dictinary of data in DAC/ADC bits
193 pos_out2V function that converts output DAC bits to Volts
194 def_in2V funtion that converts input ADC bits to Volts
195 zpGain zpiezo applied voltage per output Volt
196 zpSensitivity nm zpiezo response per applied Volt
198 photoSensitivity (Vphoto/Zcant) in Volts/nm
199 Checks for strong correlation (r-value) and low randomness chance (p-value)
201 With the current implementation, the data is regressed in DAC/ADC bits
202 and THEN converted, so we're assuming that both conversions are LINEAR.
203 if they aren't, rewrite to convert before the regression.
205 scale_Vzp_bits2V = Vzp_out2V(1) - Vzp_out2V(0)
206 scale_Vphoto_bits2V = Vzp_in2V(1) - Vphoto_in2V(0)
207 Vphoto2Vzp_out_bit, intercept = \
208 linfit.linregress(x=data["Z piezo output"],
209 y=data["Deflection input"],
210 plotVerbose=plotVerbose)
211 Vphoto2Vzp_out = Vphoto2Vzp_out * scale_Vphoto_bits2V / scale_Vzp_bits2V
213 # 1 / (Vzp/Vzp_out * Zp/Vzp * Zcant/Zp )
214 Vzp_out2Zcant = 1.0/ (zpiezo_gain * zpiezo_sensitivity) # * 1
215 return Vphoto2Vzp_out * Vzp_out2Zcant
217 def bump_plot(data, plotVerbose) :
218 "Plot the bump (Vphoto vs Vzp) if plotVerbose or PYLAB_VERBOSE == True"
219 if plotVerbose or PYLAB_VERBOSE :
221 _pylab.figure(BASE_FIGNUM)
222 _pylab.plot(data["Z piezo output"], data["Deflection input"],
224 _pylab.title("bump surface")
225 _pylab.legend(loc='upper left')
228 def bump(zpiezo, push_depth, npoints=1024, freq=100e3,
232 Wrapper around bump_aquire(), bump_save(), bump_analyze(), bump_plot()
234 data = bump_aquire(zpiezo, push_depth, npoints, freq)
235 bump_save(data, log_dir)
236 photoSensitivity = bump_analyze(data, zpiezo.gain, zpiezo.sensitivity,
237 zpiezo.pos_out2V, zpiezo.def_in2V)
238 bump_plot(data, plotVerbose)
239 return photoSensitivity
241 def bump_load_analyze_tweaked(tweak_file, zpGain=_usual_zpGain,
242 zpSensitivity=_usual_zpSensitivity,
243 Vzp_out2V=_usual_Vzp_out2V,
244 Vphoto_in2V=_usual_Vphoto_in2V,
246 "Load the output file of tweak_calib_bump.sh, return an array of slopes"
247 photoSensitivity = []
248 for line in file(tweak_file, 'r') :
249 parsed = line.split()
250 path = parsed[0].split('\n')[0]
252 full_data = bump_load(path)
253 if len(parsed) == 1 :
254 data = full_data # use whole bump
256 # use the listed sections
259 for rng in parsed[1:] :
263 zp.extend(full_data['Z piezo output'][starti:stopi])
264 df.extend(full_data['Deflection input'][starti:stopi])
265 data = {'Z piezo output': array(zp),
266 'Deflection input':array(df)}
267 pSi = bump_analyze(data, zpGain, zpSensitivity,
268 Vzp_out2V, Vphoto_in2V, plotVerbose)
269 photoSensitivity.append(pSi)
270 bump_plot(data, plotVervose)
271 return array(photoSensitivity, dtype=numpy.float)
274 # Fairly stubby, since a one shot Temp measurement is a common thing.
275 # We just wrap that to provide a consistent interface.
277 def T_aquire(get_T=None) :
279 Measure the current temperature of the sample,
280 or, if get_T == None, fake it by returning DEFAULT_TEMP
284 print "Fake temperature %g" % DEFAULT_TEMP
288 print "Measure temperature"
291 def T_save(T, log_dir=None) :
293 Save either a single T (if you are crazy :p),
294 or an array of them (more reasonable)
296 T = numpy.array(T, dtype=numpy.float)
298 log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
300 log.write_binary(T.tostring())
301 if LOG_DATA != None :
302 log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
304 log.write_binary(T.tostring())
306 def T_load(datafile=None) :
308 Load the saved T array (possibly of length 1), and return it.
309 If datafile == None, return an array of one DEFAULT_TEMP instead.
311 if datafile == None :
312 return numpy.array([DEFAULT_TEMP], dtype=numpy.float)
314 return fromfile(file=cleanpath, dtype=numpy.float)
316 def T_analyze(T, convert_to_K=C_to_K) :
318 Not much to do here, just convert to Kelvin.
319 Uses C_to_K (defined above) by default.
321 try : # if T is an array, convert element by element
322 for i in range(len(T)) :
323 T[i] = convert_to_K(T[i])
328 def T_plot(T, plotVerbose) :
330 Umm, just print the temperature?
332 if plotVerbose or PYLAB_VERBOSE or TEXT_VERBOSE :
333 print "Temperature ", T
335 def T(get_T=None, convert_to_K=C_to_K, log_dir=None, plotVerbose=None) :
337 Wrapper around T_aquire(), T_save(), T_analyze(), T_plot()
339 T_raw = T_aquire(get_T)
340 T_save(T_raw, log_dir)
341 T_ret = T_analyze(T_raw, convert_to_K)
342 T_plot(T_raw, plotVerbose)
345 def T_load_analyze_tweaked(tweak_file=None, convert_to_K=C_to_K,
347 "Read the T files listed in tweak_file, return an array of Ts in Kelvin"
348 if tweak_file != None :
350 for filepath in file(datafile, 'r') :
351 cleanpath = filepath.split('\n')[0]
352 Tlist.extend(T_load(cleanpath))
353 Tlist = numpy.array(Tlist, dtype=numpy.float)
354 T_raw = array(Tlist, dtype=numpy.float)
356 T_plot(T, plotVerbose)
359 T = T_analyze(T_raw, convert_to_K)
364 def vib_aquire(zpiezo, time, freq) :
366 Record data for TIME seconds at FREQ Hz from ZPIEZO at it's current position.
368 # round up to the nearest power of two, for efficient FFT-ing
369 nsamps = ceil_pow_of_two(time*freq)
371 # take some data, keeping the position voltage at it's current value
372 out = ones((nsamps,), dtype=uint16) * zpiezo.curPos()
374 print "get %g seconds of data" % time
375 data = zpiezo.ramp(out, freq)
377 log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
378 log_name="vib_%gHz" % freq)
379 log.write_dict_of_arrays(data)
382 def vib_save(data, log_dir=None) :
383 "Save the dictionary data, using data_logger.data_log()"
385 log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
386 log_name="vib_%gHz" % freq)
387 log.write_dict_of_arrays(data)
389 def vib_load(datafile=None) :
390 "Load the dictionary data, using data_logger.date_load()"
391 dl = data_logger.data_load()
392 data = dl.read_dict_of_arrays(path)
395 def vib_plot(data, freq,
396 chunk_size=2048, overlap=True, window=FFT_tools.window_hann,
399 If plotVerbose or PYLAB_VERBOSE == True, plot 3 subfigures:
400 Time series (Vphoto vs sample index) (show any major drift events),
401 A histogram of Vphoto, with a gaussian fit (demonstrate randomness), and
402 FFTed Vphoto data (Vphoto vs Frequency) (show major noise components).
404 if plotVerbose or PYLAB_VERBOSE :
407 _pylab.figure(BASE_FIGNUM+2)
411 _pylab.plot(data["Deflection input"], 'r.')
412 _pylab.title("free oscillation")
414 # plot histogram distribution and gaussian fit
417 _pylab.hist(data["Deflection input"], bins=30,
418 normed=1, align='center')
419 gaus = numpy.zeros((len(bins),))
420 mean = data["Deflection input"].mean()
421 std = data["Deflection input"].std()
424 for i in range(len(bins)) :
425 gaus[i] = (2*pi)**(-.5)/std * exp(-0.5*((bins[i]-mean)/std)**2)
427 _pylab.plot(bins, gaus, 'r-');
432 AC_data = data["Deflection input"] - mean
433 (freq_axis, power) = \
434 FFT_tools.unitary_avg_power_spectrum(AC_data, freq, chunk_size,
436 _pylab.semilogy(freq_axis, power, 'r.-')
438 filtered_power = FFT_tools.windowed_filter(power,
439 FFT_tools.windowed_median_point, s=5)
440 _pylab.semilogy(freq_axis, filtered_power, 'b.-')
443 def vib_analyze_naive(data, Vphoto_in2V=_usual_Vphoto_in2V,
444 zeroVphoto_bits=_usual_zeroVphoto_bits,
447 Calculate the variance of the raw data, and convert it to Volts**2.
448 This method is simple and easy to understand, but it highly succeptible
449 to noise, drift, etc.
451 Vphoto_in2V is a function converting Vphoto values from bits to Volts.
452 This function is assumed linear, see bump_analyze().
453 zeroVphoto_bits is the bit value corresponding to Vphoto = zero Volts.
455 Vphoto_std_bits = data["Deflection input"].std()
456 Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zpiezo_zeroV_in)
459 def vib_analyze(data, freq, minFreq=500, maxFreq=7000,
460 chunk_size=4096, overlap=True, window=FFT_tools.window_hann,
461 median_filter_width=5,
462 Vphoto_in2V=_usual_Vphoto_in2V, plotVerbose=False) :
464 Calculate the variance in the raw data due to the cantilever's
465 thermal oscillation and convert it to Volts**2.
467 Improves on vib_analyze_naive() by first converting Vphoto(t) to
468 frequency space, and fitting a Lorentzian in the relevant frequency
469 range (see cantilever_calib.pdf for derivation).
471 The conversion to frequency space generates an average power spectrum
472 by breaking the data into windowed chunks and averaging the power
473 spectrums for the chunks together.
474 See FFT_tools.avg_power_spectrum
476 Vphoto_in2V is a function converting Vphoto values from bits to Volts.
477 This function is assumed linear, see bump_analyze().
478 zeroVphoto_bits is the bit value corresponding to Vphoto = zero Volts.
480 Vphoto_t = numpy.zeros((len(data["Deflection input"]),),
482 # convert the data from bits to volts
484 print "Converting %d bit values to voltages" % len(Vphoto_t)
485 for i in range(len(data["Deflection input"])) :
486 Vphoto_t[i] = Vphoto_in2V(data["Deflection input"][i])
488 # Compute the average power spectral density per unit time (in uV**2/Hz)
490 print "Compute the averaged power spectral density in uV**2/Hz"
491 (freq_axis, power) = \
492 FFT_tools.unitary_avg_power_spectrum((Vphoto_t - Vphoto_t.mean())*1e6,
493 freq, chunk_size,overlap, window)
495 # cut out the relevent frequency range
497 print "Cut the frequency range %g to %g Hz" % (minFreq, maxFreq)
499 while freq_axis[imin] < minFreq : imin += 1
501 while freq_axis[imax] < maxFreq : imax += 1
502 assert imax >= imin + 10 , \
503 "Less than 10 points in freq range (%g,%g)" % (minFreq, maxFreq)
505 # We don't need to filter, because taking data is so cheap...
506 # median filter the relevent range
507 #filtered_power = FFT_tools.windowed_filter(power[imin:imax],
508 # FFT_tools.windowed_median_point,
509 # s=median_filter_width)
510 filtered_power = power[imin:imax]
513 if (plotVerbose or PYLAB_VERBOSE) and False :
515 print "Plot the FFT data for checking"
516 vib_plot(data, freq, chunk_size, overlap, window, plotVerbose)
519 #_pylab.semilogy(freq_axis, power, 'r.-')
520 #print "len ", len(freq_axis), len(power), freq_axis.min(), filtered_power.min()
521 _pylab.semilogy(freq_axis, power, 'r.-',
522 freq_axis[imin:imax], power[imin:imax], 'b.-',
523 freq_axis[imin:imax], filtered_power, 'g.-')
526 # save about-to-be-fitted data to a temporary file
528 print "Save data in temp file for fitting"
529 datafilename = "%s_%d" % (GNUFIT_DATA_BASE, time.time())
530 fd = file(datafilename, 'w')
531 for i in range(imin, imax) :
532 fd.write("%g\t%g\n" % (freq_axis[i], filtered_power[i-imin]))
535 # generate guesses for Lorentzian parameters A,B,C
537 print "Fit lorentzian"
538 max_fp_index = numpy.argmin(filtered_power)
539 max_fp = filtered_power[max_fp_index]
541 while filtered_power[half_max_index] < max_fp/2 :
543 # Lorentzian L(x) = c / ((a**2-x**2)**2 + (b*x)**2)
544 # peak centered at a, so
545 A_guess = freq_axis[max_fp_index+imin]
546 # Half width w on lower side when L(a-w) = L(a)/2
547 # (a**2 - (a-w)**2)**2 + (b*(a-w))**2 = 2*(b*a)**2
548 # Let W=(a-w)**2, A=a**2, and B=b**2
549 # (A - W)**2 + BW = 2*AB
550 # W**2 - 2AW + A**2 + BW = 2AB
551 # W**2 + (B-2A)W + (A**2-2AB) = 0
552 # W = (2A-B)/2 * [1 +/- sqrt(1 - 4(A**2-2AB)/(B-2A)**2]
553 # = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
554 # (a-w)**2 = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
555 # so w is a disaster ;)
556 # For some values of A and B (non-underdamped), W is imaginary or negative.
558 # for underdamped, max value at x ~= a, where f(a) = c/(ab)**2, so
559 C_guess = max_fp * (A_guess*B_guess)**2
561 # fit Lorentzian using Gnuplot's 'fit' command
562 g = GnuplotBiDir.Gnuplot()
563 # The Lorentzian is the predicted one-sided power spectral density
564 # For an oscillator whose position x obeys
565 # m x'' + gamma x' + kx = F_thermal(t)
566 # A is the ?driving noise?
567 # B is gamma, the damping term
568 # C is omega_0, the resonant frequency
569 # Fitting the PSD of V = photoSensitivity*x just rescales the parameters
570 g("f(x) = C / ((A**2-x**2)**2 + (x*B)**2)")
571 ## The mass m is given by m = k_B T / (pi A)
572 ## The spring constant k is given by k = m (omega_0)**2
573 ## The quality factor Q is given by Q = omega_0 m / gamma
574 g("A=%g; B=%g; C=%g" % (A_guess, B_guess, C_guess))
575 g("fit f(x) '%s' via A,B,C" % datafilename)
576 A=float(g.getvar('A'))
577 B=float(g.getvar('B'))
578 C=float(g.getvar('C'))
580 print "Fitted f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with"
581 print "A = %g, \t B = %g, \t C = %g" % (A, B, C)
582 if plotVerbose or GNUPLOT_VERBOSE :
584 print "Fitted f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with"
585 print "A = %g, \t B = %g, \t C = %g" % (A, B, C)
586 # write all the ft data now
587 fd = file(datafilename, 'w')
588 for i in range(len(freq_axis)) :
589 fd.write("%g\t%g\n" % (freq_axis[i], power[i]))
590 fd.write("\n") # blank line to separate fit data.
592 for i in range(imin, imax) :
594 fd.write("%g\t%g\n" % (freq_axis[i],
595 C / ((A**2-x**2)**2 + (x*B)**2) ) )
597 print Vphoto_t.var(), A, B, C, ((numpy.pi * C)/(2 * B * A**2) * 1e-12)
598 g("set terminal epslatex color solid")
599 g("set output 'calibration.tex")
600 g("set size 2,2") # multipliers on default 5"x3.5"
601 #g("set title 'Thermal calibration'")
603 g("set xlabel 'Frequency (Hz)'")
604 g("set ylabel 'Power ($\mu$V$^2$/Hz)'")
605 g("set xrange [0:20000]")
606 g("g(x) = x < %g || x > %g ? 1/0 : f(x)" % (minFreq, maxFreq)) # only plot on fitted region
607 g("plot '%s' title 'Vibration data', g(x) title 'Lorentzian fit'" % datafilename)
610 g("pause mouse 'click with mouse'")
611 g.getvar("MOUSE_BUTTON")
614 # Integrating the the power spectral density per unit time (power) over the
615 # frequency axis [0, fN] returns the total signal power per unit time
616 # int_0^fN power(f)df = 1/T int_0^T |x(t)**2| dt
617 # = <V(t)**2>, the variance for our AC signal.
618 # The variance from our fitted Lorentzian is the area under the Lorentzian
619 # <V(t)**2> = (pi*C) / (2*B*A**2)
620 # Our A is in uV**2, so convert back to Volts
621 return (numpy.pi * C) / (2 * B * A**2) * 1e-12
623 def vib(zpiezo, time=1, freq=50e3, log_dir=None, plotVerbose=False) :
625 Wrapper around vib_aquire(), vib_save(), vib_analyze(), vib_plot()
627 data = vib_aquire(zpiezo, time, freq)
628 vib_save(data, log_dir)
629 Vphoto_var = vib_analyze(data, freq)
630 vib_plot(data, plotVerbose)
633 def vib_load_analyze_tweaked(tweak_file, minFreq=1000, maxFreq=7000,
634 chunk_size=2048, overlap=True,
635 window=FFT_tools.window_hann,
636 median_filter_width=5,
637 Vphoto_in2V=_usual_Vphoto_in2V,
640 Read the vib files listed in tweak_file.
641 Return an array of Vphoto variances (due to the cantilever) in Volts**2
643 dl = data_logger.data_load()
645 for path in file(tweak_file, 'r') :
647 path = path.split('\n')[0]
649 data = dl.read_dict_of_arrays(path)
650 freq = float(path.split('_')[-1].split('Hz')[0])
651 if TEXT_VERBOSE and False :
652 print "Analyzing '%s' at %g Hz" % (path, freq)
654 Vphoto_var.append(vib_analyze(data, freq, minFreq, maxFreq,
655 chunk_size, overlap, window,
657 Vphoto_in2V, plotVerbose))
658 return numpy.array(Vphoto_var, dtype=numpy.float)
661 # A few positioning functions, so we can run bump_aquire() and vib_aquire()
662 # with proper spacing relative to the surface.
664 def move_just_onto_surface(stepper, zpiezo, Depth_nm=100, setpoint=2) :
666 Uses z_piezo_utils.getSurfPos() to pinpoint the position of the surface.
667 Adjusts the stepper position as required to get within stepper_tol nm
669 Then set Vzp to place the cantilever Depth_nm onto the surface.
671 If getSurfPos() fails to find the surface, backs off (for safety)
672 and steps in (without moving the zpiezo) until Vphoto > setpoint.
674 stepper_tol = 250 # nm, generous estimate of the fullstep stepsize
677 print "moving just onto surface"
680 print "zero the z piezo output"
681 zpiezo.jumpToPos(zpiezo.pos_nm2out(0))
682 # See if we're near the surface already
684 print "See if we're starting near the surface"
686 dist = zpiezo.pos_out2nm( \
687 z_piezo_utils.getSurfPos(zpiezo, zpiezo.def_V2in(setpoint))
689 except (z_piezo_utils.tooClose, z_piezo_utils.poorFit), string :
691 print "distance failed with: ", string
692 print "Back off 200 half steps"
693 # Back away 200 steps
694 stepper.step_rel(-400)
695 stepper.step_rel(200)
696 sp = zpiezo.def_V2in(setpoint) # sp = setpoint in bits
697 zpiezo.updateInputs()
698 cd = zpiezo.curDef() # cd = current deflection in bits
700 print "Single stepping approach"
703 print "deflection %g < setpoint %g. step closer" % (cd, sp)
704 stepper.step_rel(2) # Full step in
705 zpiezo.updateInputs()
707 # Back off two steps (protecting against backlash)
709 print "Step back 4 half steps to get off the setpoint"
710 stepper.step_rel(-200)
711 stepper.step_rel(196)
712 # get the distance to the surface
713 zpiezo.updateInputs()
715 print "get surf pos, with setpoint %g (%d)" % (setpoint, zpiezo.def_V2in(setpoint))
716 for i in range(20) : # HACK, keep stepping back until we get a distance
718 dist = zpiezo.pos_out2nm(getSurfPos(zpiezo, zpiezo.def_V2in(setpoint)))
719 except (tooClose, poorFit), string :
720 stepper.step_rel(-200)
721 stepper.step_rel(198)
725 print "tried %d times, still too close! bailing" % i
726 print "probably an invalid setpoint."
727 raise Exception, "weirdness"
729 print 'distance to surface ', dist, ' nm'
730 # fine tune the stepper position
731 while dist < -stepper_tol : # step back if we need to
732 stepper.step_rel(-200)
733 stepper.step_rel(198)
734 dist = zpiezo.pos_out2nm(getSurfPos(zpiezo, zpiezo.def_V2in(setpoint)))
736 print 'distance to surface ', dist, ' nm, step back'
737 while dist > stepper_tol : # and step forward if we need to
739 dist = zpiezo.pos_out2nm(getSurfPos(zpiezo, zpiezo.def_V2in(setpoint)))
741 print 'distance to surface ', dist, ' nm, step closer'
742 # now adjust the zpiezo to place us just onto the surface
743 target = dist + Depth_nm
744 zpiezo.jumpToPos(zpiezo.pos_nm2out(target))
747 print "We're %g nm into the surface" % Depth_nm
749 def move_far_from_surface(stepper, um_back=50) :
751 Step back a specified number of microns.
752 (uses very rough estimate of step distance at the moment)
755 steps = int(um_back*1000/step_nm)
756 print "step back %d steps" % steps
757 stepper.step_rel(-steps)
760 # and finally, the calib family
762 def calib_aquire(stepper, zpiezo, get_T=None,
764 bump_start_depth=100, bump_push_depth=200,
765 bump_npoints=1024, bump_freq=100e3,
766 T_convert_to_K=C_to_K,
767 vib_time=1, vib_freq=50e3,
768 num_bumps=10, num_Ts=10, num_vibs=20,
769 log_dir=None, plotVerbose=False) :
771 Aquire data for calibrating a cantilever in one function.
772 return (bump, T, vib), each of which is an array.
774 stepper a stepper.stepper_obj for coarse Z positioning
775 zpiezo a z_piezo.z_piezo for fine positioning and deflection readin
776 setpoint maximum allowed deflection (in Volts) during approaches
777 num_bumps number of 'a's (see Outputs)
778 push_depth_nm depth of each push when generating a
779 num_temps number of 'T's (see Outputs)
780 num_vibs number of 'vib's (see Outputs)
781 log_dir directory to log data to. Default 'None' disables logging.
782 Outputs (all are arrays of recorded data) :
783 bumps measured (V_photodiode / nm_tip) proportionality constant
784 Ts measured temperature (K)
785 vibs measured V_photodiode variance in free solution
788 move_just_onto_surface(stepper, zpiezo,
789 bump_start_depth, approach_setpoint)
790 bumps=zeros((num_bumps,))
791 for i in range(num_bumps) :
792 bumps[i] = bump(zpiezo, bump_push_depth, bump_npoints, bump_freq,
793 log_dir, plotVerbose)
799 get_T = lambda : DEFAULT_TEMP
800 assert T_convert_to_K == C_to_K, "You didn't define get_T()!"
801 Ts=zeros((num_Ts,), dtype=float)
802 for i in range(num_Ts) :
803 Ts[i] = T(get_T, T_convert_to_K, log_dir, plotVerbose)
804 time.sleep(1) # wait a bit to get an independent temperature measure
808 move_far_from_surface(stepper)
809 vibs=zeros((num_vibs,))
810 for i in range(num_vibs) :
811 vibs[i] = vib(zpiezo, vib_time, vib_freq, log_dir=log_dir)
814 if LOG_DATA != None :
815 data = {'bump':bumps, 'T':Ts, 'vib':vibs}
816 log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
818 log.write_dict_of_arrays(data)
820 return (bumps, Ts, vibs)
822 def calib_save(bumps, Ts, vibs, log_dir) :
824 Save a dictonary with the bump, T, and vib data.
827 data = {'bump':bumps, 'T':Ts, 'vib':vibs}
828 log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
830 log.write_dict_of_arrays(data)
832 def calib_load(datafile) :
833 "Load the dictionary data, using data_logger.date_load()"
834 dl = data_logger.data_load()
835 data = dl.read_dict_of_arrays(path)
836 return (data['bump'], data['T'], data['vib'])
838 def calib_save_analysis(k, k_s,
839 photoSensitivity2_m, photoSensitivity2_s,
840 T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s,
843 log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
844 log_name="calib_analysis_text")
846 "k (N/m) : %g +/- %g\n" % (k, k_s)
847 + "photoSensitivity**2 (V/nm)**2 : %g +/- %g\n" % \
848 (photoSensitivity2_m, photoSensitivity2_s)
849 + "T (K) : %g +/- %g\n" % (T_m, T_s)
850 + "1/Vp**2 (1/V)**2 : %g +/- %g\n" % \
851 (one_o_Vphoto2_m, one_o_Vphoto2_s)
854 def calib_plot(bumps, Ts, vibs, plotVerbose) :
855 if plotVerbose or PYLAB_VERBOSE :
857 _pylab.figure(BASE_FIGNUM+4)
859 _pylab.plot(bumps, 'g.')
860 _pylab.title('Photodiode sensitivity (V/nm)')
862 _pylab.plot(Ts, 'r.')
863 _pylab.title('Temperature (K)')
865 _pylab.plot(vibs, 'b.')
866 _pylab.title('Thermal deflection variance (Volts**2)')
869 def calib(stepper, zpiezo, tempController=None,
871 num_bumps=10, push_depth_nm=300,
876 Calibrate a cantilever in one function.
877 The I-don't-care-about-the-details black box version :p.
882 k cantilever spring constant (in N/m, or equivalently nN/nm)
883 k_s standard deviation in our estimate of k
885 See get_calibration_data() for the data aquisition code
886 See analyze_calibration_data() for the analysis code
888 a, T, vib = calib_aquire(stepper, zpiezo, tempController,
891 push_depth_nm=push_depth_nm,
895 calib_save(a, T, vib, log_dir)
896 k,k_s,ps2_m, ps2_s,T_m,T_s,one_o_Vp2_m,one_o_Vp2_s = \
897 calib_analyze(a, T, vib, log_dir=log_dir)
898 calib_save_analysis(k,k_s,ps2_m, ps2_s,T_m,T_s,one_o_Vp2_m,one_o_Vp2_s,
900 calib_plot(a, T, vib, plotVerbose)
903 def calib_load_analyze_tweaks(bump_tweaks, vib_tweaks,
907 a = read_tweaked_bumps(bump_tweaks)
908 vib = V_photo_variance_from_file(vib_tweaks)
909 return analyze_calibration_data(a, T, vib, log_dir=log_dir)