Note
Go to the end 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).
9 import os
10
11 import matplotlib.pyplot as plt
12 import numpy as np
13
14 import pyfstat
15
16 label = "PyFstatExampleGridSearchF0"
17 outdir = os.path.join("PyFstat_example_data", label)
18 logger = pyfstat.set_up_logger(label=label, outdir=outdir)
19
20 # Properties of the GW data
21 sqrtS = "1e-23"
22 IFOs = "H1"
23 # IFOs = "H1,L1"
24 sqrtSX = ",".join(np.repeat(sqrtS, len(IFOs.split(","))))
25 tstart = 1000000000
26 duration = 100 * 86400
27 tend = tstart + duration
28 tref = 0.5 * (tstart + tend)
29
30 # parameters for injected signals
31 depth = 70
32 inj = {
33 "tref": tref,
34 "F0": 30.0,
35 "F1": 0,
36 "F2": 0,
37 "Alpha": 1.0,
38 "Delta": 1.5,
39 "h0": float(sqrtS) / depth,
40 "cosi": 0.0,
41 }
42
43 data = pyfstat.Writer(
44 label=label,
45 outdir=outdir,
46 tstart=tstart,
47 duration=duration,
48 sqrtSX=sqrtSX,
49 detectors=IFOs,
50 **inj,
51 )
52 data.make_data()
53
54 m = 0.001
55 dF0 = np.sqrt(12 * m) / (np.pi * duration)
56 DeltaF0 = 800 * dF0
57 F0s = [inj["F0"] - DeltaF0 / 2.0, inj["F0"] + DeltaF0 / 2.0, dF0]
58 F1s = [inj["F1"]]
59 F2s = [inj["F2"]]
60 Alphas = [inj["Alpha"]]
61 Deltas = [inj["Delta"]]
62 search = pyfstat.GridSearch(
63 label=label,
64 outdir=outdir,
65 sftfilepattern=data.sftfilepath,
66 F0s=F0s,
67 F1s=F1s,
68 F2s=F2s,
69 Alphas=Alphas,
70 Deltas=Deltas,
71 tref=tref,
72 minStartTime=tstart,
73 maxStartTime=tend,
74 )
75 search.run()
76
77 # report details of the maximum point
78 max_dict = search.get_max_twoF()
79 logger.info(
80 "max2F={:.4f} from GridSearch, offsets from injection: {:s}.".format(
81 max_dict["twoF"],
82 ", ".join(
83 [
84 "{:.4e} in {:s}".format(max_dict[key] - inj[key], key)
85 for key in max_dict.keys()
86 if not key == "twoF"
87 ]
88 ),
89 )
90 )
91 search.generate_loudest()
92
93 logger.info("Plotting 2F(F0)...")
94 fig, ax = plt.subplots()
95 frequencies = search.data["F0"]
96 twoF = search.data["twoF"]
97 # mismatch = np.sign(x-inj["F0"])*(duration * np.pi * (x - inj["F0"]))**2 / 12.0
98 ax.plot(frequencies, twoF, "k", lw=1)
99 DeltaF = frequencies - inj["F0"]
100 sinc = np.sin(np.pi * DeltaF * duration) / (np.pi * DeltaF * duration)
101 A = np.abs((np.max(twoF) - 4) * sinc**2 + 4)
102 ax.plot(frequencies, A, "-r", lw=1)
103 ax.set_ylabel("$\\widetilde{2\\mathcal{F}}$")
104 ax.set_xlabel("Frequency")
105 ax.set_xlim(F0s[0], F0s[1])
106 dF0 = np.sqrt(12 * 1) / (np.pi * duration)
107 xticks = [inj["F0"] - 10 * dF0, inj["F0"], inj["F0"] + 10 * dF0]
108 ax.set_xticks(xticks)
109 xticklabels = ["$f_0 {-} 10\\Delta f$", "$f_0$", "$f_0 {+} 10\\Delta f$"]
110 ax.set_xticklabels(xticklabels)
111 plt.tight_layout()
112 fig.savefig(os.path.join(outdir, label + "_1D.png"), dpi=300)