Ns pulse shape in Vela X-1¶

In this example we will use the function period_sliding_window to obtain both the NS period of the source and observe its evolution and observe the pulse shape. For it we will use a 0.5 s light curve from a XMM-Newton observation.

In [1]:

Copied!

from xraybinaryorbit import *
from xraybinaryorbit import *


HELLO, nice to see you! :)

PLEASE READ THIS, IT'S VERY IMPORTANT:

These are the units that must be used within this package:

- Rstar: Solar radius
- Mstar: Solar masses
- Inclination: Sexagesimal degrees
- Periapsis: Sexagesimal degrees
- Semimajor: Stellar radius
- Periods: Days (Periods in the case of the period_sliding_window function will support any units)
- Iphase: Radians

A list of the functions contained in this package will be displayed by runing the function list_functions().

As these functions use a lot of parameters, which can sometimes be difficult to handle, we have implemented a user-friendly method for parameter input:
A form will be displayed, and the parameters will be saved in the directory for further interactions. These saved parameters will be used if new parameters are not provided.
For the function to work, the submit button must be pressed.
If the parameters are already saved within the working directory, setting "load_directly=True" no form will be displayed and that parameters will be used within the function.
Alternatively, the parameters can be provided as lists in the following format: parameter_list=[parameters] for the theoretical functions or as bound_list = [lower_bounds], [upper_bounds] for the fitting functions.

Please, take into account that fits in general will take A LOT of time to complete.

If you need help, contact graciela.sanjurjo@ua.es.

Load and visualize the high energy range (3-10 keV) and low energy range (0.3-3 keV) light curves¶

In [2]:

Copied!





# Load the light curves
high = pd.read_csv(
    "vela_high.txt",
    comment="!",
    sep=r"\s+",
    skiprows=1,
    header=None,
)

low = pd.read_csv(
    "vela_low.txt",
    comment="!",
    sep=r"\s+",
    skiprows=1,
    header=None,
)

# Columns:
# 0 -> count rate
# 1 -> count-rate uncertainty
# 2 -> time

th, ch, sch = preprocess_data(
    high[2].to_numpy(),
    high[0].to_numpy(),
    high[1].to_numpy(),
)

tl, cl, scl = preprocess_data(
    low[2].to_numpy(),
    low[0].to_numpy(),
    low[1].to_numpy(),
)


# Visualize
plt.figure(figsize=(10, 4))

plt.errorbar(
    tl,
    cl,
    yerr=scl,
    fmt=".",
    markersize=3,
    label="0.2–3 keV",
)

plt.errorbar(
    th,
    ch,
    yerr=sch,
    fmt=".",
    markersize=3,
    label="3–10 keV",
)

plt.xlabel("Time (s)")
plt.ylabel("Count rate")
plt.legend()
plt.tight_layout()
plt.show()
# Load the light curves
high = pd.read_csv(
    "vela_high.txt",
    comment="!",
    sep=r"\s+",
    skiprows=1,
    header=None,
)

low = pd.read_csv(
    "vela_low.txt",
    comment="!",
    sep=r"\s+",
    skiprows=1,
    header=None,
)

# Columns:
# 0 -> count rate
# 1 -> count-rate uncertainty
# 2 -> time

th, ch, sch = preprocess_data(
    high[2].to_numpy(),
    high[0].to_numpy(),
    high[1].to_numpy(),
)

tl, cl, scl = preprocess_data(
    low[2].to_numpy(),
    low[0].to_numpy(),
    low[1].to_numpy(),
)


# Visualize
plt.figure(figsize=(10, 4))

plt.errorbar(
    tl,
    cl,
    yerr=scl,
    fmt=".",
    markersize=3,
    label="0.2–3 keV",
)

plt.errorbar(
    th,
    ch,
    yerr=sch,
    fmt=".",
    markersize=3,
    label="3–10 keV",
)

plt.xlabel("Time (s)")
plt.ylabel("Count rate")
plt.legend()
plt.tight_layout()
plt.show()

No description has been provided for this image

We calculate the sliding window in the next cell, we will chose to have the folded pulses to the periods calculated¶

