Add vibration_analyze.breit_wigner_resonant_frequency().

[calibcant.git] / calibcant / vibration_analyze.py

diff --git a/calibcant/vibration_analyze.py b/calibcant/vibration_analyze.py
index 1b8fec256e9ec6d2cc295fb341482e491e68e5df..a1125b6a5c15746a5efab3771b862850a2337e49 100644 (file)
--- a/calibcant/vibration_analyze.py
+++ b/calibcant/vibration_analyze.py
@@ -454,6 +454,11 @@ def breit_wigner_area(A, B, C):
     #  <V(t)**2> = (pi*C) / (2*B*A**2)
     return (_numpy.pi * C) / (2 * B * A**2)
 
+def breit_wigner_resonant_frequency(A, B):
+    if (B**2 >= 2*A**2):
+        return 0  # over- or critically-damped
+    return _numpy.sqrt(A**2 - B**2/2)
+
 def save(filename, group='/', raw=None, config=None,
          deflection_channel_config=None, processed=None):
     specs = [
@@ -479,7 +484,7 @@ def load(filename=None, group='/'):
     return _load(filename=filename, group=group, specs=specs)
 
 def plot(deflection=None, freq_axis=None, power=None, A=None, B=None,
-             C=None, D=0, config=None, analyze=False):
+         C=None, D=0, config=None, analyze=False):
     """Plot 3 subfigures displaying vibration data and analysis.
 
      Time series (Vphoto vs sample index) (show any major drift events),
@@ -520,10 +525,9 @@ def plot(deflection=None, freq_axis=None, power=None, A=None, B=None,
         std = deflection.std()
         pi = _numpy.pi
         exp = _numpy.exp
-        gauss = _numpy.sqrt(2*pi)/std * exp(-0.5*((bins-mean)/std)**2)
+        gauss = exp(-0.5*((bins-mean)/std)**2) / (_numpy.sqrt(2*pi) * std)
         # Matplotlib's normed histogram uses bin heights of n/(len(x)*dbin)
-        dbin = bins[1]-bins[0]
-        hist_axes.plot(bins, gauss/dbin, 'r-')
+        hist_axes.plot(bins, gauss, 'r-')
         hist_axes.autoscale(tight=True)
     if power is not None:
         freq_axes.hold(True)
@@ -544,9 +548,8 @@ def plot(deflection=None, freq_axis=None, power=None, A=None, B=None,
             freq_axes.plot(freq_axis, fitdata, 'b-')
             noisefloor = D + 0*freq_axis;
             freq_axes.plot(freq_axis, noisefloor)
-
-            if B**2 < 2*A**2:
-                res_freq = _numpy.sqrt(A**2 - B**2/2)
+            res_freq = breit_wigner_resonant_frequency(A=A, B=B)
+            if res_freq > 0:
                 freq_axes.axvline(res_freq, color='b', zorder=-1)
 
         freq_axes.set_title('power spectral density %s' % timestamp)