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