Note
Click here to download the full example code
Directed grid search: Monochromatic sourceΒΆ
Search for a monochromatic (no spindown) signal using a parameter space grid (i.e. no MCMC).
8 import pyfstat
9 import numpy as np
10 import matplotlib.pyplot as plt
11 import os
12
13 label = "PyFstat_example_grid_search_F0"
14 outdir = os.path.join("PyFstat_example_data", label)
15
16 F0 = 30.0
17 F1 = 0
18 F2 = 0
19 Alpha = 1.0
20 Delta = 1.5
21
22 # Properties of the GW data
23 depth = 70
24 sqrtS = "1e-23"
25 h0 = float(sqrtS) / depth
26 cosi = 0
27 IFOs = "H1"
28 # IFOs = "H1,L1"
29 sqrtSX = ",".join(np.repeat(sqrtS, len(IFOs.split(","))))
30 tstart = 1000000000
31 duration = 100 * 86400
32 tend = tstart + duration
33 tref = 0.5 * (tstart + tend)
34
35 data = pyfstat.Writer(
36 label=label,
37 outdir=outdir,
38 tref=tref,
39 tstart=tstart,
40 F0=F0,
41 F1=F1,
42 F2=F2,
43 duration=duration,
44 Alpha=Alpha,
45 Delta=Delta,
46 h0=h0,
47 cosi=cosi,
48 sqrtSX=sqrtSX,
49 detectors=IFOs,
50 )
51 data.make_data()
52
53 m = 0.001
54 dF0 = np.sqrt(12 * m) / (np.pi * duration)
55 DeltaF0 = 800 * dF0
56 F0s = [F0 - DeltaF0 / 2.0, F0 + DeltaF0 / 2.0, dF0]
57 F1s = [F1]
58 F2s = [F2]
59 Alphas = [Alpha]
60 Deltas = [Delta]
61 search = pyfstat.GridSearch(
62 label,
63 outdir,
64 os.path.join(outdir, "*" + label + "*sft"),
65 F0s,
66 F1s,
67 F2s,
68 Alphas,
69 Deltas,
70 tref,
71 tstart,
72 tend,
73 )
74 search.run()
75
76 print("Plotting 2F(F0)...")
77 fig, ax = plt.subplots()
78 frequencies = search.data["F0"]
79 twoF = search.data["twoF"]
80 # mismatch = np.sign(x-F0)*(duration * np.pi * (x - F0))**2 / 12.0
81 ax.plot(frequencies, twoF, "k", lw=1)
82 DeltaF = frequencies - F0
83 sinc = np.sin(np.pi * DeltaF * duration) / (np.pi * DeltaF * duration)
84 A = np.abs((np.max(twoF) - 4) * sinc ** 2 + 4)
85 ax.plot(frequencies, A, "-r", lw=1)
86 ax.set_ylabel("$\\widetilde{2\\mathcal{F}}$")
87 ax.set_xlabel("Frequency")
88 ax.set_xlim(F0s[0], F0s[1])
89 dF0 = np.sqrt(12 * 1) / (np.pi * duration)
90 xticks = [F0 - 10 * dF0, F0, F0 + 10 * dF0]
91 ax.set_xticks(xticks)
92 xticklabels = ["$f_0 {-} 10\\Delta f$", "$f_0$", "$f_0 {+} 10\\Delta f$"]
93 ax.set_xticklabels(xticklabels)
94 plt.tight_layout()
95 fig.savefig(os.path.join(outdir, label + "_1D.png"), dpi=300)
Total running time of the script: ( 0 minutes 0.000 seconds)