In [3]:

Copied!





## Choose step configuration
window_sec=6000
step_sec=600
snr_pulse=10

resulth, pulsesh = period_sliding_window(
    th,
    ch,
    sch,
    window_sec=window_sec,
    step_sec=step_sec,
    min_period=260,
    max_period=310,
    folded_pulses=True,
    snr_pulse=snr_pulse,
    nbin_pulse=None,
)

resultl, pulsesl = period_sliding_window(
    tl,
    cl,
    scl,
    window_sec=window_sec,
    step_sec=step_sec,
    min_period=260,
    max_period=310,
    folded_pulses=True,
    snr_pulse=snr_pulse,
    nbin_pulse=None,
)
## Choose step configuration
window_sec=6000
step_sec=600
snr_pulse=10

resulth, pulsesh = period_sliding_window(
    th,
    ch,
    sch,
    window_sec=window_sec,
    step_sec=step_sec,
    min_period=260,
    max_period=310,
    folded_pulses=True,
    snr_pulse=snr_pulse,
    nbin_pulse=None,
)

resultl, pulsesl = period_sliding_window(
    tl,
    cl,
    scl,
    window_sec=window_sec,
    step_sec=step_sec,
    min_period=260,
    max_period=310,
    folded_pulses=True,
    snr_pulse=snr_pulse,
    nbin_pulse=None,
)

All periods are returned, even if not reliable, we need to filter them to meet our requirements¶

In [12]:

Copied!





#Requirements
false_alarm_threshold= 1e-10
snr_threshold = 10

###############################
def filter_period_results(
    result,
    pulses,
    period_min,
    period_max,
    max_false_alarm,
    min_periodogram_snr,
):
    """Filter period-search results without failing when no peaks were found."""
    if result is None or result.empty:
        return pd.DataFrame(), {}

    mask = (
        result["Period"].between(period_min, period_max, inclusive="neither")
        & (result["False_alarm"] < max_false_alarm)
        & (result["snr"] >= min_periodogram_snr)
    )

    filtered_result = result.loc[mask].copy()

    filtered_pulses = {
        index: pulses[index]
        for index in filtered_result.index
        if index in pulses
    }

    return filtered_result, filtered_pulses


filt_resulth, filt_pulsesh = filter_period_results(
    resulth,
    pulsesh,
    period_min=282,
    period_max=285,
    max_false_alarm=false_alarm_threshold,
    min_periodogram_snr=snr_threshold,
)

filt_resultl, filt_pulsesl = filter_period_results(
    resultl,
    pulsesl,
    period_min=282,
    period_max=285,
    max_false_alarm=false_alarm_threshold,
    min_periodogram_snr=snr_threshold,
)

filt_resultl.head()
#Requirements
false_alarm_threshold= 1e-10
snr_threshold = 10

###############################
def filter_period_results(
    result,
    pulses,
    period_min,
    period_max,
    max_false_alarm,
    min_periodogram_snr,
):
    """Filter period-search results without failing when no peaks were found."""
    if result is None or result.empty:
        return pd.DataFrame(), {}

    mask = (
        result["Period"].between(period_min, period_max, inclusive="neither")
        & (result["False_alarm"] < max_false_alarm)
        & (result["snr"] >= min_periodogram_snr)
    )

    filtered_result = result.loc[mask].copy()

    filtered_pulses = {
        index: pulses[index]
        for index in filtered_result.index
        if index in pulses
    }

    return filtered_result, filtered_pulses


filt_resulth, filt_pulsesh = filter_period_results(
    resulth,
    pulsesh,
    period_min=282,
    period_max=285,
    max_false_alarm=false_alarm_threshold,
    min_periodogram_snr=snr_threshold,
)

filt_resultl, filt_pulsesl = filter_period_results(
    resultl,
    pulsesl,
    period_min=282,
    period_max=285,
    max_false_alarm=false_alarm_threshold,
    min_periodogram_snr=snr_threshold,
)

filt_resultl.head()

