IBSEn: IntraBinary Shock Emission n-not-meaning-anything is a Python package for modeling and visualizing observational data of the gamma-ray binariy pulsars, especially PSR B1259-63. This project includes tools to compute intrabinary shock (IBS) emissivity: spectra and light curves, and (somewhere in future) fit theoretical models to light curves and spectra.
Clone and install manually. Normal installation:
git clone https://github.com/Alvkuzin/IBSEn
cd IBSEn
pip install .For the installation and contributing, replace the last line with
pip install -e .After the installation, install the required packages with
pip install -r requirements.txtThis project requires:
- Python 3 (please install yourself),
- Naima,
- And some standard scientific and graphic libraries. Install them, including Naima, with:
pip install -r requirements.txtNote on Naima. If you have both Naima and numpy installed and they are not in conflict, you are good to go. It is a good idea to create a separate environment for IBSEn where the numpy<2 and Naima will be installed.
Scripts now are mainly not suited for running from the command line... And I don't have proper tests yet, so to find out if the installation works, try running a very basic python script that simply initializes a lot of classes with more or less default parameters and stores the output in a file:
python test_ibsen.py --testall True --ndim -1 > ibsen_test_out.txt
Compare the output with the file ibsen_test_results_0_5_5.txt
There is a poor attempt at the graphical interface: run it with
ibsenYou can explore how the IBS, SEDs, and light curves change for different parameters there, as well as display some data together with a model, and perform a simple normalization fitting to SEDs / LCs.
Mainly, though, the package is meant to be ran as a part of a python script.
The package consists of six main classes, you can find their description and usage examples in tutorials. The general idea is this: you define
an orbit of the binary, then the outflows ("winds") properties, then the intrabinary shock (IBS), then compute the electrons evolution/transport over IBS,
then calculate the photon spectrum. There is also a possibility of doing all these steps at once by calculating a light curve.
Orbit(self-explanatory); Let's initiate an Orbit object and plot an orbit:
from ibsen import Orbit
import matplotlib.pyplot as plt
import numpy as np
DAY = 86400
orb = Orbit(T=25*DAY, e=0.7, M=30*2e33,
nu_los=90*np.pi/180, incl_los=20*np.pi/180)
t = np.linspace(-20*DAY, 70*DAY, 1000)
plt.plot(orb.x(t), orb.y(t))Windsin which currently all information about NS, optical star, and their winds are stored. Here you can calculate the magnetic/photon fields in any point, winds pressure, the position of the equilibrium between pulsar and optical stars winds as a function of time, etc. InitiateWindswith a decretion disk pressure 100 times stronger than the polar wind (see tutorials for further info) and plot the star-to-emission zone distance VS time.
from ibsen import Pulsar, OpticalStar, Winds
star = OpticalStar(f_d=100, Ropt=7e11, Mopt=28*2e33, Topt=4e4)
pulsar = Pulsar(b_ref=10, r_b_ref=star.Ropt, r_p_ref=star.Ropt)
winds = Winds(orbit=orb, star=star, pulsar=pulsar)
t1 = np.linspace(-3*DAY, 3*DAY, 1000)
plt.plot(t1/DAY, winds.dist_se_1d(t1))
td1, td2 = winds.times_of_disk_passage
plt.axvline(td1/DAY)
plt.axvline(td2/DAY)IBS/IBS3D: Intrabinary shock - an object with stuff about IBS geometry and the bulk motion along it; Initiate IBS at a time of 2 days after the periastron and take a look at it, colorcoding it with the doppler factor of matter bulk motion.
from ibsen import IBS
ibs = IBS(winds=winds, t_to_calculate_beta_eff=2*DAY)
ibs.peek(showtime=(-3*DAY, 3*DAY), show_winds=True, ibs_color='doppler')ElectronsOnIBSdescribes the population of ultra-relativistic electrons on the IBS, allows to calculate stationary e-spectra in each point of IBS; Take a look at the electron spectrum over IBS if the apex magnetic field is 10 [G] (was specified inWinds):
from ibsen import ElectronsOnIBS
elev = ElectronsOnIBS(ibs=ibs, cooling='stat_mimic')
elev.calculate()
elev.peek()SpectrumIBScomputes the non-thermal photon spectrum emitted by the ultra-relativistic electrons; Calculate the synchrotron + inverse Compton spectrum from the population of electrons we just found:
from ibsen import SpectrumIBS
spec = SpectrumIBS(els=elev, method='simple', mechanisms=['syn', 'ic',], distance=3e3*3e18)
spec.calculate(e_ph = np.geomspace(3e2, 1e13, 1001))
spec.peek()LightCurveperforms the spectrum calculation for a number of moments of time, returning the light curve and all intermediate calculation results.
from ibsen import LightCurve
t_lc = np.linspace(-3*DAY, 3*DAY, 100)
lc = LightCurve(times = t_lc,
to_parall=True, n_cores=4,
T=25*DAY, e=0.7, M=30*2e33,
nu_los=90*np.pi/180, incl_los=20*np.pi/180,
Ropt=7e11, Mopt=28*2e33, Topt=4e4,
distance=3e3*3e18,
bands = ([300, 1e4], ), cooling='stat_mimic',
f_d=100,
puls_b_ref=10, puls_r_ref=star.Ropt,
method='simple', mechanisms=['syn', 'ic'])
lc.calculate()
lc.peek()See tutorials in tutorials folder for the complete description of these models.
- How to visualize winds in 3D?
- Write fitting utils for LC/SEDs. It should take several datasets and fit to the theoretical model using the same sets of parameters except normalizations.
- If one is bored so much that one feels the need to make the graphical interface better, it'd be cool to add all gamma-ray binaies systems in the drop-down windows, as well as allow for the basic system parameters variation (such as Torb, Mopt, e...) instead of keeping them strictly fixed. Also, the default values and ranges of sliders should be system-dependent.
- The spectrum calculation utils are currently very close to being a general function for calculating the non-thermal spectra from an extended emission zone. It'd be a good idea to actually make it so.
- Do something with occasional NaNs in IC spectrum.