Directed grid search: Monochromatic sourceΒΆ

Search for a monochromatic (no spindown) signal using a parameter space grid (i.e. no MCMC).

 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
 import pyfstat
 import numpy as np
 import matplotlib.pyplot as plt
 import os

 label = "PyFstat_example_grid_search_F0"
 outdir = os.path.join("PyFstat_example_data", label)

 F0 = 30.0
 F1 = 0
 F2 = 0
 Alpha = 1.0
 Delta = 1.5

 # Properties of the GW data
 depth = 70
 sqrtS = "1e-23"
 h0 = float(sqrtS) / depth
 cosi = 0
 IFOs = "H1"
 # IFOs = "H1,L1"
 sqrtSX = ",".join(np.repeat(sqrtS, len(IFOs.split(","))))
 tstart = 1000000000
 duration = 100 * 86400
 tend = tstart + duration
 tref = 0.5 * (tstart + tend)

 data = pyfstat.Writer(
     label=label,
     outdir=outdir,
     tref=tref,
     tstart=tstart,
     F0=F0,
     F1=F1,
     F2=F2,
     duration=duration,
     Alpha=Alpha,
     Delta=Delta,
     h0=h0,
     cosi=cosi,
     sqrtSX=sqrtSX,
     detectors=IFOs,
 )
 data.make_data()

 m = 0.001
 dF0 = np.sqrt(12 * m) / (np.pi * duration)
 DeltaF0 = 800 * dF0
 F0s = [F0 - DeltaF0 / 2.0, F0 + DeltaF0 / 2.0, dF0]
 F1s = [F1]
 F2s = [F2]
 Alphas = [Alpha]
 Deltas = [Delta]
 search = pyfstat.GridSearch(
     label,
     outdir,
     os.path.join(outdir, "*" + label + "*sft"),
     F0s,
     F1s,
     F2s,
     Alphas,
     Deltas,
     tref,
     tstart,
     tend,
 )
 search.run()

 print("Plotting 2F(F0)...")
 fig, ax = plt.subplots()
 frequencies = search.data["F0"]
 twoF = search.data["twoF"]
 # mismatch = np.sign(x-F0)*(duration * np.pi * (x - F0))**2 / 12.0
 ax.plot(frequencies, twoF, "k", lw=1)
 DeltaF = frequencies - F0
 sinc = np.sin(np.pi * DeltaF * duration) / (np.pi * DeltaF * duration)
 A = np.abs((np.max(twoF) - 4) * sinc ** 2 + 4)
 ax.plot(frequencies, A, "-r", lw=1)
 ax.set_ylabel("$\\widetilde{2\\mathcal{F}}$")
 ax.set_xlabel("Frequency")
 ax.set_xlim(F0s[0], F0s[1])
 dF0 = np.sqrt(12 * 1) / (np.pi * duration)
 xticks = [F0 - 10 * dF0, F0, F0 + 10 * dF0]
 ax.set_xticks(xticks)
 xticklabels = ["$f_0 {-} 10\\Delta f$", "$f_0$", "$f_0 {+} 10\\Delta f$"]
 ax.set_xticklabels(xticklabels)
 plt.tight_layout()
 fig.savefig(os.path.join(outdir, label + "_1D.png"), dpi=300)

Total running time of the script: ( 0 minutes 0.000 seconds)

Gallery generated by Sphinx-Gallery