Out[12]:

	window	window_start	window_end	min_time	max_time	Frequency	Frequency_Lower	Frequency_Upper	Freq_Error	Freq_Error_Minus	...	Period_Upper	Period_Error	Period_Error_Minus	Period_Error_Plus	Power	Power_Noise	Power_Error	False_alarm	snr	Delta_Chi2
8	62.0	2.649787e+08	2.649847e+08	2.649787e+08	2.649847e+08	0.003511	0.003505	0.003518	0.000007	0.000007	...	285.310839	0.532775	0.536075	0.529475	0.009378	0.000619	0.000619	3.059815e-20	13.703299	1.0
9	63.0	2.649793e+08	2.649853e+08	2.649793e+08	2.649853e+08	0.003516	0.003510	0.003522	0.000006	0.000006	...	284.874796	0.481358	0.483014	0.479702	0.009338	0.000833	0.000833	3.124064e-20	10.090100	1.0

2 rows × 23 columns

OBSERVE RESULT¶

Rebin light curves for visualization purposes.¶

In [13]:

Copied!





high_reb = rebin_snr(
    th,
    ch,
    sch,
    min_snr=50,
    keep_partial=True,
)

low_reb = rebin_snr(
    tl,
    cl,
    scl,
    min_snr=10,
    keep_partial=True,
)


fig, ax = plt.subplots(figsize=(10, 5))


# ---------------------------------------------------------------------
# Light curves scaled only for visual comparison with the spin periods
# ---------------------------------------------------------------------

if len(high_reb[0]) > 0 and not filt_resulth.empty:
    ax.plot(
        high_reb[0],
        scale(
            high_reb[1],
            np.asarray(filt_resulth["Period"], dtype=float),
        ),
        color="red",
        alpha=0.20,
        label="3–10 keV light curve",
    )

if len(low_reb[0]) > 0 and not filt_resultl.empty:
    ax.plot(
        low_reb[0],
        scale(
            low_reb[1],
            np.asarray(filt_resultl["Period"], dtype=float),
        ),
        color="green",
        alpha=0.20,
        label="0.2–3 keV light curve",
    )


# ---------------------------------------------------------------------
# Helper for plotting periods with symmetric time errors and asymmetric
# period errors
# ---------------------------------------------------------------------

def plot_period_results(
    ax,
    result,
    color,
    label,
    alpha=1.0,
):
    """Plot sliding-window periods with asymmetric period uncertainties."""

    if result is None or result.empty:
        return

    x = (
        np.asarray(result["min_time"], dtype=float)
        + np.asarray(result["max_time"], dtype=float)
    ) / 2.0

    y = np.asarray(result["Period"], dtype=float)

    # Half-width of each temporal window.
    xerr = (
        np.asarray(result["max_time"], dtype=float)
        - np.asarray(result["min_time"], dtype=float)
    ) / 2.0

    # Matplotlib expects asymmetric errors with shape (2, N):
    # first row = lower errors; second row = upper errors.
    yerr = np.vstack(
        [
            np.asarray(
                result["Period_Error_Minus"],
                dtype=float,
            ),
            np.asarray(
                result["Period_Error_Plus"],
                dtype=float,
            ),
        ]
    )

    finite = (
        np.isfinite(x)
        & np.isfinite(y)
        & np.isfinite(xerr)
        & np.isfinite(yerr[0])
        & np.isfinite(yerr[1])
        & (xerr >= 0)
        & (yerr[0] >= 0)
        & (yerr[1] >= 0)
    )

    if not np.any(finite):
        return

    ax.errorbar(
        x[finite],
        y[finite],
        #xerr=xerr[finite],
        yerr=yerr[:, finite],
        fmt="o",
        color=color,
        ecolor=color,
        markerfacecolor=color,
        markeredgecolor=color,
        markersize=5,
        elinewidth=1.2,
        capsize=3,
        alpha=alpha,
        label=label,
        linestyle="none",
        zorder=3,
    )


plot_period_results(
    ax,
    filt_resulth,
    color="blue",
    label="3–10 keV periods",
    alpha=1.0,
)

plot_period_results(
    ax,
    filt_resultl,
    color="green",
    label="0.2–3 keV periods",
    alpha=0.7,
)


