create_your_own_source_model.py 2.27 KB
Newer Older
1 2
#!/bin/python
"""
3
A script to demonstrate how to use your own source model
4 5 6 7 8
"""
from __future__ import division, print_function
import tupak
import numpy as np

9
# First set up logging and some output directories and labels
10
tupak.utils.setup_logger()
11 12 13 14
outdir = 'outdir'
label = 'create_your_own_source_model'
sampling_frequency = 4096
time_duration = 1
15 16


17 18
# Here we define out source model - this is the sine-Gaussian model in the
# frequency domain.
19 20 21 22 23 24 25
def sine_gaussian(f, A, f0, tau, phi0, geocent_time, ra, dec, psi):
    arg = -(np.pi * tau * (f-f0))**2 + 1j * phi0
    plus = np.sqrt(np.pi) * A * tau * np.exp(arg) / 2.
    cross = plus * np.exp(1j*np.pi/2)
    return {'plus': plus, 'cross': cross}


26
# We now define some parameters that we will inject and then a waveform generator
27
injection_parameters = dict(A=1e-23, f0=100, tau=1, phi0=0, geocent_time=0,
28
                            ra=0, dec=0, psi=0)
29 30 31 32
waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=time_duration,
                                                                sampling_frequency=sampling_frequency,
                                                                frequency_domain_source_model=sine_gaussian,
                                                                parameters=injection_parameters)
33 34 35 36 37 38 39 40 41 42 43 44
hf_signal = waveform_generator.frequency_domain_strain()

# Set up interferometers.
IFOs = [tupak.detector.get_interferometer_with_fake_noise_and_injection(
    name, injection_polarizations=hf_signal,
    injection_parameters=injection_parameters, time_duration=time_duration,
    sampling_frequency=sampling_frequency, outdir=outdir)
    for name in ['H1', 'L1', 'V1']]

# Here we define the priors for the search. We use the injection parameters
# except for the amplitude, f0, and geocent_time
prior = injection_parameters.copy()
45 46
prior['A'] = tupak.prior.PowerLaw(alpha=-1, minimum=1e-25, maximum=1e-21, name='A')
prior['f0'] = tupak.prior.Uniform(90, 110, 'f')
47

48
likelihood = tupak.likelihood.GravitationalWaveTransient(IFOs, waveform_generator)
49 50

result = tupak.sampler.run_sampler(
51
    likelihood, prior, sampler='dynesty', outdir=outdir, label=label,
52
    resume=False, sample='unif', injection_parameters=injection_parameters)
53 54 55 56
result.plot_walks()
result.plot_distributions()
result.plot_corner()
print(result)