XRB Lightcurve Simulation - Python Version

This is a Python implementation of the XRB (X-Ray Binary) Lightcurve simulation, migrated from the original R code. The simulation calculates column densities for eclipsing binary systems as the compact object eclipses its companion star.

Features

• Vectorized Operations: Uses NumPy for efficient array operations instead of nested loops
• Command Line Interface: Full argparse support with all parameters configurable
• Flexible Output: Configurable output file with default naming
• Modular Design: Clean, well-documented functions for each component
• Type Hints: Full type annotations for better code maintainability

Installation

1. Ensure you have Python 3.7+ installed
2. Install dependencies:

   pip install -r requirements.txt

Usage

Basic Usage (Default Parameters)

python xrb_lightcurve.py

This will run the simulation with default parameters and save results to xrb_lightcurve_output.csv.

Custom Parameters

python xrb_lightcurve.py --r 0.001 --R 2.0 --d1 11.0 --d2 8.0 --gma0 -90.0 --i0 26.0 --dth 1.0 --output my_results.csv

Available Parameters

Parameter	Default	Description
--r	0.001	Radius of smaller star B (compact object) in solar radii
--R	2.0	Radius of larger star A (companion) in solar radii
--d1	11.0	Distance of star B from COM in solar radii
--d2	8.0	Distance of star A from COM in solar radii
--gma0	-90.0	Starting phase angle in degrees
--i0	26.0	Orbital inclination in degrees
--dth	1.0	Orbital increment in degrees
--d2h	6.0	Angular cell size for polar grid (degrees)
--dz	0.1	Step size along line of sight (solar radii)
--n_jobs	1	Number of parallel workers (1 = serial)
--flux_method	legacy	Method for nH to flux conversion (see below)
--flux_csv	None	Path to flux vs nH CSV file (required for interpolate/refit)
--lam	0.589537	Scaling parameter for nH (in 1e22 cm^-2 units)
--lam2	0.589537	Scaling parameter for constant velocity wind
--output	xrb_lightcurve_output.csv	Output file name for results

Flux Conversion Methods

The simulation computes column densities, which are then converted to observable fluxes. Three methods are available:

1. legacy (default): Uses hardcoded exponential fits derived from earlier XSPEC modeling:

• Hard band (2.0-7.0 keV): flux = 9.524 * exp(-0.057 * nH) (legacy coefficients; update via XSPEC CSV for exact band)
  • Soft band (0.3-2 keV): flux = 9.3923 * exp(-2.5062 * nH)
  • Where nH is in units of 1e22 cm^-2

2. interpolate: Directly interpolates from CSV file generated by compute_flux_vs_nH.py:

   • Requires --flux_csv pointing to a CSV with columns: nH_1e22, flux_soft_ph, flux_hard_ph
   • Uses log-log interpolation for smooth behavior
   • Most accurate but requires XSPEC-generated CSV file

3. refit: Fits new exponential functions to CSV data:

   • Requires --flux_csv as above
   • Fits A * exp(-B * nH) to the CSV data
   • Combines accuracy of XSPEC with simplicity of exponential form
   • Prints fitted coefficients for reference

Examples

Basic simulation with default (legacy) flux conversion:

python xrb_lightcurve.py

High-resolution simulation:

python xrb_lightcurve.py --dth 0.5 --output high_res_simulation.csv

Different binary configuration:

python xrb_lightcurve.py --r 0.002 --R 3.0 --d1 15.0 --d2 10.0 --output large_binary.csv

Using XSPEC-computed flux via interpolation:

python xrb_lightcurve.py --flux_method interpolate --flux_csv data_flux_vs_nH.csv --output xspec_interp_results.csv

Fitting new exponentials to XSPEC data:

python xrb_lightcurve.py --flux_method refit --flux_csv data_flux_vs_nH.csv --output xspec_refit_results.csv

Custom nH scaling:

python xrb_lightcurve.py --lam 1.0 --lam2 1.0 --output custom_nh_scale.csv

Output

The simulation generates a CSV file with the following columns:

• deg: Phase angle in degrees
• ph: Phase angle in radians
• phase: Normalized phase (0-1)
• A2: Area calculations
• flx: Column density integral (atoms/solar_radius^4) - accelerated wind
• flx2: Column density integral (atoms/solar_radius^4) - constant velocity wind
• icd: Integrated column density
• time: Time calculations
• l3, L3, h3: Geometric parameters
• fl, fl2: Hydrogen column density nH in units of 1e22 cm^-2
• nfl_hard_av, nfl_hard_cv: Hard band (2.0-7.0 keV) photon fluxes (photons/cm^2/s)
• nfl_soft_av, nfl_soft_cv: Soft band (0.3-2 keV) photon fluxes (photons/cm^2/s)
• pho_count_hard_av, pho_count_soft_av: Photon counts (only in legacy mode)

Understanding Units

The key unit conversions in the simulation:

1. Geometric column density (flx): Computed from the wind integral with units of atoms/(solar radius)^4
2. Scaled nH (fl): fl = flx * lam, where lam is chosen so mean(fl) = 0.589537 by default
   • Units: 1e22 cm^-2
   • Example: fl = 1.0 means nH = 1.0 × 10^22 cm^-2
3. Photon flux: Converted from nH using one of the three methods (legacy/interpolate/refit)
   • Units: photons/cm^2/s
   • Hard band: 2.0-7.0 keV
   • Soft band: 0.3-2 keV

Generating Flux vs nH CSV

To generate a CSV file for use with --flux_method interpolate or --flux_method refit, use the companion script:

python compute_flux_vs_nH.py \
    --specdir ./data/IC10X1_spec \
    --out_csv data_flux_vs_nH.csv \
    --out_png data_flux_vs_nH.png \
    --nH_min 1e20 --nH_max 1e24 --nH_points 60

This requires XSPEC to be available in your Python environment.

Code Structure

Main Functions

1. simulate_lightcurve(): Main simulation function
2. create_grid(): Creates polar grid for wind integral calculations
3. density_function(): Calculates wind density along line of sight
4. wind_los_integral(): Calculates column integral along line of sight

Key Improvements Over R Version

1. Vectorization: Uses NumPy arrays and vectorized operations instead of loops
2. Modularity: Each component is a separate, well-documented function
3. Error Handling: Better handling of edge cases and empty arrays
4. Type Safety: Full type hints for better code maintainability
5. Command Line Interface: Easy parameter configuration via argparse
6. Flexible Output: Configurable output file with default naming

Performance

The Python version is significantly faster than the R version due to:

• Vectorized operations using NumPy
• Reduced nested loops
• More efficient array handling
• Better memory management

Dependencies

• numpy>=1.21.0: For numerical computations and array operations
• pandas>=1.3.0: For data manipulation and CSV output

Troubleshooting

Memory Issues: For very high-resolution simulations (small dth values), consider reducing the grid resolution or using a larger dth value.

Convergence Issues: If the simulation doesn't converge, try adjusting the geometric parameters (r, R, d1, d2) to more physically reasonable values.

Output File Issues: Ensure you have write permissions in the output directory.

License

This code is provided as-is for educational and research purposes.