22 from numpy.testing
import assert_allclose
28 from lalpulsar
import simulateCW
as simCW
31 earth_ephem_file = os.path.join(
32 os.environ[
"LAL_TEST_PKGDATADIR"],
"earth00-40-DE405.dat.gz"
34 sun_ephem_file = os.path.join(
35 os.environ[
"LAL_TEST_PKGDATADIR"],
"sun00-40-DE405.dat.gz"
37 ephemerides = lalpulsar.InitBarycenter(earth_ephem_file, sun_ephem_file)
69 fsmin = freq - (Nfs - 1) * 0.5 * dfreq
75 dphi = lal.TWOPI * (freq * dt + f1dot * 0.5 * dt**2)
76 ap = h0 * (1.0 + cosi**2) / 2.0
86 sft_ts = lalpulsar.MakeTimestamps(tstart, Tdata, Tsft, 0)
87 sft_catalog = lalpulsar.AddToFakeSFTCatalog(
None, detector, sft_ts)
90 Fstat_opt_args = lalpulsar.FstatOptionalArgs(lalpulsar.FstatOptionalArgsDefaults)
91 Fstat_opt_args.FstatMethod = lalpulsar.FMETHOD_DEMOD_BEST
94 Fstat_signal = lalpulsar.CreatePulsarParamsVector(1)
95 Fstat_signal.data[0].Amp.aPlus = 0.5 * h0 * (1.0 + cosi * cosi)
96 Fstat_signal.data[0].Amp.aCross = h0 * cosi
97 Fstat_signal.data[0].Amp.psi = psi
98 Fstat_signal.data[0].Amp.phi0 = phi0
99 Fstat_signal.data[0].Doppler.refTime = tref
100 Fstat_signal.data[0].Doppler.Alpha = alpha
101 Fstat_signal.data[0].Doppler.Delta = delta
102 Fstat_signal.data[0].Doppler.fkdot[0] = freq
103 Fstat_signal.data[0].Doppler.fkdot[1] = f1dot
104 Fstat_opt_args.injectSources = Fstat_signal
105 Fstat_assume_noise = lalpulsar.MultiNoiseFloor()
106 Fstat_assume_noise.length = 1
107 Fstat_assume_noise.sqrtSn[0] = assume_sqrtSh
108 Fstat_opt_args.assumeSqrtSX = Fstat_assume_noise
111 Fstat_input = lalpulsar.CreateFstatInput(
112 sft_catalog, fcmin, fcmax, dfreq, ephemerides, Fstat_opt_args
115 doppler = lalpulsar.PulsarDopplerParams(Fstat_signal.data[0].Doppler)
116 doppler.fkdot[0] = fsmin
119 Fstat_res = lalpulsar.ComputeFstat(
120 Fstat_res, Fstat_input, doppler, Nfs, lalpulsar.FSTATQ_2F
123 return Fstat_res.twoF
132 S = simCW.CWSimulator(
143 earth_ephem_file=earth_ephem_file,
144 sun_ephem_file=sun_ephem_file,
149 for path, i, N
in S.write_sft_files(fmax, Tsft,
"simCWTestFstatSpindown"):
150 allpaths.append(path)
153 sft_catalog = lalpulsar.SFTdataFind(
154 "V-1_V1_1800SFT_simCWTestFstatSpindown-*.sft",
None
158 Fstat_opt_args = lalpulsar.FstatOptionalArgs(lalpulsar.FstatOptionalArgsDefaults)
159 Fstat_opt_args.FstatMethod = lalpulsar.FMETHOD_DEMOD_BEST
160 Fstat_assume_noise = lalpulsar.MultiNoiseFloor()
161 Fstat_assume_noise.length = 1
162 Fstat_assume_noise.sqrtSn[0] = assume_sqrtSh
163 Fstat_opt_args.assumeSqrtSX = Fstat_assume_noise
166 Fstat_input = lalpulsar.CreateFstatInput(
167 sft_catalog, fcmin, fcmax, dfreq, ephemerides, Fstat_opt_args
170 doppler = lalpulsar.PulsarDopplerParams()
171 doppler.refTime = tref
172 doppler.Alpha = alpha
173 doppler.Delta = delta
174 doppler.fkdot[0] = fsmin
175 doppler.fkdot[1] = f1dot
178 Fstat_res = lalpulsar.ComputeFstat(
179 Fstat_res, Fstat_input, doppler, Nfs, lalpulsar.FSTATQ_2F
183 for path
in allpaths:
186 return Fstat_res.twoF
193 "F-statistic of reference spindown signal:\n %s\n"
194 % (
", ".join([
"%0.4f"] * len(Fstat_ref)) % tuple(Fstat_ref))
200 "F-statistic of spindown signal from simulateCW:\n %s\n"
201 % (
" ".join([
"%0.4f"] * len(Fstat_simCW)) % tuple(Fstat_simCW))
205 Fstat_rmserr = np.sqrt(np.mean((1 - (Fstat_simCW / Fstat_ref)) ** 2))
210 err_msg=
"Root-mean-square error between F-statistics: %0.2e" % Fstat_rmserr,
222 S = simCW.CWSimulator(
226 waveform(h0=0, cosi=0, freq=0, f1dot=0),
233 earth_ephem_file=earth_ephem_file,
234 sun_ephem_file=sun_ephem_file,
240 paths[
"None"] = [[], []]
241 paths[
"1"] = [[], []]
242 for path, i, N
in S.write_sft_files(
243 fmax, Tsft, noise_sqrt_Sh=asd, comment=
"simCWTestSeedNoneRound1"
245 allpaths.append(path)
246 paths[
"None"][0].append(path)
248 for path, i, N
in S.write_sft_files(
249 fmax, Tsft, noise_sqrt_Sh=asd, comment=
"simCWTestSeedNoneRound2"
251 allpaths.append(path)
252 paths[
"None"][1].append(path)
254 for path, i, N
in S.write_sft_files(
255 fmax, Tsft, noise_sqrt_Sh=asd, comment=
"simCWTestSeed1Round1", noise_seed=1
257 allpaths.append(path)
258 paths[
"1"][0].append(path)
260 for path, i, N
in S.write_sft_files(
261 fmax, Tsft, noise_sqrt_Sh=asd, comment=
"simCWTestSeed1Round2", noise_seed=1
263 allpaths.append(path)
264 paths[
"1"][1].append(path)
267 return_codes[
"None"] = np.zeros(2)
268 return_codes[
"1"] = np.zeros(2)
269 for seed
in [
"None",
"1"]:
271 catalog0 = lalpulsar.SFTdataFind(paths[seed][0][i],
None)
272 sft0 = lalpulsar.LoadSFTs(catalog0, -1, -1)
273 catalog1 = lalpulsar.SFTdataFind(paths[seed][1][i],
None)
274 sft1 = lalpulsar.LoadSFTs(catalog1, -1, -1)
276 np.mean(abs(sft0.data[0].data.data - sft1.data[0].data.data) ** 2)
278 return_codes[seed][i] = rms > 0
279 for code
in return_codes[
"None"]:
282 for code
in return_codes[
"1"]:
287 for path
in allpaths:
291 if __name__ ==
"__main__":
292 args = sys.argv[1:]
or [
"-v",
"-rs",
"--junit-xml=junit-simulateCW.xml"]
293 sys.exit(pytest.main(args=[__file__] + args))
def waveform(h0, cosi, freq, f1dot)
def compute_Fstat_spindown_reference()
def compute_Fstat_spindown_simulateCW()