# ---------------------------------------------------------------------
# Figure formatting
# ---------------------------------------------------------------------

all_times = np.concatenate(
    [
        np.asarray(th, dtype=float),
        np.asarray(tl, dtype=float),
    ]
)

finite_times = all_times[np.isfinite(all_times)]

if len(finite_times) > 0:
    ax.set_xlim(
        np.min(finite_times),
        np.max(finite_times),
    )

ax.set_xlabel("Time (s)", fontsize=20)
ax.set_ylabel("NS spin period (s)", fontsize=20)

ax.tick_params(
    axis="both",
    which="major",
    labelsize=14,
)

ax.legend()
fig.tight_layout()

fig.savefig(
    "velax1_spin.png",
    dpi=300,
    bbox_inches="tight",
)

plt.show()
high_reb = rebin_snr(
    th,
    ch,
    sch,
    min_snr=50,
    keep_partial=True,
)

low_reb = rebin_snr(
    tl,
    cl,
    scl,
    min_snr=10,
    keep_partial=True,
)


fig, ax = plt.subplots(figsize=(10, 5))


# ---------------------------------------------------------------------
# Light curves scaled only for visual comparison with the spin periods
# ---------------------------------------------------------------------

if len(high_reb[0]) > 0 and not filt_resulth.empty:
    ax.plot(
        high_reb[0],
        scale(
            high_reb[1],
            np.asarray(filt_resulth["Period"], dtype=float),
        ),
        color="red",
        alpha=0.20,
        label="3–10 keV light curve",
    )

if len(low_reb[0]) > 0 and not filt_resultl.empty:
    ax.plot(
        low_reb[0],
        scale(
            low_reb[1],
            np.asarray(filt_resultl["Period"], dtype=float),
        ),
        color="green",
        alpha=0.20,
        label="0.2–3 keV light curve",
    )


# ---------------------------------------------------------------------
# Helper for plotting periods with symmetric time errors and asymmetric
# period errors
# ---------------------------------------------------------------------

def plot_period_results(
    ax,
    result,
    color,
    label,
    alpha=1.0,
):
    """Plot sliding-window periods with asymmetric period uncertainties."""

    if result is None or result.empty:
        return

    x = (
        np.asarray(result["min_time"], dtype=float)
        + np.asarray(result["max_time"], dtype=float)
    ) / 2.0

    y = np.asarray(result["Period"], dtype=float)

    # Half-width of each temporal window.
    xerr = (
        np.asarray(result["max_time"], dtype=float)
        - np.asarray(result["min_time"], dtype=float)
    ) / 2.0

    # Matplotlib expects asymmetric errors with shape (2, N):
    # first row = lower errors; second row = upper errors.
    yerr = np.vstack(
        [
            np.asarray(
                result["Period_Error_Minus"],
                dtype=float,
            ),
            np.asarray(
                result["Period_Error_Plus"],
                dtype=float,
            ),
        ]
    )

    finite = (
        np.isfinite(x)
        & np.isfinite(y)
        & np.isfinite(xerr)
        & np.isfinite(yerr[0])
        & np.isfinite(yerr[1])
        & (xerr >= 0)
        & (yerr[0] >= 0)
        & (yerr[1] >= 0)
    )

    if not np.any(finite):
        return

    ax.errorbar(
        x[finite],
        y[finite],
        #xerr=xerr[finite],
        yerr=yerr[:, finite],
        fmt="o",
        color=color,
        ecolor=color,
        markerfacecolor=color,
        markeredgecolor=color,
        markersize=5,
        elinewidth=1.2,
        capsize=3,
        alpha=alpha,
        label=label,
        linestyle="none",
        zorder=3,
    )


plot_period_results(
    ax,
    filt_resulth,
    color="blue",
    label="3–10 keV periods",
    alpha=1.0,
)

plot_period_results(
    ax,
    filt_resultl,
    color="green",
    label="0.2–3 keV periods",
    alpha=0.7,
)


# ---------------------------------------------------------------------
# Figure formatting
# ---------------------------------------------------------------------

