calibcant/vib_analyze.py

   1 #!/usr/bin/python
   2 #
   3 # calibcant - tools for thermally calibrating AFM cantilevers
   4 #
   5 # Copyright (C) 2007-2009 William Trevor King
   6 #
   7 # This program is free software; you can redistribute it and/or
   8 # modify it under the terms of the GNU General Public License as
   9 # published by the Free Software Foundation; either version 3 of the
  10 # License, or (at your option) any later version.
  11 #
  12 # This program is distributed in the hope that it will be useful, but
  13 # WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  15 # See the GNU General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU General Public License
  18 # along with this program; if not, write to the Free Software
  19 # Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  20 # 02111-1307, USA.
  21 #
  22 # The author may be contacted at <wking@drexel.edu> on the Internet, or
  23 # write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
  24 # Philadelphia PA 19104, USA.
  25
  26 """
  27 Separate the more general vib_analyze() from the other vib_*()
  28 functions in calibcant.  Also provide a command line interface
  29 for analyzing data acquired through other workflows.
  30
  31 The relevent physical quantities are :
  32  Vphoto   The photodiode vertical deflection voltage (what we measure)
  33 """
  34
  35 import os, time, numpy
  36 import GnuplotBiDir   # can be used for fitting lorentzian
  37 import scipy.optimize # can be used for fitting lorentzian
  38
  39 import data_logger
  40 import FFT_tools
  41 from splittable_kwargs import splittableKwargsFunction, \
  42     make_splittable_kwargs_function
  43
  44 import .common
  45 import .config
  46
  47
  48 PLOT_GUESSED_LORENTZIAN=False
  49
  50 @splittableKwargsFunction()
  51 def vib_analyze_naive(deflection_bits, Vphoto_in2V=config.Vphoto_in2V,
  52                       zeroVphoto_bits=config.zeroVphoto_bits) :
  53     """
  54     Calculate the variance of the raw data, and convert it to Volts**2.
  55     This method is simple and easy to understand, but it highly succeptible
  56     to noise, drift, etc.
  57
  58     Vphoto_in2V is a function converting Vphoto values from bits to Volts.
  59     This function is assumed linear, see bump_analyze().
  60     """
  61     if config.TEXT_VERBOSE :
  62         print "Calculating the naive (raw) variance"
  63     Vphoto_std_bits = deflection_bits.std()
  64     Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zeroVphoto_bits)
  65     if config.TEXT_VERBOSE :
  66         print "Average deflection (Volts):", Vphoto_in2V(deflection_bits.mean())
  67     return Vphoto_std**2
  68
  69 #@splittableKwargsFunction((fit_psd, 'freq_axis', 'psd_data'),
  70 #                          (vib_plot, 'deflection_bits', 'freq_axis', 'power',
  71 #                           'A', 'B', 'C', 'minFreq', 'maxFreq'))
  72 # Some of the child functions aren't yet defined, so postpone
  73 # make-splittable until later in the module.
  74 def vib_analyze(deflection_bits, freq,
  75                 chunk_size=4096, overlap=True, window=FFT_tools.window_hann,
  76                 Vphoto_in2V=config.Vphoto_in2V, **kwargs) :
  77     """
  78     Calculate the variance in the raw data due to the cantilever's
  79     thermal oscillation and convert it to Volts**2.
  80
  81     Improves on vib_analyze_naive() by first converting Vphoto(t) to
  82     frequency space, and fitting a Lorentzian in the relevant
  83     frequency range (see cantilever_calib.pdf for derivation).
  84     However, there may be cases where the fit is thrown off by noise
  85     spikes in frequency space.  To protect from errors, the fitted
  86     variance is compared to the naive variance (where all noise is
  87     included), and the minimum variance is returned.
  88
  89     Vphoto_in2V is a function converting Vphoto values from bits to Volts.
  90     This function is assumed linear, see bump_analyze().
  91     """
  92     fit_psd_kwargs,vib_plot_kwargs = vib_analyze._splitargs(vib_analyze,kwargs)
  93     if 'minFreq' in fit_psd_kwargs:
  94         vib_plot_kwargs['minFreq'] = fit_psd_kwargs['minFreq']
  95     if 'maxFreq' in fit_psd_kwargs:
  96         vib_plot_kwargs['maxFreq'] = fit_psd_kwargs['maxFreq']
  97
  98     Vphoto_t = numpy.zeros((len(deflection_bits),),
  99                            dtype=numpy.float)
 100     # convert the data from bits to volts
 101     if config.TEXT_VERBOSE :
 102         print "Converting %d bit values to voltages" % len(Vphoto_t)
 103     for i in range(len(deflection_bits)) :
 104         Vphoto_t[i] = Vphoto_in2V(deflection_bits[i])
 105     if config.TEXT_VERBOSE :
 106         print "Average deflection (Volts):", Vphoto_t.mean()
 107
 108     # Compute the average power spectral density per unit time (in uV**2/Hz)
 109     if config.TEXT_VERBOSE :
 110         print "Compute the averaged power spectral density in uV**2/Hz"
 111     (freq_axis, power) = \
 112         FFT_tools.unitary_avg_power_spectrum((Vphoto_t - Vphoto_t.mean())*1e6,
 113                                              freq, chunk_size,overlap, window)
 114
 115     A,B,C,D = fit_psd(freq_axis, power, **kwargs)
 116
 117     if config.TEXT_VERBOSE :
 118         print "Fitted f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with"
 119         print "A = %g, \t B = %g, \t C = %g, \t D = %g" % (A, B, C, D)
 120
 121     vib_plot(deflection_bits, freq_axis, power, A, B, C, D, **kwargs)
 122
 123     # Our A is in uV**2, so convert back to Volts**2
 124     fitted_variance = lorentzian_area(A,B,C) * 1e-12
 125
 126     naive_variance = vib_analyze_naive(deflection_bits,
 127                                        Vphoto_in2V=Vphoto_in2V)
 128     if config.TEXT_VERBOSE:
 129         print "Fitted variance: %g V**2" % fitted_variance
 130         print "Naive variance : %g V**2" % naive_variance
 131
 132     return min(fitted_variance, naive_variance)
 133
 134 @splittableKwargsFunction()
 135 def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=25000) :
 136     """
 137     freq_axis : array of frequencies in Hz
 138     psd_data  : array of PSD amplitudes for the frequencies in freq_axis
 139
 140     Calculate the variance in the raw data due to the cantilever's
 141     thermal oscillation and convert it to Volts**2.
 142
 143     The conversion to frequency space generates an average power spectrum
 144     by breaking the data into windowed chunks and averaging the power
 145     spectrums for the chunks together.
 146     See FFT_tools.unitary_avg_power_spectrum().
 147     """
 148     # cut out the relevent frequency range
 149     if config.TEXT_VERBOSE :
 150         print "Cut the frequency range %g to %g Hz" % (minFreq, maxFreq)
 151     imin = 0
 152     while freq_axis[imin] < minFreq : imin += 1
 153     imax = imin
 154     while freq_axis[imax] < maxFreq : imax += 1
 155     assert imax >= imin + 10 , \
 156         "Less than 10 points in freq range (%g,%g)" % (minFreq, maxFreq)
 157
 158     # generate guesses for Lorentzian parameters A,B,C
 159     max_psd_index = numpy.argmax(psd_data[imin:imax]) + imin
 160     max_psd = psd_data[max_psd_index]
 161     half_max_index = imin
 162     while psd_data[half_max_index] < max_psd/2 :
 163         half_max_index += 1
 164     res_freq = freq_axis[max_psd_index]
 165     half_width = freq_axis[max_psd_index] - freq_axis[half_max_index]
 166     # Lorentzian L(x) = C / ((A**2-x**2)**2 + (B*x)**2)
 167     # is expected power spectrum for
 168     # x'' + B x' + A^2 x'' = F_external(t)/m
 169     # (A = omega_0)
 170     # C = (2 k_B T B) / (pi m)
 171     #
 172     # A = resonant frequency
 173     # peak at  x_res = sqrt(A^2 - B^2/2)
 174     #  which could be complex if there isn't a peak (overdamped)
 175     # peak height    = C / (3 x_res^4 + A^4)
 176     # Q = A/B
 177     #
 178     # guess Q = 1, and adjust from there
 179     Q_guess = 1
 180     # so x_res^2 = B^2 Q^2 - B^2/2 = (Q^2-1/2)B^2
 181     #    B = x_res / sqrt(Q^2-1/2)
 182     if config.TEXT_VERBOSE :
 183         print "params : %g, %g" % (res_freq, max_psd)
 184     B_guess = res_freq / numpy.sqrt(Q_guess**2-0.5)
 185     A_guess = Q_guess*B_guess
 186     C_guess = max_psd * (-res_freq**4 + A_guess**4)
 187     D_guess = 0 # allow a noise floor offset (DISABLED in fitting below)
 188     #
 189     # Half width w on lower side when L(a-w) = L(a)/2
 190     #  (a**2 - (a-w)**2)**2 + (b*(a-w))**2 = 2*(b*a)**2
 191     # Let W=(a-w)**2, A=a**2, and B=b**2
 192     #  (A - W)**2 + BW = 2*AB
 193     #  W**2 - 2AW + A**2 + BW = 2AB
 194     #  W**2 + (B-2A)W + (A**2-2AB) = 0
 195     #  W = (2A-B)/2 * [1 +/- sqrt(1 - 4(A**2-2AB)/(B-2A)**2]
 196     #    = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
 197     #  (a-w)**2 = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
 198     #  so w is a disaster ;)
 199     # For some values of A and B (non-underdamped), W is imaginary or negative.
 200     #
 201     # The mass m is given by m = k_B T / (pi A)
 202     # The spring constant k is given by k = m (omega_0)**2
 203     # The quality factor Q is given by Q = omega_0 m / gamma
 204     if config.TEXT_VERBOSE :
 205         print "guesses : %g, %g, %g" % (A_guess, B_guess, C_guess)
 206     if PLOT_GUESSED_LORENTZIAN:
 207         vib_plot(None, freq_axis, psd_data, A_guess, B_guess, C_guess,
 208                  minFreq, maxFreq, plotVerbose=True)
 209
 210     # Fitting the PSD of V = photoSensitivity*x just rescales the parameters
 211
 212     if False: # Gnuplot fit worse than scipy.optimize.leastsq fit.
 213         # fit Lorentzian using Gnuplot's 'fit' command
 214
 215         # save about-to-be-fitted data to a temporary file
 216         datafilename = "%s_%d" % (config.GNUFIT_DATA_BASE, time.time())
 217         fd = file(datafilename, 'w')
 218         for i in range(imin, imax) :
 219             fd.write("%g\t%g\n" % (freq_axis[i], psd_data[i]))
 220         fd.close()
 221
 222         g = GnuplotBiDir.Gnuplot()
 223         g("f(x) = C / ((A**2-x**2)**2 + (B*x)**2)")
 224         g("A=%g; B=%g; C=%g" % (A_guess, B_guess, C_guess))
 225         g("fit f(x) '%s' via A,B,C" % datafilename)
 226         A=abs(float(g.getvar('A'))) # A and B only show up as squares in f(x)
 227         B=abs(float(g.getvar('B'))) # so ensure we get positive values
 228         C=float(g.getvar('C'))
 229
 230         os.remove(datafilename)
 231     else:
 232         # fit Lorenzian using scipy.optimize.leastsq
 233         def residual(p, y, x):
 234             A = p[0]
 235             B = p[1]
 236             C = p[2]
 237             Y = C / ((A**2-x**2)**2 + (B*x)**2)
 238             return Y-y
 239         leastsq = scipy.optimize.leastsq
 240         p,cov,info,mesg,ier = leastsq(residual, (A_guess, B_guess, C_guess),
 241                                       args=(psd_data[imin:imax],
 242                                             freq_axis[imin:imax]),
 243                                       full_output=True, maxfev=10000)
 244         if config.TEXT_VERBOSE:
 245             print "Fitted params:",p
 246             print "Covariance mx:",cov
 247             #print "Info:", info  # too much information :p
 248             print "mesg:", mesg
 249             if ier == 1 :
 250                 print "Solution converged"
 251             else :
 252                 print "Solution did not converge"
 253         A,B,C = p
 254         D = 0 # do not fit the noise floor (fit doesn't look very convincing)
 255         A=abs(A) # A and B only show up as squares in f(x)
 256         B=abs(B) # so ensure we get positive values
 257     return (A, B, C, D)
 258
 259 def lorentzian_area(A, B, C) :
 260     # Integrating the the power spectral density per unit time (power) over the
 261     # frequency axis [0, fN] returns the total signal power per unit time
 262     #  int_0^fN power(f)df = 1/T int_0^T |x(t)**2| dt
 263     #                      = <V(t)**2>, the variance for our AC signal.
 264     # The variance from our fitted Lorentzian is the area under the Lorentzian
 265     #  <V(t)**2> = (pi*C) / (2*B*A**2)
 266     return (numpy.pi * C) / (2 * B * A**2)
 267
 268 def gnuplot_check_fit(datafilename, A, B, C) :
 269     """
 270     return a string containing a gnuplot script displaying the fit.
 271     """
 272     string = 'f(x) = C / ((A**2-x**2)**2 + (x*B)**2)\n'
 273     string += 'A= %g ; B= %g ; C= %g\n' % (A, B, C)
 274     string += 'set logscale y\n'
 275     string += "plot '%s', f(x)\n" % datafilename
 276     return string
 277
 278 @splittableKwargsFunction()
 279 def vib_save(data, log_dir=None) :
 280     """Save the dictionary data, using data_logger.data_log()
 281     data should be a dictionary of arrays with the fields
 282       'Z piezo output'
 283       'Z piezo input'
 284       'Deflection input'
 285     """
 286     if log_dir :
 287         freq = data['sample frequency Hz']
 288         log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
 289                                    log_name="vib_%gHz" % freq)
 290         log.write_dict_of_arrays(data)
 291
 292 def vib_load(datafile=None) :
 293     """Load the dictionary data, using data_logger.date_load()"
 294     data should be a dictionary of arrays with the fields
 295       'Z piezo output'
 296       'Z piezo input'
 297       'Deflection input'
 298     """
 299     dl = data_logger.data_load()
 300     data = dl.read_dict_of_arrays(datafile)
 301     return data
 302
 303 @splittableKwargsFunction()
 304 def vib_plot(deflection_bits, freq_axis, power, A, B, C, D,
 305              minFreq=None, maxFreq=None, plotVerbose=False) :
 306     """
 307     If plotVerbose or config.PYLAB_VERBOSE == True, plot 3 subfigures:
 308      Time series (Vphoto vs sample index) (show any major drift events),
 309      A histogram of Vphoto, with a gaussian fit (demonstrate randomness), and
 310      FFTed Vphoto data (Vphoto vs Frequency) (show major noise components).
 311     """
 312     if plotVerbose or config.PYLAB_VERBOSE :
 313         common._import_pylab()
 314         common._pylab.figure(config.BASE_FIGNUM+2)
 315
 316         if deflection_bits != None:
 317             # plot time series
 318             common._pylab.subplot(311)
 319             common._pylab.hold(False)
 320             common._pylab.plot(deflection_bits, 'r.')
 321             common._pylab.title("free oscillation")
 322
 323             # plot histogram distribution and gaussian fit
 324             common._pylab.subplot(312)
 325             common._pylab.hold(False)
 326             n, bins, patches = \
 327                 common._pylab.hist(deflection_bits, bins=30,
 328                                    normed=1, align='center')
 329             gaus = numpy.zeros((len(bins),), dtype=numpy.float)
 330             mean = deflection_bits.mean()
 331             std = deflection_bits.std()
 332             pi = numpy.pi
 333             exp = numpy.exp
 334             for i in range(len(bins)) :
 335                 gaus[i] = (2*pi)**(-.5)/std * exp(-0.5*((bins[i]-mean)/std)**2)
 336             common._pylab.hold(True)
 337             common._pylab.plot(bins, gaus, 'r-');
 338             common._pylab.hold(False)
 339
 340             # plot FFTed data
 341             axes = common._pylab.subplot(313)
 342         else:
 343             # use a nice big subplot just for the FFTed data
 344             axes = common._pylab.subplot(111)
 345         common._pylab.hold(False)
 346         common._pylab.semilogy(freq_axis, power, 'r.-')
 347         common._pylab.hold(True)
 348         xmin,xmax = axes.get_xbound()
 349         ymin,ymax = axes.get_ybound()
 350
 351         if minFreq is not None and maxFreq is not None:
 352             # highlight the region we're fitting in
 353             common._pylab.axvspan(minFreq, maxFreq, facecolor='g', alpha=0.1,
 354                                   zorder = -1)
 355
 356         fitdata = C / ((A**2-freq_axis**2)**2 + (B*freq_axis)**2) + D
 357         common._pylab.plot(freq_axis, fitdata, 'b-')
 358         noisefloor = D + 0*freq_axis;
 359         common._pylab.plot(freq_axis, noisefloor)
 360         common._pylab.hold(False)
 361         axes.axis([xmin,xmax,ymin,ymax])
 362
 363         common._flush_plot()
 364     if (plotVerbose or config.GNUPLOT_VERBOSE) and False: # TODO: cleanup and test
 365         # write all the ft data now
 366         fd = file(datafilename, 'w')
 367         for i in range(len(freq_axis)) :
 368             fd.write("%g\t%g\n" % (freq_axis[i], power[i]))
 369         fd.write("\n") # blank line to separate fit data.
 370         # write the fit data
 371         for i in range(imin, imax) :
 372             x = freq_axis[i]
 373             fd.write("%g\t%g\n" % (freq_axis[i],
 374                                    C / ((A**2-x**2)**2 + (x*B)**2) ) )
 375         fd.close()
 376         print Vphoto_t.var(), A, B, C, ((numpy.pi * C)/(2 * B * A**2) * 1e-12)
 377         g("set terminal epslatex color solid")
 378         g("set output 'calibration.tex'")
 379         g("set size 2,2") # multipliers on default 5"x3.5"
 380         #g("set title 'Thermal calibration'")
 381         g("set logscale y")
 382         g("set xlabel 'Frequency (Hz)'")
 383         g("set ylabel 'Power ($\mu$V$^2$/Hz)'")
 384         g("set xrange [0:20000]")
 385         g("g(x) = x < %g || x > %g ? 1/0 : f(x)" % (minFreq, maxFreq)) # only plot on fitted region
 386         g("plot '%s' title 'Vibration data', g(x) title 'Lorentzian fit'" % datafilename)
 387
 388         g("set mouse")
 389         g("pause mouse 'click with mouse'")
 390         g.getvar("MOUSE_BUTTON")
 391         del(g)
 392
 393 # Make vib_analyze splittable (was postponed until the child functions
 394 # were defined).
 395 make_splittable_kwargs_function(vib_analyze,
 396     (fit_psd, 'freq_axis', 'psd_data'),
 397     (vib_plot, 'deflection_bits', 'freq_axis', 'power', 'A', 'B', 'C', 'D',
 398      'minFreq', 'maxFreq'))
 399
 400 @splittableKwargsFunction((vib_analyze, 'deflection_bits', 'freq'))
 401 def vib_load_analyze_tweaked(tweak_file, analyze=vib_analyze, **kwargs) :
 402     """
 403     Read the vib files listed in tweak_file.
 404     Return an array of Vphoto variances (due to the cantilever) in Volts**2
 405     """
 406     analyze_kwargs, = \
 407         vib_load_analyze_tweaked._splitargs(vib_load_analyze_tweaked, kwargs)
 408     dl = data_logger.data_load()
 409     Vphoto_var = []
 410     for path in file(tweak_file, 'r') :
 411         path = path.strip()
 412         if path[0] == '#': # a comment
 413             continue
 414         # read the data
 415         data = vib_load(path)
 416         deflection_bits = data['Deflection input']
 417         freq = float(path.split('_')[-1].split('Hz')[0])
 418         if config.TEXT_VERBOSE :
 419             print "Analyzing '%s' at %g Hz" % (path, freq)
 420         # analyze
 421         if 'freq' in analyze._kwargs(analyze):
 422             var = analyze(deflection_bits, freq, **analyze_kwargs)
 423         else:
 424             var = analyze(deflection_bits, **analyze_kwargs)
 425         Vphoto_var.append(var)
 426     return numpy.array(Vphoto_var, dtype=numpy.float)
 427
 428 # commandline interface functions
 429 import scipy.io, sys
 430
 431 def read_data(ifile):
 432     """
 433     ifile can be a filename string or open (seekable) file object.
 434     returns an array of data values (1 column)
 435     """
 436     #returns (column 1 array, column 2 array)
 437     #"""
 438     if ifile == None :  ifile = sys.stdin
 439     unlabeled_data=scipy.io.read_array(ifile)
 440     return unlabeled_data #(unlabeled_data[:,0], unlabeled_data[:,1])
 441
 442 if __name__ == '__main__' :
 443     # command line interface
 444     from optparse import OptionParser
 445
 446     usage_string = ('%prog <input-file>\n'
 447                     '2007-2009 W. Trevor King.\n'
 448                     '\n'
 449                     'There are several operation modes, each handling a different <input-file>\n'
 450                     'type.  In each case, the output is the fitted variance (in square Volts).\n'
 451                     '\n'
 452                     'Single file mode without datalogger mode (the default) :\n'
 453                     '<input-file> should be a 1 column ASCII without a header line of cantilever\n'
 454                     'deflection (in bits)\n'
 455                     'e.g: "<deflection bits>\\n..."\n'
 456                     '\n'
 457                     'Single file mode with datalogger mode :\n'
 458                     'Same as without datalogger mode, except the input should be a datalogger file\n'
 459                     'with data stored in bits (e.g. as saved by the unfold module).\n'
 460                     '\n'
 461                     'Tweak file mode:\n'
 462                     'Runs the same analysis as in single file mode for each vib file in a tweak\n'
 463                     'file.  Each line in the tweak file specifies a single vib file.  Blank lines\n'
 464                     'and those beginning with a pound sign (#) are ignored.  Vib files are valid\n'
 465                     'datalogger files as per single-file-with-datalogger-mode.\n'
 466                     )
 467     parser = OptionParser(usage=usage_string, version='%prog 0.1')
 468     parser.add_option('-f', '--sample-freq', dest='freq', default=100e3,
 469                       help='sample frequency in Hz (default %default)',
 470                       type='float', metavar='FREQ')
 471     parser.add_option('-m', '--min-freq', dest='min_freq', default=500,
 472                       help='minimum frequency in Hz (default %default)',
 473                       type='float', metavar='FREQ')
 474     parser.add_option('-M', '--max-freq', dest='max_freq', default=7000,
 475                       help='maximum frequency in Hz (default %default)',
 476                       type='float', metavar='FREQ')
 477     parser.add_option('-o', '--output-file', dest='ofilename',
 478                       help='write output to FILE (default stdout)',
 479                       type='string', metavar='FILE')
 480     parser.add_option('-c', '--comma-out', dest='comma_out', action='store_true',
 481                       help='Output comma-seperated values (default %default)',
 482                       default=False)
 483     parser.add_option('-g', '--gnuplot', dest='gnuplot', action='store_true',
 484                       help='Print gnuplot fit check script to stderr',
 485                       default=False)
 486     parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
 487                       help='Produce pylab fit checks during execution',
 488                       default=False)
 489     parser.add_option('-G', '--plot-guess', dest='plot_guess', action='store_true',
 490                       help='Produce plots of the pre-fit Lorentzian',
 491                       default=False)
 492     parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
 493                       help='Run in tweak-file mode',
 494                       default=False)
 495     parser.add_option('-d', '--datalogger-mode', dest='datalogger_mode', action='store_true',
 496                       help='Read input files with datalogger.read_dict_of_arrays().  This is useful, for example, to test a single line from a tweakfile.',
 497                       default=False)
 498     parser.add_option('-s', '--disable-spectrum-fitting', dest='spectrum_fitting',
 499                       action='store_false',
 500                       help='Disable PSD fitting, just use the raw variance',
 501                       default=True)
 502     parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
 503                       help='Print lots of debugging information',
 504                       default=False)
 505
 506     options,args = parser.parse_args()
 507     parser.destroy()
 508     assert len(args) >= 1, "Need an input file"
 509
 510     ifilename = args[0]
 511
 512     if options.ofilename != None :
 513         ofile = file(options.ofilename, 'w')
 514     else :
 515         ofile = sys.stdout
 516     config.TEXT_VERBOSE = options.verbose
 517     config.PYLAB_INTERACTIVE = False
 518     config.PYLAB_VERBOSE = options.pylab
 519     config.GNUPLOT_VERBOSE = options.gnuplot
 520     PLOT_GUESSED_LORENTZIAN = options.plot_guess
 521     if options.spectrum_fitting == True:
 522         analyze = vib_analyze
 523         analyze_args = {'freq':options.freq,
 524                         'minFreq':options.min_freq,
 525                         'maxFreq':options.max_freq}
 526     else:
 527         analyze = vib_analyze_naive
 528         analyze_args = dict()
 529         make_splittable_kwargs_function(vib_load_analyze_tweaked,
 530                                         (vib_analyze, 'deflection_bits'))
 531
 532     if options.tweakmode == False : # single file mode
 533         if options.datalogger_mode:
 534             data = vib_load(ifilename)
 535             deflection_bits = data['Deflection input']
 536         else:
 537             deflection_bits = read_data(ifilename)
 538         Vphoto_var = analyze(deflection_bits, **analyze_args)
 539         print >> ofile, Vphoto_var
 540     else : # tweak mode
 541         if 'freq' in analyze_args:
 542             analyze_args.pop('freq') # freq extracted from vib file names
 543         Vphoto_vars = vib_load_analyze_tweaked(ifilename, analyze=analyze,
 544                                                **analyze_args)
 545         if options.comma_out :
 546             sep = ','
 547         else :
 548             sep = '\n'
 549         common.write_array(ofile, Vphoto_vars, sep)
 550
 551     if common._final_flush_plot != None:
 552         common._final_flush_plot()
 553
 554     if options.ofilename != None :
 555         ofile.close()
 556
 557 #  LocalWords:  calibcant AFM