calibcant/vib_analyze.py

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