all_times = np.concatenate(
    [
        np.asarray(th, dtype=float),
        np.asarray(tl, dtype=float),
    ]
)

finite_times = all_times[np.isfinite(all_times)]

if len(finite_times) > 0:
    ax.set_xlim(
        np.min(finite_times),
        np.max(finite_times),
    )

ax.set_xlabel("Time (s)", fontsize=20)
ax.set_ylabel("NS spin period (s)", fontsize=20)

ax.tick_params(
    axis="both",
    which="major",
    labelsize=14,
)

ax.legend()
fig.tight_layout()

fig.savefig(
    "velax1_spin.png",
    dpi=300,
    bbox_inches="tight",
)

plt.show()

We observe the filtered high and low energy folded pulses¶

In [9]:

Copied!





# Number of pulses to plot
num_pulses = len(filt_pulsesh)

# Define number of rows and columns
n_cols = 4
n_rows = (num_pulses + n_cols - 1) // n_cols  # Ensure enough rows for all pulses

# Create subplots
fig, axes = plt.subplots(n_rows, n_cols, figsize=(n_cols * 4+1, n_rows * 3+1))

# Flatten axes in case of multiple rows
axes = axes.flatten()

# Loop over the filtered pulses and plot
for i, (idx, pulse_data) in enumerate(filt_pulsesh.items()):
    ax = axes[i]
    ph_pulse = pulse_data['ph_pulse']
    c_pulse = pulse_data['c_pulse']
    sc_pulse = pulse_data['sc_pulse']
    
    concatenated_c_pulse = np.concatenate([c_pulse[np.argmin(c_pulse):len(c_pulse)], c_pulse[0:np.argmin(c_pulse)]])
    concatenated_ph_pulse = np.concatenate([ph_pulse[np.argmin(c_pulse):len(c_pulse)], np.array(ph_pulse[0:np.argmin(c_pulse)]) + 1])
    concatenated_sc_pulse = np.concatenate([sc_pulse[np.argmin(c_pulse):len(c_pulse)], sc_pulse[0:np.argmin(c_pulse)]])
    
    # Plot pulse data
    ax.errorbar(concatenated_ph_pulse, concatenated_c_pulse, yerr=concatenated_sc_pulse, fmt=':', label=f'Pulse {i+1}',alpha=1)
    ax.errorbar(concatenated_ph_pulse, concatenated_c_pulse,  fmt='b')
    ax.set_xlabel('NS Spin Phase')
    ax.set_ylabel('Counts')
    ax.legend()


# Turn off empty subplots if any
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

plt.savefig("high_folded_pulses_velax1.png")
# Number of pulses to plot
num_pulses = len(filt_pulsesh)

# Define number of rows and columns
n_cols = 4
n_rows = (num_pulses + n_cols - 1) // n_cols  # Ensure enough rows for all pulses

# Create subplots
fig, axes = plt.subplots(n_rows, n_cols, figsize=(n_cols * 4+1, n_rows * 3+1))

# Flatten axes in case of multiple rows
axes = axes.flatten()

# Loop over the filtered pulses and plot
for i, (idx, pulse_data) in enumerate(filt_pulsesh.items()):
    ax = axes[i]
    ph_pulse = pulse_data['ph_pulse']
    c_pulse = pulse_data['c_pulse']
    sc_pulse = pulse_data['sc_pulse']
    
    concatenated_c_pulse = np.concatenate([c_pulse[np.argmin(c_pulse):len(c_pulse)], c_pulse[0:np.argmin(c_pulse)]])
    concatenated_ph_pulse = np.concatenate([ph_pulse[np.argmin(c_pulse):len(c_pulse)], np.array(ph_pulse[0:np.argmin(c_pulse)]) + 1])
    concatenated_sc_pulse = np.concatenate([sc_pulse[np.argmin(c_pulse):len(c_pulse)], sc_pulse[0:np.argmin(c_pulse)]])
    
    # Plot pulse data
    ax.errorbar(concatenated_ph_pulse, concatenated_c_pulse, yerr=concatenated_sc_pulse, fmt=':', label=f'Pulse {i+1}',alpha=1)
    ax.errorbar(concatenated_ph_pulse, concatenated_c_pulse,  fmt='b')
    ax.set_xlabel('NS Spin Phase')
    ax.set_ylabel('Counts')
    ax.legend()


