20Generate strain time series of a continuous-wave signal in the detector frame,
21given a function which computes signal phase and amplitudes as functions of time.
23The given function should compute quantities \f$\Delta\phi(t)\f$, \f$a_+(t)\f$,
24and \f$a_\times(t)\f$ such that the plus and cross gravitational-wave
25polarisations \f$h_+(t)\f$ and \f$h_\times(t)\f$ are given by:
27h_+(t) & = & a_+(t) \cos[\phi_0 + \Delta\phi(t)] \; , \\
28h_\times(t) & = & a_\times(t)\sin[\phi_0 + \Delta\phi(t)] \; .
32CWSimulator() can also directly write frame files
and SFT files. Example usage:
36from lalpulsar
import simulateCW
38def waveform(h0, cosi, freq, f1dot):
40 dphi = lal.TWOPI * (freq * dt + f1dot * 0.5 * dt**2)
41 ap = h0 * (1.0 + cosi**2) / 2.0
60wf = waveform(h0, cosi, freq, f1dot)
64for file, i, N
in S.write_sft_files(fmax=32, Tsft=1800, comment=
"simCW"):
65 print(
'Generated SFT file %s (%i of %i)' % (file, i+1, N))
68for file, i, N
in S.write_frame_files(fs=1, Tframe=1800, comment=
"simCW"):
69 print(
'Generated frame file %s (%i of %i)' % (file, i+1, N))
75from __future__ import division, print_function
85from . import git_version
87__author__ = "Karl Wette <karl.wette@ligo.org>"
88__version__ = git_version.id
89__date__ = git_version.date
92class CWSimulator(object):
105 earth_ephem_file="earth00-40-DE405.dat.gz",
106 sun_ephem_file="sun00-40-DE405.dat.gz",
111 Initialise a continuous-wave signal simulator.
113 @param tref: reference time of signal phase at Solar System barycentre,
in GPS seconds
114 (but see
@b tref_at_det)
115 @param tstart: start time of signal,
in GPS seconds
116 @param Tdata: total duration of signal,
in seconds
117 @param waveform: function which computes signal phase
and amplitudes
as functions of time:
118 @b dphi,
@b aplus,
@b across =
@b waveform(dt), where:
119 @b dt = time since reference time
@b tref;
120 @b dphi = phase of signal at time
@b dt relative to reference time
@b tref,
in radians;
121 @b aplus = strain amplitude of plus polarisation at time
@b dt;
122 @b across = strain amplitude of cross polarisation at time
@b dt
123 @param dt_wf: sampling time of the function
@c waveform; this need only be small enough to ensure
124 that
@c dphi,
@c aplus,
and @c across are smoothly interpolated,
and does
not need to make
125 the desired sampling frequency of the output strain time series
126 @param phi0: initial phase of the gravitational-wave signal at
@b tstart,
in radians
127 @param psi: polarisation angle of the gravitational-wave source,
in radians
128 @param alpha: right ascension of the gravitational-wave source,
in radians
129 @param delta: declination of the gravitational-wave source,
in radians
130 @param det_name: name of the gravitational-wave detector to simulate a response
for;
131 e.g.
@c "H1" for LIGO Hanford,
@c "L1" for LIGO Livingston,
@c "V1" for Virgo
132 @param earth_ephem_file: name of file to load Earth ephemeris
from
133 @param sun_ephem_file: name of file to load Sun ephemeris
from
134 @param tref_at_det: default
False;
if True, shift reference time
@b tref so that
@b dt = 0
is
135 @b tref
in @e detector frame instead of Solar System barycentre frame, useful
if e.g. one
136 wants to turn on signal only
for @b dt > 0, one can use
@b tref
as the turn-on time
137 @param extra_comment: additional text to add to comment string
in frame/SFT headers
138 (
not filenames), e.g.
for wrapper script commandlines
156 _, self.
__site_index = lalpulsar.FindCWDetector(det_name,
True)
159 raise ValueError(
"Invalid detector name det_name='%s'" % det_name)
163 self.
__ephemerides = lalpulsar.InitBarycenter(earth_ephem_file, sun_ephem_file)
167 bary_state = lalpulsar.EarthState()
168 bary_input = lalpulsar.BarycenterInput()
169 bary_emit = lalpulsar.EmissionTime()
170 bary_input.tgps = tref
171 bary_input.site = self.
__site
172 for i
in range(0, 3):
173 bary_input.site.location[i] /= lal.C_SI
174 bary_input.alpha = alpha
175 bary_input.delta = delta
176 bary_input.dInv = 0.0
177 lalpulsar.BarycenterEarth(bary_state, bary_input.tgps, self.
__ephemerides)
178 lalpulsar.Barycenter(bary_emit, bary_input, bary_state)
181 tref += bary_emit.deltaT
187 Tpad_wf = 2.0 * lal.AU_SI / lal.C_SI
188 tstart_wf = tstart - Tpad_wf
189 Nwf =
int(math.ceil(
float(Tdata + 2 * Tpad_wf) /
float(dt_wf)))
192 self.
__phi = lal.CreateREAL8TimeSeries(
193 "phi", tstart_wf, 0, dt_wf, lal.DimensionlessUnit, Nwf
198 self.
__a = lal.REAL4TimeVectorSeries()
199 self.
__a.name =
"a+,ax"
200 self.
__a.epoch = tstart_wf
201 self.
__a.deltaT = dt_wf
203 self.
__a.sampleUnits = lal.StrainUnit
204 self.
__a.data = lal.CreateREAL4VectorSequence(Nwf, 2)
207 dt =
float(tstart_wf - tref)
208 for i
in range(0, Nwf):
209 dphi, aplus, across = waveform(dt)
210 self.
__phi.data.data[i] = phi0 + dphi
211 self.
__a.data.data[i][0] = aplus
212 self.
__a.data.data[i][1] = across
216 self.
__wf = lalpulsar.PulsarCoherentGW()
217 self.
__wf.position.system = lal.COORDINATESYSTEM_EQUATORIAL
218 self.
__wf.position.longitude = alpha
219 self.
__wf.position.latitude = delta
225 self.
__detector = lalpulsar.PulsarDetectorResponse()
229 def _simulate_coherent_gw(self, h, noise_sqrt_Sh, noise_seed):
234 if noise_sqrt_Sh > 0:
235 assert noise_seed
is not None
236 assert noise_seed > 0
237 noise_sigma = math.sqrt(0.5 / h.deltaT) * noise_sqrt_Sh
238 lalpulsar.AddGaussianNoise(h, noise_sigma, noise_seed)
240 def get_strain(self, fs, tmin=0, tmax=None, noise_sqrt_Sh=0, noise_seed=None):
242 Generate strain time series of a continuous-wave signal in the detector frame.
244 @param fs: sampling frequency of strain time series,
in Hz
245 @param tmin: start time
for strain time series,
as offsets
from self.
__tstart
246 @param tmax: start time
for strain time series,
as offsets
from self.
__tstart
247 @param noise_sqrt_Sh:
if >0, add Gaussian noise
with square-root single-sided power
248 spectral density given by this value,
in Hz^(-1/2)
249 @param noise_seed: use this seed
for the random number generator used to create noise;
250 required
if noise_sqrt_Sh >0
252 @return (
@b t,
@b h), where:
253 @b t = start of time strain time series,
in GPS seconds;
254 @b h = strain time series
258 if (tmin < 0)
or (tmin >= self.
__Tdata):
259 raise ValueError(
"tmin must be within [0,{}).".
format(self.
__Tdata))
262 elif (tmax <= 0)
or (tmax > self.
__Tdata):
263 raise ValueError(
"tmax must be within (0,{}].".
format(self.
__Tdata))
268 h = lal.CreateREAL4TimeSeries(
269 "h", self.
__tstart + tmin, 0, 1.0 / fs, lal.DimensionlessUnit, Nh
274 return (h.epoch, h.data.data)
277 self, fs, Tblock, noise_sqrt_Sh=0, noise_seed=None, prog_bar=None
280 Generate strain time series of a continuous-wave signal in the detector frame,
in contiguous blocks.
282 @param fs: sampling frequency of strain time series,
in Hz
283 @param Tblock: length of each block,
in seconds; should divide evenly into
@b Tdata
284 @param noise_sqrt_Sh:
if >0, add Gaussian noise
with square-root single-sided power
285 spectral density given by this value,
in Hz^(-1/2)
286 @param noise_seed: use this seed
for the random number generator used to create noise;
287 if None, generate a random seed using os.urandom()
288 @param prog_bar: use to display a progress bar, e.g. prog_bar=tqdm
290 @return (
@b t,
@b h,
@b i,
@b N), where:
291 @b t = start of time strain time series,
in GPS seconds;
292 @b h = strain time series;
293 @b i = block index, starting
from zero;
294 @b N = number of blocks
296 This
is a Python generator function
and so should be called
as follows:
299 for t, h, i, N
in S.get_strain_blocks(...):
306 if Tblock * Nblock > self.
__Tdata:
308 "Length of block Tblock=%g does not divide evenly into Tdata=%g"
315 while not (0 < noise_seed
and noise_seed < lal.LAL_INT4_MAX - Nblock):
316 noise_seed = int.from_bytes(os.urandom(4), sys.byteorder, signed=
False)
319 blocks =
range(0, Nblock)
320 if prog_bar
is not None:
327 for iblock
in blocks:
332 noise_sqrt_Sh=noise_sqrt_Sh,
333 noise_seed=noise_seed + iblock,
335 yield epoch, hoft, iblock, Nblock
349 Write frame files [1] containing strain time series of a continuous-wave signal.
351 The strain time series is written
as double-precision post-processed data (ProcData) channel named
352 <tt><detector>:SIMCW-STRAIN</tt>,
353 where <tt><detector></tt>
is the 2-character detector prefix (e.g. <tt>H1</tt>
for LIGO Hanford,
354 <tt>L1</tt>
for LIGO Livingston, <tt>V1</tt>
for Virgo).
356 @param fs: sampling frequency of strain time series,
in Hz
357 @param Tframe: length of each frame,
in seconds; should divide evenly into
@b Tdata
358 @param comment: frame file name comment, may only contain A-Z, a-z, 0-9, _, +,
359 @param out_dir: output directory to write frame files into
360 @param noise_sqrt_Sh:
if >0, add Gaussian noise
with square-root single-sided power
361 spectral density given by this value,
in Hz^(-1/2)
362 @param noise_seed: use this seed
for the random number generator used to create noise;
363 if None, generate a random seed using os.urandom()
364 @param prog_bar: use to display a progress bar, e.g. prog_bar=tqdm
366 @return (
@b file,
@b i,
@b N), where:
367 @b file = name of frame file just written;
368 @b i = frame file index, starting
from zero;
369 @b N = number of frame files
371 This
is a Python generator function
and so should be called
as follows:
374 for file, i, N
in S.write_frame_files(...):
378 [1] https://dcc.ligo.org/LIGO-T970130/public
384 raise ImportError(
"SWIG wrappings of LALFrame cannot be imported")
387 valid_comment = re.compile(
r"^[A-Za-z0-9_+#]+$")
388 if not valid_comment.match(comment):
390 "Frame file comment='%s' may only contain A-Z, a-z, 0-9, _, +, # characters"
399 noise_sqrt_Sh=noise_sqrt_Sh,
400 noise_seed=noise_seed,
405 frame_h_channel =
"%s:SIMCW-STRAIN" % self.
__site.frDetector.prefix
406 frame_h = lal.CreateREAL8TimeSeries(
407 frame_h_channel, t, 0, 1.0 / fs, lal.DimensionlessUnit, len(h)
410 frame_h.data.data = h
413 frame_src = self.
__site.frDetector.prefix[0]
415 frame_t0 =
int(t.gpsSeconds)
416 frame_Tdata =
int(math.ceil(
float(t + Tframe)) - math.floor(
float(t)))
417 frame_name =
"%s-%s-%u-%u.gwf" % (
423 frame_path = os.path.join(out_dir, frame_name)
427 frame = lalframe.FrameNew(
428 t, Tframe, self.__class__.__name__, -1, i, frame_det_bits
432 lalframe.FrameAddREAL8TimeSeriesProcData(frame, frame_h)
435 lalframe.FrameAddFrHistory(frame,
"origin", self.
__origin_str)
438 lalframe.FrameWrite(frame, frame_path)
441 yield frame_path, i, N
449 window="rectangular",
455 Generate SFTs [2] containing strain time series of a continuous-wave signal.
457 @param fmax: maximum SFT frequency,
in Hz
458 @param Tsft: length of each SFT,
in seconds; should divide evenly into
@b Tdata
459 @param noise_sqrt_Sh:
if >0, add Gaussian noise
with square-root single-sided power
460 spectral density given by this value,
in Hz^(-1/2)
461 @param noise_seed: use this seed
for the random number generator used to create noise;
462 if None, generate a random seed using os.urandom()
463 @param window:
if not None, window the time series before performing the FFT, using
465 @param window_param: parameter
for the window function given by
@b window,
if needed
466 @param prog_bar: use to display a progress bar, e.g. prog_bar=tqdm
467 @param fmin:
if >0, minimum SFT frequency,
in Hz
468 (note that this only modifies the output SFT bandwidth; internally a full-band SFT
is still produced by sampling at 2*fmax)
470 @return (
@b sft,
@b i,
@b N), where:
472 @b i = SFT file index, starting
from zero;
473 @b N = number of SFTs
475 This
is a Python generator function
and so should be called
as follows:
478 for sft, i, N
in S.get_sfts(...):
482 [2] https://dcc.ligo.org/LIGO-T040164/public
486 sft_ts = lalpulsar.CreateTimestampVector(1)
495 noise_sqrt_Sh=noise_sqrt_Sh,
496 noise_seed=noise_seed,
501 sft_name = self.
__site.frDetector.prefix
502 sft_h = lal.CreateREAL8TimeSeries(
503 sft_name, t, 0, 1.0 / sft_fs, lal.DimensionlessUnit, len(h)
510 sft_vect = lalpulsar.MakeSFTsFromREAL8TimeSeries(
511 sft_h, sft_ts, window, window_param
515 sft_vect = lalpulsar.ExtractStrictBandFromSFTVector(
516 sft_vect, fMin=fmin, Band=fmax - fmin
520 yield sft_vect.data[0], i, N
530 window="rectangular",
536 Write SFT files [2] containing strain time series of a continuous-wave signal.
538 @param fmax: maximum SFT frequency,
in Hz
539 @param Tsft: length of each SFT,
in seconds; should divide evenly into
@b Tdata
540 @param comment: SFT file name comment, may only contain A-Z, a-z, 0-9 characters
541 @param out_dir: output directory to write SFT files into
542 @param noise_sqrt_Sh:
if >0, add Gaussian noise
with square-root single-sided power
543 spectral density given by this value,
in Hz^(-1/2)
544 @param noise_seed: use this seed
for the random number generator used to create noise;
545 if None, generate a random seed using os.urandom()
546 @param window:
if not None, window the time series before performing the FFT, using
548 @param window_param: parameter
for the window function given by
@b window,
if needed
549 @param prog_bar: use to display a progress bar, e.g. prog_bar=tqdm
550 @param fmin:
if >0, minimum SFT frequency,
in Hz
551 (note that this only modifies the output SFT bandwidth; internally a full-band SFT
is still produced by sampling at 2*fmax)
553 @return (
@b file,
@b i,
@b N), where:
554 @b file = name of SFT file just written;
555 @b i = SFT file index, starting
from zero;
556 @b N = number of SFT files
558 This
is a Python generator function
and so should be called
as follows:
561 for file, i, N
in S.write_sft_files(...):
565 [2] https://dcc.ligo.org/LIGO-T040164/public
569 valid_comment = re.compile(
r"^[A-Za-z0-9]+$")
570 if not valid_comment.match(comment):
572 "SFT file comment='%s' may only contain A-Z, a-z, 0-9 characters"
577 spec = lalpulsar.SFTFilenameSpec()
579 spec.window_type = window
580 spec.window_param = window_param
581 spec.privMisc = comment
587 noise_sqrt_Sh=noise_sqrt_Sh,
588 noise_seed=noise_seed,
590 window_param=window_param,
595 lalpulsar.WriteSFT2StandardFile(sft, spec, self.
__origin_str)
598 sft_path = lalpulsar.BuildSFTFilenameFromSpec(spec)
def write_sft_files(self, fmax, Tsft, comment="simCW", out_dir=".", noise_sqrt_Sh=0, noise_seed=None, window="rectangular", window_param=0, prog_bar=None, fmin=0)
Write SFT files [2] containing strain time series of a continuous-wave signal.
def get_strain(self, fs, tmin=0, tmax=None, noise_sqrt_Sh=0, noise_seed=None)
Generate strain time series of a continuous-wave signal in the detector frame.
def get_sfts(self, fmax, Tsft, noise_sqrt_Sh=0, noise_seed=None, window="rectangular", window_param=0, prog_bar=None, fmin=0)
Generate SFTs [2] containing strain time series of a continuous-wave signal.
def _simulate_coherent_gw(self, h, noise_sqrt_Sh, noise_seed)
def write_frame_files(self, fs, Tframe, comment="simCW", out_dir=".", noise_sqrt_Sh=0, noise_seed=None, prog_bar=None)
Write frame files [1] containing strain time series of a continuous-wave signal.
def __init__(self, tref, tstart, Tdata, waveform, dt_wf, phi0, psi, alpha, delta, det_name, earth_ephem_file="earth00-40-DE405.dat.gz", sun_ephem_file="sun00-40-DE405.dat.gz", tref_at_det=False, extra_comment=None)
Initialise a continuous-wave signal simulator.
def get_strain_blocks(self, fs, Tblock, noise_sqrt_Sh=0, noise_seed=None, prog_bar=None)
Generate strain time series of a continuous-wave signal in the detector frame, in contiguous blocks.
REAL8Window * XLALCreateNamedREAL8Window(const char *windowName, REAL8 beta, UINT4 length)