Added guessed Lorentzian plot option to vib_analyze.

This explains why my Lorentzians are occasionally too sharp...
For example, consider the 2nd and 4th fit in
  cd ~/rsrch/data/z_piezo_calib/
  cc_vib_analyze.py  -ct 20080919_cantA_tweak.vib -pG

(I know, I know, should be rsrch/data/calibcant...)

diff --git a/calibcant/vib_analyze.py b/calibcant/vib_analyze.py
index dc4af3b494684117127e619e3098ddbefd87024f..6fb16154b7766ad9d6c997437df656f02ad95e7e 100755 (executable)
--- a/calibcant/vib_analyze.py
+++ b/calibcant/vib_analyze.py
@@ -39,6 +39,8 @@ import config # config module for the calibcant package
 import data_logger
 import FFT_tools
 
 import data_logger
 import FFT_tools
 

+PLOT_GUESSED_LORENTZIAN=False
+

 def vib_analyze_naive(deflection_bits, Vphoto_in2V=config.Vphoto_in2V,
                       zeroVphoto_bits=config.zeroVphoto_bits) :
     """
 def vib_analyze_naive(deflection_bits, Vphoto_in2V=config.Vphoto_in2V,
                       zeroVphoto_bits=config.zeroVphoto_bits) :
     """


@@ -178,6 +180,9 @@ def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) :
     # For some values of A and B (non-underdamped), W is imaginary or negative.
     if config.TEXT_VERBOSE : 
         print "guesses : %g, %g, %g" % (A_guess, B_guess, C_guess)
     # For some values of A and B (non-underdamped), W is imaginary or negative.
     if config.TEXT_VERBOSE : 
         print "guesses : %g, %g, %g" % (A_guess, B_guess, C_guess)

+    if PLOT_GUESSED_LORENTZIAN:
+        vib_plot(None, freq_axis, psd_data, A_guess, B_guess, C_guess,
+                 plotVerbose=True)

 
     # fit Lorentzian using Gnuplot's 'fit' command
     g = GnuplotBiDir.Gnuplot()
 
     # fit Lorentzian using Gnuplot's 'fit' command
     g = GnuplotBiDir.Gnuplot()


@@ -251,35 +256,38 @@ def vib_plot(deflection_bits, freq_axis, power, A, B, C,
      FFTed Vphoto data (Vphoto vs Frequency) (show major noise components).
     """
     if plotVerbose or config.PYLAB_VERBOSE :
      FFTed Vphoto data (Vphoto vs Frequency) (show major noise components).
     """
     if plotVerbose or config.PYLAB_VERBOSE :

-        print "plotting"

         common._import_pylab()
         common._pylab.figure(config.BASE_FIGNUM+2)
 
         common._import_pylab()
         common._pylab.figure(config.BASE_FIGNUM+2)
 

-        # plot time series
-        common._pylab.subplot(311)
-        common._pylab.hold(False)
-        common._pylab.plot(deflection_bits, 'r.')
-        common._pylab.title("free oscillation")
-
-        # plot histogram distribution and gaussian fit
-        common._pylab.subplot(312)
-        common._pylab.hold(False)
-        n, bins, patches = \
-            common._pylab.hist(deflection_bits, bins=30,
-                               normed=1, align='center')
-        gaus = numpy.zeros((len(bins),), dtype=numpy.float)
-        mean = deflection_bits.mean()
-        std = deflection_bits.std()
-        pi = numpy.pi
-        exp = numpy.exp
-        for i in range(len(bins)) :
-            gaus[i] = (2*pi)**(-.5)/std * exp(-0.5*((bins[i]-mean)/std)**2)
-        common._pylab.hold(True)
-        common._pylab.plot(bins, gaus, 'r-');
-        common._pylab.hold(False)
-
-        # plot FFTed data
-        axes = common._pylab.subplot(313)
+        if deflection_bits != None:
+            # plot time series
+            common._pylab.subplot(311)
+            common._pylab.hold(False)
+            common._pylab.plot(deflection_bits, 'r.')
+            common._pylab.title("free oscillation")
+    
+            # plot histogram distribution and gaussian fit
+            common._pylab.subplot(312)
+            common._pylab.hold(False)
+            n, bins, patches = \
+                common._pylab.hist(deflection_bits, bins=30,
+                                   normed=1, align='center')
+            gaus = numpy.zeros((len(bins),), dtype=numpy.float)
+            mean = deflection_bits.mean()
+            std = deflection_bits.std()
+            pi = numpy.pi
+            exp = numpy.exp
+            for i in range(len(bins)) :
+                gaus[i] = (2*pi)**(-.5)/std * exp(-0.5*((bins[i]-mean)/std)**2)
+            common._pylab.hold(True)
+            common._pylab.plot(bins, gaus, 'r-');
+            common._pylab.hold(False)
+    
+            # plot FFTed data
+            axes = common._pylab.subplot(313)
+        else:
+            # use a nice big subplot just for the FFTed data
+            axes = common._pylab.subplot(111)

         common._pylab.hold(False)
         common._pylab.semilogy(freq_axis, power, 'r.-')
         fitdata = C / ((A**2-freq_axis**2)**2 + (B*freq_axis)**2)
         common._pylab.hold(False)
         common._pylab.semilogy(freq_axis, power, 'r.-')
         fitdata = C / ((A**2-freq_axis**2)**2 + (B*freq_axis)**2)


@@ -291,7 +299,7 @@ def vib_plot(deflection_bits, freq_axis, power, A, B, C,
         axes.axis([xmin,xmax,ymin,ymax])
 
         common._flush_plot()
         axes.axis([xmin,xmax,ymin,ymax])
 
         common._flush_plot()

-    if (plotVerbose or config.GNUPLOT_VERBOSE): # TODO: cleanup and test
+    if (plotVerbose or config.GNUPLOT_VERBOSE) and False: # TODO: cleanup and test

         # write all the ft data now
         fd = file(datafilename, 'w')
         for i in range(len(freq_axis)) :
         # write all the ft data now
         fd = file(datafilename, 'w')
         for i in range(len(freq_axis)) :


@@ -394,6 +402,9 @@ if __name__ == '__main__' :
     parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
                       help='Produce pylab fit checks during execution',
                       default=False)
     parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
                       help='Produce pylab fit checks during execution',
                       default=False)

+    parser.add_option('-G', '--plot-guess', dest='plot_guess', action='store_true',
+                      help='Produce plots of the pre-fit Lorentzian',
+                      default=False)

     parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
                       help='Run in tweak-file mode',
                       default=False)
     parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
                       help='Run in tweak-file mode',
                       default=False)


@@ -419,6 +430,7 @@ if __name__ == '__main__' :
     config.PYLAB_INTERACTIVE = False
     config.PYLAB_VERBOSE = options.pylab
     config.GNUPLOT_VERBOSE = options.gnuplot
     config.PYLAB_INTERACTIVE = False
     config.PYLAB_VERBOSE = options.pylab
     config.GNUPLOT_VERBOSE = options.gnuplot

+    PLOT_GUESSED_LORENTZIAN = options.plot_guess

 
     if options.tweakmode == False :
         data = read_data(ifilename)
 
     if options.tweakmode == False :
         data = read_data(ifilename)