# Turn off empty subplots if any
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

plt.savefig("high_folded_pulses_velax1.png")

In [10]:

Copied!





# Number of pulses to plot
num_pulses = len(filt_pulsesl)

# Define number of rows and columns
n_cols = 4
n_rows = (num_pulses + n_cols - 1) // n_cols  # Ensure enough rows for all pulses

# Create subplots
fig, axes = plt.subplots(n_rows, n_cols, figsize=(n_cols * 4+1, n_rows * 3+1))

# Flatten axes in case of multiple rows
axes = axes.flatten()

# Loop over the filtered pulses and plot
for i, (idx, pulse_data) in enumerate(filt_pulsesl.items()):
    ax = axes[i]
    ph_pulse = pulse_data['ph_pulse']
    c_pulse = pulse_data['c_pulse']
    sc_pulse = pulse_data['sc_pulse']
    
    concatenated_c_pulse = np.concatenate([c_pulse[np.argmin(c_pulse):len(c_pulse)], c_pulse[0:np.argmin(c_pulse)]])
    concatenated_ph_pulse = np.concatenate([ph_pulse[np.argmin(c_pulse):len(c_pulse)], np.array(ph_pulse[0:np.argmin(c_pulse)]) + 1])
    concatenated_sc_pulse = np.concatenate([sc_pulse[np.argmin(c_pulse):len(c_pulse)], sc_pulse[0:np.argmin(c_pulse)]])
    
    # Plot pulse data
    ax.errorbar(concatenated_ph_pulse, concatenated_c_pulse, yerr=concatenated_sc_pulse, fmt=':', label=f'Pulse {i+1}',alpha=0.2)
    ax.errorbar(concatenated_ph_pulse, concatenated_c_pulse,  fmt='g')
    ax.set_xlabel('NS Spin Phase')
    ax.set_ylabel('Counts')
    ax.legend()


# Turn off empty subplots if any
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

plt.savefig("low_folded_pulses_velax1.png")
# Number of pulses to plot
num_pulses = len(filt_pulsesl)

# Define number of rows and columns
n_cols = 4
n_rows = (num_pulses + n_cols - 1) // n_cols  # Ensure enough rows for all pulses

# Create subplots
fig, axes = plt.subplots(n_rows, n_cols, figsize=(n_cols * 4+1, n_rows * 3+1))

# Flatten axes in case of multiple rows
axes = axes.flatten()

# Loop over the filtered pulses and plot
for i, (idx, pulse_data) in enumerate(filt_pulsesl.items()):
    ax = axes[i]
    ph_pulse = pulse_data['ph_pulse']
    c_pulse = pulse_data['c_pulse']
    sc_pulse = pulse_data['sc_pulse']
    
    concatenated_c_pulse = np.concatenate([c_pulse[np.argmin(c_pulse):len(c_pulse)], c_pulse[0:np.argmin(c_pulse)]])
    concatenated_ph_pulse = np.concatenate([ph_pulse[np.argmin(c_pulse):len(c_pulse)], np.array(ph_pulse[0:np.argmin(c_pulse)]) + 1])
    concatenated_sc_pulse = np.concatenate([sc_pulse[np.argmin(c_pulse):len(c_pulse)], sc_pulse[0:np.argmin(c_pulse)]])
    
    # Plot pulse data
    ax.errorbar(concatenated_ph_pulse, concatenated_c_pulse, yerr=concatenated_sc_pulse, fmt=':', label=f'Pulse {i+1}',alpha=0.2)
    ax.errorbar(concatenated_ph_pulse, concatenated_c_pulse,  fmt='g')
    ax.set_xlabel('NS Spin Phase')
    ax.set_ylabel('Counts')
    ax.legend()


# Turn off empty subplots if any
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

plt.savefig("low_folded_pulses_velax1.png")

In [ ]: