Variable sources

The example source is pulsed and has a sweeping PA. Both the pulse and sweep are at a frequency of 10 Hz. The simulated light curve is sin(phase * pi)**2 and the simulated PA sweep is fixed-rate counterclockwise motion, with PD = 50%.

import numpy as np
import leakagelib
from ixpeobssim.irf import load_arf

>>> PyXSPEC is not installed, you will no be able to use it.

We’ll load the data as usual, and cut to within 280 arcsec of the center

source = leakagelib.Source.no_image(False)
datas = [leakagelib.IXPEData(source, (
    "data/pulse/event_l2/ixpepulse_det1_evt2_v00.fits",
    "data/pulse/hk/ixpepulse_det1_att_v00.fits",
), energy_cut=(2,8))]

for data in datas:
    data.iterative_centroid_center()
    data.retain(np.sqrt(data.evt_xs**2 + data.evt_ys**2) < 280)

However, we plan to implement phase weights. LeakageLib does not automatically record phases. You should calculate them yourself (e.g. with PINT), then replace the evt_times field with these phases. For this mock data, the frequency is exactly 10 Hz, so the phases are just \(10 t\).

for data in datas:
    # It's convenient to overwrite the "times" list with a list of phases. This source sweeps with frequency of 10 Hz, so multiplying by 10 gives the phase
    data.evt_times *= 10
    data.evt_times = np.fmod(data.evt_times, 1)

Now we create the point source and background sources as usual, setting their spectra. For this fit, we’ll fix the background polarization

settings = leakagelib.FitSettings(datas)
settings.add_point_source("src")
settings.fix_flux("src", 1)

settings.add_background("bkg")
settings.fix_qu("bkg", (0, 0))
settings.set_initial_flux("bkg", 1)

arf = load_arf()
settings.set_spectrum("bkg", lambda e: arf(e) * e**-2.5)
settings.set_spectrum("src", lambda e: arf(e) * e**-1.5)

settings.apply_circular_roi(280)

>>> Reading (in memory) /opt/homebrew/anaconda3/lib/python3.12/site-packages/ixpeobssim/caldb/ixpe/gpd/cpf/arf/ixpe_d1_obssim20240101_v013.arf...

To implement phase weights, we need to tell LeakageLib the source’s light curve. Note that normalization doesn’t matter.

Warning

If you cut based on phase e.g. for an on-off fit, you will need to supply the duty_cycle argument in set_lightcurve.

settings.set_lightcurve("src", lambda ph: np.sin(ph * np.pi)**2)

Now we can perform the fit and see what we get.

fitter = leakagelib.Fitter(datas, settings)
print(fitter)
const_result = fitter.fit()
const_result

FITTED PARAMETERS:
Source  Param
src:    q
src:    u
bkg:    f

FIXED PARAMETERS:
Source  Param   Value
src:    f       1
bkg:    q       0
bkg:    u       0

FitResult:
        q (src) = 0.0636 +/- 0.0724
        u (src) = -0.0575 +/- 0.0713
        f (bkg) = 2.1429 +/- 0.0410

Polarization:
        PD (src): 0.0857 +/- 0.0719 (1.2 sig)
        PA (src): -21.0491 deg +/- 23.9862
Likelihood 20065.77004059819, dof 15765
Optimization terminated successfully.

The best-fit polarization degree was quite low. This is because the source’s true PA sweeps, and we modeled it as constant. We need to make a sweeping PA model.

This model will have two parameters: a constant PD and a PA at phase zero. Let’s create those parameters. The add_param function allows us to do this, and give initial fit values and bounds.

settings.add_param("sweep-PD", 0.1, [0, 1])
settings.add_param("sweep-PA", 0, [-np.pi, np.pi])

Now we need to create a function which takes in the phase and gives the expected Q and U polarization, and tell the fitter to use this function.

def model_fn(ph, fit_data, param_array):
    pd = fit_data.param_to_value(param_array, "sweep-PD")
    pa = fit_data.param_to_value(param_array, "sweep-PA")
    q = pd * np.cos(2 * (ph * 2*np.pi + pa))
    u = pd * np.sin(2 * (ph * 2*np.pi + pa))
    return q, u
settings.set_model_fn("src", model_fn)

Note that set_model_fn automatically tells the fitter to no longer use the src q and src u params. When we re-do the fit, it will use the sweep-PD and sweep-PA params instead.

fitter = leakagelib.Fitter(datas, settings)
print(fitter)
sweep_result = fitter.fit()
sweep_result

FITTED PARAMETERS:
Source  Param
bkg:    f
None:   sweep-PD
None:   sweep-PA

FIXED PARAMETERS:
Source  Param   Value
src:    q       0
src:    u       0
src:    f       1
bkg:    q       0
bkg:    u       0

FitResult:
        f (bkg) = 2.1439 +/- 0.0410
        sweep-PD = 0.4620 +/- 0.0703
        sweep-PA = 0.0202 +/- 0.0784

Polarization:
Likelihood 20086.330839012127, dof 15765
Optimization terminated successfully.

The PD of the sweeping fit was higher. The likelihood was also higher, indicating the sweep model provided a better fit. The true source polarization had PD = 0.5 and PA = 0, which the fit agrees with.