rewords the author's note

comprehensive.tex

@@ -488,7 +488,7 @@ The observation of gravity waves was a great moment in scientific history, and t
 \vspace{5pt}
 \hrule
 \vspace{5pt}
-\emph{Author's note}: Shortly after the writing of this paper, the LIGO Scientific Collaboration announced a third detection of a binary black hole system designated as GW170104\cite{GW170104}. The event was reconstructed with a total mass of 50$M_{\astrosun}$. The new data updates the coelescence rate to \err{103}{110}{63}\si{\per\giga\parsec\cubed\per\year}, with the $m_1$ mass distribution parameter, $\alpha$, updated to \err{2.3}{1.3}{1.4}.
+\emph{Author's note}: Shortly after the writing of this paper, the LIGO Scientific Collaboration announced a third detection of a binary black hole system designated as GW170104\cite{GW170104}. The event was reconstructed with a total mass of 50$M_{\astrosun}$. The new data updates the coelescence rate to \err{103}{110}{63}\si{\per\giga\parsec\cubed\per\year}, and the $m_1$ mass distribution parameter, $\alpha$, to \err{2.3}{1.3}{1.4}.
 \vspace{5pt}
 \hrule
 \vspace{-10pt}

scripts/pdh_demo2.py

@@ -5,22 +5,22 @@ from matplotlib import animation
 import numpy as np
 from numpy import cos, pi
 
-matplotlib.rcParams.update({'font.size': 15})
+# matplotlib.rcParams.update({'font.size': 15})
 
 c = 3E8  # m/s
 L = 1.0  # m
-FSR = 2*pi*c/(2*L)  # Hz
+FSR = 2*pi*c/(2*L)  # (rad/s)
 
 A0 = 1.0
-w0 = 1E9*FSR
+w0 = 1E1*FSR
 f0 = w0 / 2 / pi
-a_mod = 200
-w_mod = 0.0002 * w0
+a_mod = 2
+w_mod = w0 / 3000
 f_mod = w_mod / 2 / pi
 
 
 def R_e(w):
-    T = 0.05
+    T = 0.25
     R = 1-T
     return 1 - T**2 / (1 + R**2 - 2*R*cos(w*2*L/c))
 
@@ -29,9 +29,84 @@ def get_ref_power(freqs, powers):
     return sum(power*R_e(freq) for freq, power in zip(freqs, powers))
 
 
+def main3():
+    N = 2000000
+    ts = np.linspace(0, 10000*(1/f0), N)
+    mod = a_mod*np.sin(w_mod*ts)
+    Es = A0*np.cos((w0 + (-0.012)*FSR)*ts + mod)
+    maxima = ts[np.r_[True, Es[1:] < Es[:-1]] &
+                np.r_[Es[:-1] < Es[1:], True]]
+    maxima = maxima[1:-1]
+
+    n_periods = 100
+    periods = (maxima[n_periods:] - maxima[:-n_periods])/n_periods
+    centers = maxima[:-n_periods] + (periods/2)
+    powers = R_e(2*pi/periods)
+    freqs = (2*pi/periods - w0)/FSR
+
+
+    def moving_average(a, n=3):
+        ret = np.cumsum(a, dtype=float)
+        ret[n:] = ret[n:] - ret[:-n]
+        return ret[n - 1:] / n
+
+    smooth = 25
+    freqs = moving_average(freqs, n=smooth)
+    powers = moving_average(powers, n=smooth)
+    centers = centers[smooth-1:]
+
+    Rs_ws = np.linspace(w0-0.5*FSR, w0+0.5*FSR, 300)
+    Rs_Rs = R_e(Rs_ws)
+    Rs_ws = (Rs_ws - w0)/FSR
+
+    fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, ncols=1)
+
+    n_frames = 300
+    p_max = 0.16
+
+    def do_plot(i):
+        print(i)
+        idx = int(i/n_frames * len(centers))
+        t = centers[idx]
+
+        f = freqs[idx]
+        p = powers[idx]
+
+        ax1.clear()
+        ax1.plot(Rs_ws, Rs_Rs, 'b')
+        ax1.plot([f], [p], 'r.')
+        ax1.set_xlim((-0.1, 0.1))
+        ax1.set_xlabel(r"$\Delta\omega(\mathrm{FSR})$")
+        ax1.set_ylim((0, p_max))
+
+        ax2.clear()
+        ax2.plot(centers*f_mod, powers, 'b')
+        ax2.set_ylabel('Refl. Power')
+        ax2.set_ylim((0, p_max))
+        ax2.set_xlim((0, ts[-1]*f_mod))
+        ax2.set_xticks([])
+        rect = plt.Rectangle((0, 0), t*f_mod, 1, facecolor='#aaaaaa')
+        ax2.add_patch(rect)
+
+        ax3.clear()
+        ax3.plot(ts*f_mod, mod, 'r')
+        ax3.set_ylabel('Modulation')
+        ax3.set_ylim((-a_mod, a_mod))
+        ax3.set_xlim((0, ts[-1]*f_mod))
+        ax3.set_xlabel(r"$t(\tau_{\mathrm{mod}})$")
+        rect = plt.Rectangle((0, -a_mod),
+                             t*f_mod, 2*a_mod, facecolor='#aaaaaa')
+        ax3.add_patch(rect)
+
+    # do_plot(0)
+    # plt.show()
+    anim = animation.FuncAnimation(fig, do_plot, n_frames)
+    anim.save('output.mp4', fps=30, dpi=270)
+
+
 def main():
-    sample_n = 2000000
-    ts = np.linspace(0, 10000*(1/f0), sample_n)
+    N = 2000000
+    ts = np.linspace(0, 10000*(1/f0), N)
     dt = ts[1]-ts[0]
     mod = a_mod*np.sin(w_mod*ts)
     Es = A0*np.cos(w0*ts + mod)
@@ -40,42 +115,55 @@ def main():
 
     n_groups = 100
     win_size = int((1/f_mod) / dt / 12)
-    win_spacing = (sample_n - win_size) // (n_groups+2)
+    win_spacing = (N - win_size) // (n_groups-1)
     fig, (ax0, ax1, ax2) = plt.subplots(nrows=3, ncols=1)
+    fig.set_tight_layout(True)
+
+    Rs_ws = np.linspace(w0-0.5*FSR, w0+0.5*FSR, 300)
+    Rs_Rs = R_e(Rs_ws)
+    Rs_ws = (Rs_ws - w0)/FSR
 
     def do_plot(i):
         print(i)
-        # ax0.clear()
-        ax1.clear()
+
         n = ts.size // n_groups
-        win_start = (win_size//2) + i*win_spacing
+        win_start = i*win_spacing
         win_end = win_start + win_size
-        print(f"{win_start},{win_end}")
+        # print(f"{win_start},{win_end}")
 
         # sec_ts = ts[win_start:win_end]
         sec = Es[win_start:win_end]
         # mod_sec = mod[win_start:win_end]
         sec_ft = np.fft.rfftn(sec)
-        freq = np.fft.fftfreq(sec.size, d=dt)[:n//2]
-        power = (sec_ft[:n//2].real**2 + sec_ft[:n//2].imag**2) / win_size
+        freq = 2*pi*np.fft.fftfreq(sec.size, d=dt)[:n//2]
+        power = (sec_ft[:n//2]*sec_ft[:n//2].conj()).real / win_size
+        # power = (sec_ft[:n//2].real**2 + sec_ft[:n//2].imag**2) / win_size
 
-        # ax0.plot(sec_ts, mod_sec)
-        ax0.plot(ts[(win_start+win_end) // 2], mod[(win_start+win_end) // 2], 'k.')
+        ax0.plot(ts[(win_start+win_end) // 2],
+                 mod[(win_start+win_end) // 2], 'k.')
         ax0.set_ylim(-1.1*a_mod, 1.1*a_mod)
         ax0.set_xlim(0, ts[-1])
 
-        ax1.plot(2*pi*freq/w0, power)
+        ax1.clear()
+        ax1.plot(freq/w0, power)
         ax1.set_xlim(0.80, 1.20)
 
-        ref_power = get_ref_power(freq, power)
-        ax2.plot(ts[(win_start+win_end) // 2], ref_power, 'r.')
-        ax2.set_xlim(0, ts[-1])
+        mean = np.average(freq, weights=power)
+        print('{:2.2f} : {:2.2E}'.format((mean-w0)/FSR, R_e(mean*w0)))
 
-        plt.tight_layout()
+        ax2.clear()
+        ax2.plot(Rs_ws, Rs_Rs, 'b')
+        ax2.plot((mean-w0)/FSR, R_e(mean), 'r.')
+        ax2.set_xlim(-.5, .5)
+        # ref_power = get_ref_power(freq, power)
+        # ax2.plot(ts[(win_start+win_end) // 2], ref_power, 'r.')
+        # ax2.set_xlim(0, ts[-1])
 
     anim = animation.FuncAnimation(fig, do_plot, n_groups)
     anim.save("output.mp4")
 
 
 if __name__ == '__main__':
-    main()
+    # main()
+    # main2()
+    main3()