calibcant/calibrate_cantilever.py

   1 #!/usr/bin/python
   2
   3 """
   4 Aquire and analyze cantilever calibration data.
   5
   6 W. Trevor King Dec. 2007-Jan. 2008
   7
   8 GPL BOILERPLATE
   9
  10
  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
  21
  22 Which are related by the parameters :
  23  zpGain           Vzp_out / Vzp
  24  zpSensitivity    Zp / Vzp
  25  photoSensitivity Vphoto / Zcant
  26  k_cant           Fcant / Zcant
  27
  28 Cantilever calibration assumes a pre-calibrated z-piezo (zpSensitivity) and
  29 a amplifier (zpGain).
  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
  32 on the Nanoscope.
  33
  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)
  39
  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.
  54
  55 We do all these measurements a few times to estimate statistical errors.
  56
  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
  66  load_analyze_tweaked
  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().
  73
  74 Also defines the two positioning functions:
  75  move_just_onto_surface() and move_far_from_surface()
  76 """
  77
  78 import numpy
  79 import time
  80 import data_logger
  81 import z_piezo_utils
  82 import FFT_tools
  83 import linfit
  84 import GnuplotBiDir # used for fitting lorentzian
  85
  86 kb = 1.3806504e-23 # Boltzmann's constant in J/K
  87 DEFAULT_TEMP = 22  # assume 22 deg C
  88
  89 LOG_DATA = True  # quietly grab all real-world data and log to LOG_DIR
  90 LOG_DIR = '/home/wking/rsrch/data/calibrate_cantilever'
  91
  92 GNUFIT_DATA_BASE='./calibrate_cantilever_fitdata'
  93
  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
  99
 100 # handle extra verbose input modules, only imported if we need them
 101 _flush_plot = None
 102 _final_flush_plot = None
 103 _pylab = None
 104 def _import_pylab() :
 105     "Import pylab plotting functions for when we need to plot"
 106     global _pylab
 107     global _flush_plot
 108     global _final_flush_plot
 109     if _pylab == None :
 110         import pylab as _pylab
 111     if _flush_plot == None :
 112         if PYLAB_INTERACTIVE :
 113             _flush_plot = _pylab.draw
 114         else :
 115             def _flush_plot () : pass
 116     if _final_flush_plot == None :
 117         if PYLAB_INTERACTIVE :
 118             _final_flush_plot = _pylab.draw
 119         else :
 120             _final_flush_plot = _pylab.show
 121
 122
 123 # HACK
 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
 126 _usual_zpGain=20
 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)
 134
 135 def C_to_K(celsius) :
 136     "Convert Celsius -> Kelvin."
 137     return celsius + 273.15
 138
 139 # bump family
 140
 141 def bump_aquire(zpiezo, push_depth, npoints, freq) :
 142     """
 143     Ramps closer push_depth and returns to the original position.
 144     Inputs:
 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.
 152     """
 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
 161     if TEXT_VERBOSE :
 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)
 169     return data
 170
 171 def bump_save(data, log_dir) :
 172     "Save the dictionary data, using data_logger.data_log()"
 173     if log_dir != None :
 174         log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
 175                                    log_name="bump")
 176         log.write_dict_of_arrays(data)
 177
 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)
 182     return data
 183
 184 def bump_analyze(data, zpGain=_usual_zpGain,
 185                  zpSensitivity=_usual_zpSensitivity,
 186                  Vzp_out2V=_usual_Vzp_out2V,
 187                  Vphoto_in2V=_usual_Vphoto_in2V,
 188                  plotVerbose=False) :
 189     """
 190     Return the slope of the bump ;).
 191     Inputs:
 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
 197     Returns:
 198      photoSensitivity (Vphoto/Zcant) in Volts/nm
 199     Checks for strong correlation (r-value) and low randomness chance (p-value)
 200
 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.
 204     """
 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
 212
 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
 216
 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 :
 220         _import_pylab()
 221         _pylab.figure(BASE_FIGNUM)
 222         _pylab.plot(data["Z piezo output"], data["Deflection input"],
 223                     '.', label='bump')
 224         _pylab.title("bump surface")
 225         _pylab.legend(loc='upper left')
 226         _flush_plot()
 227
 228 def bump(zpiezo, push_depth, npoints=1024, freq=100e3,
 229          log_dir=None,
 230          plotVerbose=False) :
 231     """
 232     Wrapper around bump_aquire(), bump_save(), bump_analyze(), bump_plot()
 233     """
 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
 240
 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,
 245                               plotVerbose=False) :
 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]
 251         # read the data
 252         full_data = bump_load(path)
 253         if len(parsed) == 1 :
 254             data = full_data # use whole bump
 255         else :
 256             # use the listed sections
 257             zp = []
 258             df = []
 259             for rng in parsed[1:] :
 260                 p = rng.split(':')
 261                 starti = int(p[0])
 262                 stopi = int(p[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)
 272
 273 # T family.
 274 # Fairly stubby, since a one shot Temp measurement is a common thing.
 275 # We just wrap that to provide a consistent interface.
 276
 277 def T_aquire(get_T=None) :
 278     """
 279     Measure the current temperature of the sample,
 280     or, if get_T == None, fake it by returning DEFAULT_TEMP
 281     """
 282     if get_T == None :
 283         if TEXT_VERBOSE :
 284             print "Fake temperature %g" % DEFAULT_TEMP
 285         return DEFAULT_TEMP
 286     else :
 287         if TEXT_VERBOSE :
 288             print "Measure temperature"
 289         return get_T()
 290
 291 def T_save(T, log_dir=None) :
 292     """
 293     Save either a single T (if you are crazy :p),
 294     or an array of them (more reasonable)
 295     """
 296     T = numpy.array(T, dtype=numpy.float)
 297     if log_dir != None :
 298         log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
 299                                    log_name="T_float")
 300         log.write_binary(T.tostring())
 301     if LOG_DATA != None :
 302         log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
 303                                    log_name="T_float")
 304         log.write_binary(T.tostring())
 305
 306 def T_load(datafile=None) :
 307     """
 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.
 310     """
 311     if datafile == None :
 312         return numpy.array([DEFAULT_TEMP], dtype=numpy.float)
 313     else :
 314         return fromfile(file=cleanpath, dtype=numpy.float)
 315
 316 def T_analyze(T, convert_to_K=C_to_K) :
 317     """
 318     Not much to do here, just convert to Kelvin.
 319     Uses C_to_K (defined above) by default.
 320     """
 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])
 324     except TypeError :
 325         T = convert_to_K(T)
 326     return T
 327
 328 def T_plot(T, plotVerbose) :
 329     """
 330     Umm, just print the temperature?
 331     """
 332     if plotVerbose or PYLAB_VERBOSE or TEXT_VERBOSE :
 333         print "Temperature ", T
 334
 335 def T(get_T=None, convert_to_K=C_to_K, log_dir=None, plotVerbose=None) :
 336     """
 337     Wrapper around T_aquire(), T_save(), T_analyze(), T_plot()
 338     """
 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)
 343     return T_ret
 344
 345 def T_load_analyze_tweaked(tweak_file=None, convert_to_K=C_to_K,
 346                            plotVerbose=False) :
 347     "Read the T files listed in tweak_file, return an array of Ts in Kelvin"
 348     if tweak_file != None :
 349         Tlist=[]
 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)
 355         for Ti in T_raw :
 356             T_plot(T, plotVerbose)
 357     else :
 358         T_raw = T_load(None)
 359     T = T_analyze(T_raw, convert_to_K)
 360     return T
 361
 362 # vib family
 363
 364 def vib_aquire(zpiezo, time, freq) :
 365     """
 366     Record data for TIME seconds at FREQ Hz from ZPIEZO at it's current position.
 367     """
 368     # round up to the nearest power of two, for efficient FFT-ing
 369     nsamps = ceil_pow_of_two(time*freq)
 370     time = nsamps / freq
 371     # take some data, keeping the position voltage at it's current value
 372     out = ones((nsamps,), dtype=uint16) * zpiezo.curPos()
 373     if TEXT_VERBOSE :
 374         print "get %g seconds of data" % time
 375     data = zpiezo.ramp(out, freq)
 376     if LOG_DATA :
 377         log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
 378                                    log_name="vib_%gHz" % freq)
 379         log.write_dict_of_arrays(data)
 380     return data
 381
 382 def vib_save(data, log_dir=None) :
 383     "Save the dictionary data, using data_logger.data_log()"
 384     if log_dir :
 385         log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
 386                                    log_name="vib_%gHz" % freq)
 387         log.write_dict_of_arrays(data)
 388
 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)
 393     return data
 394
 395 def vib_plot(data, freq,
 396              chunk_size=2048, overlap=True, window=FFT_tools.window_hann,
 397              plotVerbose=True) :
 398     """
 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).
 403     """
 404     if plotVerbose or PYLAB_VERBOSE :
 405         _import_pylab()
 406         _pylab.hold(False)
 407         _pylab.figure(BASE_FIGNUM+2)
 408
 409         # plot time series
 410         _pylab.subplot(311)
 411         _pylab.plot(data["Deflection input"], 'r.')
 412         _pylab.title("free oscillation")
 413
 414         # plot histogram distribution and gaussian fit
 415         _pylab.subplot(312)
 416         n, bins, patches = \
 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()
 422         pi = numpy.pi
 423         exp = numpy.exp
 424         for i in range(len(bins)) :
 425             gaus[i] = (2*pi)**(-.5)/std * exp(-0.5*((bins[i]-mean)/std)**2)
 426         _pylab.hold(True)
 427         _pylab.plot(bins, gaus, 'r-');
 428         _pylab.hold(False)
 429
 430         # plot FFTed data
 431         _pylab.subplot(313)
 432         AC_data = data["Deflection input"] - mean
 433         (freq_axis, power) = \
 434             FFT_tools.unitary_avg_power_spectrum(AC_data, freq, chunk_size,
 435                                                  overlap, window)
 436         _pylab.semilogy(freq_axis, power, 'r.-')
 437         _pylab.hold(True)
 438         filtered_power = FFT_tools.windowed_filter(power,
 439                             FFT_tools.windowed_median_point, s=5)
 440         _pylab.semilogy(freq_axis, filtered_power, 'b.-')
 441         _flush_plot()
 442
 443 def vib_analyze_naive(data, Vphoto_in2V=_usual_Vphoto_in2V,
 444                       zeroVphoto_bits=_usual_zeroVphoto_bits,
 445                       plotVerbose=False) :
 446     """
 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.
 450
 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.
 454     """
 455     Vphoto_std_bits = data["Deflection input"].std()
 456     Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zpiezo_zeroV_in)
 457     return Vphoto_std**2
 458
 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) :
 463     """
 464     Calculate the variance in the raw data due to the cantilever's
 465     thermal oscillation and convert it to Volts**2.
 466
 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).
 470
 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
 475
 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.
 479     """
 480     Vphoto_t = numpy.zeros((len(data["Deflection input"]),),
 481                            dtype=numpy.float)
 482     # convert the data from bits to volts
 483     if TEXT_VERBOSE :
 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])
 487
 488     # Compute the average power spectral density per unit time (in uV**2/Hz)
 489     if TEXT_VERBOSE :
 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)
 494
 495     # cut out the relevent frequency range
 496     if TEXT_VERBOSE :
 497         print "Cut the frequency range %g to %g Hz" % (minFreq, maxFreq)
 498     imin = 0
 499     while freq_axis[imin] < minFreq : imin += 1
 500     imax = imin
 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)
 504
 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]
 511
 512     # check
 513     if (plotVerbose or PYLAB_VERBOSE) and False :
 514         if TEXT_VERBOSE :
 515             print "Plot the FFT data for checking"
 516         vib_plot(data, freq, chunk_size, overlap, window, plotVerbose)
 517         _import_pylab()
 518         _pylab.figure(5)
 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.-')
 524         _flush_plot()
 525
 526     # save about-to-be-fitted data to a temporary file
 527     if TEXT_VERBOSE :
 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]))
 533     fd.close()
 534
 535     # generate guesses for Lorentzian parameters A,B,C
 536     if TEXT_VERBOSE :
 537         print "Fit lorentzian"
 538     max_fp_index = numpy.argmin(filtered_power)
 539     max_fp = filtered_power[max_fp_index]
 540     half_max_index = 0
 541     while filtered_power[half_max_index] < max_fp/2 :
 542         half_max_index += 1
 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.
 557     B_guess = A_guess/2
 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
 560
 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'))
 579     if TEXT_VERBOSE :
 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 :
 583         if TEXT_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.
 591         # write the fit data
 592         for i in range(imin, imax) :
 593             x = freq_axis[i]
 594             fd.write("%g\t%g\n" % (freq_axis[i],
 595                                    C / ((A**2-x**2)**2 + (x*B)**2) ) )
 596         fd.close()
 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'")
 602         g("set logscale y")
 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)
 608
 609         g("set mouse")
 610         g("pause mouse 'click with mouse'")
 611         g.getvar("MOUSE_BUTTON")
 612         del(g)
 613
 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
 622
 623 def vib(zpiezo, time=1, freq=50e3, log_dir=None, plotVerbose=False) :
 624     """
 625     Wrapper around vib_aquire(), vib_save(), vib_analyze(), vib_plot()
 626     """
 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)
 631     return Vphoto_var
 632
 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,
 638                              plotVerbose=False) :
 639     """
 640     Read the vib files listed in tweak_file.
 641     Return an array of Vphoto variances (due to the cantilever) in Volts**2
 642     """
 643     dl = data_logger.data_load()
 644     Vphoto_var = []
 645     for path in file(tweak_file, 'r') :
 646         if path[-1] == '\n':
 647             path = path.split('\n')[0]
 648         # read the data
 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)
 653         # analyze
 654         Vphoto_var.append(vib_analyze(data, freq, minFreq, maxFreq,
 655                                       chunk_size, overlap, window,
 656                                       median_filter_width,
 657                                       Vphoto_in2V, plotVerbose))
 658     return numpy.array(Vphoto_var, dtype=numpy.float)
 659
 660
 661 # A few positioning functions, so we can run bump_aquire() and vib_aquire()
 662 # with proper spacing relative to the surface.
 663
 664 def move_just_onto_surface(stepper, zpiezo, Depth_nm=100, setpoint=2) :
 665     """
 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
 668     of the surface.
 669     Then set Vzp to place the cantilever Depth_nm onto the surface.
 670
 671     If getSurfPos() fails to find the surface, backs off (for safety)
 672     and steps in (without moving the zpiezo) until Vphoto > setpoint.
 673     """
 674     stepper_tol = 250 # nm, generous estimate of the fullstep stepsize
 675
 676     if TEXT_VERBOSE :
 677         print "moving just onto surface"
 678     # Zero the piezo
 679     if TEXT_VERBOSE :
 680         print "zero the z piezo output"
 681     zpiezo.jumpToPos(zpiezo.pos_nm2out(0))
 682     # See if we're near the surface already
 683     if TEXT_VERBOSE :
 684         print "See if we're starting near the surface"
 685     try :
 686         dist = zpiezo.pos_out2nm( \
 687             z_piezo_utils.getSurfPos(zpiezo, zpiezo.def_V2in(setpoint))
 688                                 )
 689     except (z_piezo_utils.tooClose, z_piezo_utils.poorFit), string :
 690         if TEXT_VERBOSE :
 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
 699         if TEXT_VERBOSE :
 700             print "Single stepping approach"
 701         while cd < sp :
 702             if TEXT_VERBOSE :
 703                 print "deflection %g < setpoint %g.  step closer" % (cd, sp)
 704             stepper.step_rel(2) # Full step in
 705             zpiezo.updateInputs()
 706             cd = zpiezo.curDef()
 707         # Back off two steps (protecting against backlash)
 708         if TEXT_VERBOSE :
 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()
 714         if TEXT_VERBOSE :
 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
 717             try :
 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)
 722                 continue
 723             break
 724         if i >= 19 :
 725             print "tried %d times, still too close! bailing" % i
 726             print "probably an invalid setpoint."
 727             raise Exception, "weirdness"
 728     if TEXT_VERBOSE :
 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)))
 735         if TEXT_VERBOSE :
 736             print 'distance to surface ', dist, ' nm, step back'
 737     while dist > stepper_tol : # and step forward if we need to
 738         stepper.step_rel(2)
 739         dist = zpiezo.pos_out2nm(getSurfPos(zpiezo, zpiezo.def_V2in(setpoint)))
 740         if TEXT_VERBOSE :
 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))
 745     # and we're there :)
 746     if TEXT_VERBOSE :
 747         print "We're %g nm into the surface" % Depth_nm
 748
 749 def move_far_from_surface(stepper, um_back=50) :
 750     """
 751     Step back a specified number of microns.
 752     (uses very rough estimate of step distance at the moment)
 753     """
 754     step_nm = 100
 755     steps = int(um_back*1000/step_nm)
 756     print "step back %d steps" % steps
 757     stepper.step_rel(-steps)
 758
 759
 760 # and finally, the calib family
 761
 762 def calib_aquire(stepper, zpiezo, get_T=None,
 763                  approach_setpoint=2,
 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) :
 770     """
 771     Aquire data for calibrating a cantilever in one function.
 772     return (bump, T, vib), each of which is an array.
 773     Inputs :
 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
 786     """
 787     # get bumps
 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)
 794     if TEXT_VERBOSE :
 795         print bumps
 796
 797     # get Ts
 798     if get_T == None :
 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
 805     print Ts
 806
 807     # get vibs
 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)
 812     print vibs
 813
 814     if LOG_DATA != None :
 815         data = {'bump':bumps, 'T':Ts, 'vib':vibs}
 816         log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
 817                                    log_name="calib")
 818         log.write_dict_of_arrays(data)
 819
 820     return (bumps, Ts, vibs)
 821
 822 def calib_save(bumps, Ts, vibs, log_dir) :
 823     """
 824     Save a dictonary with the bump, T, and vib data.
 825     """
 826     if log_dir != None :
 827         data = {'bump':bumps, 'T':Ts, 'vib':vibs}
 828         log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
 829                                    log_name="calib")
 830         log.write_dict_of_arrays(data)
 831
 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'])
 837
 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,
 841                         log_dir=None) :
 842     if log_dir != None :
 843         log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
 844                                    log_name="calib_analysis_text")
 845         log.write_binary((
 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)
 852                          ))
 853
 854 def calib_plot(bumps, Ts, vibs, plotVerbose) :
 855     if plotVerbose or PYLAB_VERBOSE :
 856         _import_pylab()
 857         _pylab.figure(BASE_FIGNUM+4)
 858         _pylab.subplot(311)
 859         _pylab.plot(bumps, 'g.')
 860         _pylab.title('Photodiode sensitivity (V/nm)')
 861         _pylab.subplot(312)
 862         _pylab.plot(Ts, 'r.')
 863         _pylab.title('Temperature (K)')
 864         _pylab.subplot(313)
 865         _pylab.plot(vibs, 'b.')
 866         _pylab.title('Thermal deflection variance (Volts**2)')
 867         _flush_plot()
 868
 869 def calib(stepper, zpiezo, tempController=None,
 870           setpoint=2.0,
 871           num_bumps=10, push_depth_nm=300,
 872           num_temps=10,
 873           num_vibs=10,
 874           log_dir=None) :
 875     """
 876     Calibrate a cantilever in one function.
 877     The I-don't-care-about-the-details black box version :p.
 878     return (k, k_s)
 879     Inputs:
 880      (see calib_aquire())
 881     Outputs :
 882      k    cantilever spring constant (in N/m, or equivalently nN/nm)
 883      k_s  standard deviation in our estimate of k
 884     Notes :
 885      See get_calibration_data() for the data aquisition code
 886      See analyze_calibration_data() for the analysis code
 887     """
 888     a, T, vib = calib_aquire(stepper, zpiezo, tempController,
 889                              setpoint,
 890                              num_bumps=num_bumps,
 891                              push_depth_nm=push_depth_nm,
 892                              num_temps=num_temps,
 893                              num_vibs=num_vibs,
 894                              log_dir=log_dir)
 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,
 899                         log_dir)
 900     calib_plot(a, T, vib, plotVerbose)
 901     return (k, k_s)
 902
 903 def calib_load_analyze_tweaks(bump_tweaks, vib_tweaks,
 904                               T_tweaks=None,
 905                               log_dir=None,
 906                               plotVerbose=True) :
 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)
 910
 911