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