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