Documentation for pyLCSIM
**************************
============
Introduction
============
pyLCSIM is a python package to simulate X-ray lightcurves from coherent signals and power spectrum models.
Coherent signals can be specified as a sum of one or more sinusoids, each with its frequency, pulsed fraction and phase shift; or as a series of harmonics of a fundamental frequency (each with its pulsed fraction and phase shift).
Power spectra can be simulated from a model of the power spectrum density (PSD), using as a template one or more of the built-in library functions. The user can also define his/her custom models.
Models are additive.
A PDF version of these notes is available `here <../latex/pyLCSIM.pdf>`_.
**Warning**: the current release (0.x.y) is HIGHLY EXPERIMENTAL! Use at your own risk...
=============
Prerequisites
=============
pyLCSIM requires `Numpy `_ (at least v1.8) and `Astropy `_ (at least v0.3).
`Matplotlib `_ is highly recommended if you want to plot your simulations.
============
Installation
============
The most straightforward way to install pyLCSIM is to use the `Python Package Index `_ and its pip utility (may require administrator privileges)::
$ pip install pyLCSIM
To upgrade from a previous version::
$ pip install pyLCSIM --upgrade
The package source can be also downloaded `here <../pyLCSIM-0.4.0.tar.gz>`_.
The installation in this case follows the usual steps::
$ tar xzvf pyLCSIM-0.x.y.tar.gz
$ cd pyLCSIM-0.x.y
$ python setup.py install
The last step may require administrator privileges.
=============
Mailing list
=============
A dedicated mailing list is available, to host announcements of new releases, to report issues and to get support.
If interested, subscribe `here `_.
=============
License
=============
The package is distributed under the `MIT license `_.
============
Example 1
============
Let's begin with a PSD model simulation.
We import the usual packages::
import matplotlib.pyplot as plt
import numpy as np
import pyLCSIM
We assume that our source has a rate of 30000 counts/s. Moreover, we have a 5000 counts/s background rate, and we have made a 50 s exposure. Our observation has a time resolution of 10 ms.
We want to simulate a QPO at a frequency of 10 Hz, superimposed to a continuum modelled as a smoothly-varying broken power law (spectral indices 1 and 2, with a steepness change at 1 Hz).
The required fractional RMS variation of the signal is 1%::
rate_src = 30000.0
rate_bkg = 5000.0
t_exp = 50.0
dt = 0.01
frms = 0.1
The total bins of the lightcurve are therefore::
nbins = t_exp/dt
The simulation follows as::
# Instantiate a simulation object
sim = pyLCSIM.Simulation()
# Add two PSD models: a smooth broken power law and a Lorentzian representing a QPO.
# See the documentation for details.
sim.addModel('smoothbknpo', [1., 1, 2, 1])
sim.addModel('lorentzian', [10., 1., 10, 2])
# Run the simulation
sim.run(dt, nbins, rate_src, rms=frms)
# Add Poisson noise to the light curve
sim.poissonRandomize(dt, rate_bkg)
# Get lightcurve and power spectrum as 1-D arrays
time, rate = sim.getLightCurve()
f, psd = sim.getPowerSpectrum()
Done! We can save the results as FITS files::
# Save FITS files with lightcurve and spectrum
pyLCSIM.saveFITSLC("myLC.fits", time, rate)
pyLCSIM.saveFITSPSD("myPSD.fits", f, psd)
and view the results::
# Plot the lightcurve and power spectrum
fig0 = plt.figure()
plt.plot(time, rate)
plt.xlabel("Time [s]")
plt.ylabel("Rate [counts/s]")
plt.title("Lightcurve")
fig1 = plt.figure()
plt.loglog(f, psd, drawstyle='steps-mid', color='black')
plt.xlabel("Frequency [Hz]")
plt.ylabel("PSD [Leahy normalized]")
plt.title("Power spectrum density")
plt.show()
.. image:: img/ex1_lc.png
.. image:: img/ex1_psd.png
The PSD model used in the simulation can be obtained (and plotted) following the next example code::
freq_m, model_tot, model_comp = sim.getPSDModel(dt, nbins)
plt.figure()
for mod in model_comp:
# Plot the various additive component as dashed lines
plt.loglog(freq_m, mod, ls='dashed', color='red')
# Plot the total model
plt.loglog(freq_m, model_tot)
plt.show()
.. image:: img/ex1_psdmodel.png
============
Example 2
============
The user can also define his/her own PSD models, simply as python functions.
The only caveat is that the function should be positive-valued, and it is suggested to avoid too small output values.
The following example shows a simulation using an user-defined function (a Gaussian centered at 10 Hz, in this case)::
def myFunc(f, p):
"""
Example of user-defined function: a Gaussian.
User-defined PSD models should be positive-valued!
Moreover, in this example the output is clipped at 1e-32 to avoid too small values.
"""
f = p[0]*np.exp(-(f-p[1])**2/p[2]**2)
return np.clip(f, 1e-32, np.max(f))
sim = pyLCSIM.Simulation()
sim.addModel('smoothbknpo', [1., 1, 2, 1])
sim.addModel(myFunc, [1000., 10, 1.])
# Run the simulation
sim.run(dt, nbins, rate_src, rms=frms)
.. image:: img/ex2_psd.png
Additionally, PSD models can be also defined as `astropy.modeling.models `_ functions.
The following example is equivalent to the previous one (a Gaussian centered at 10 Hz)::
from astropy.modeling import models
sim = pyLCSIM.Simulation()
sim.addModel('smoothbknpo', [1., 1, 2, 1])
# For astropy.modeling.models, you should define a dictionary of parameters instead of the usual list:
sim.addModel(models.Gaussian1D, {'amplitude':1000., 'mean':10., 'stddev':1.})
# Run the simulation
sim.run(dt, nbins, rate_src, rms=frms)
============
Example 3
============
The following example shows the simulation of a coherent signal
using a sum of sinusoids::
import matplotlib.pyplot as plt
import numpy as np
import pyLCSIM
rate_src = 30000.0
rate_bkg = 5000.0
t_exp = 1.0
dt = 0.0001
nbins = t_exp/dt
Note the different exposure time (1 s) and time resolution (100 us)::
# Instantiate a simulation object, this time as coherent
sim = pyLCSIM.Simulation(kind='coherent')
# Run the simulation, using:
# four sinusoidal frequencies: 340, 550, 883, 1032 Hz;
# with pulsed fractions 10%, 5%, 7% and 15% respectively;
# the third frequency has a 35 degree phase shift with respect to the others
sim.run(dt, nbins, rate_src, freq=[340, 550, 883, 1032], amp=[0.1, 0.05, 0.07, 0.15], phi=[0., 0, 35., 0.])
# Add Poisson noise to the light curve
sim.poissonRandomize(dt, rate_bkg)
# Get lightcurve and power spectrum as 1-D arrays
time, rate = sim.getLightCurve()
f, psd = sim.getPowerSpectrum()
# Plot the lightcurve and power spectrum
fig0 = plt.figure()
plt.plot(time, rate)
plt.xlabel("Time [s]")
plt.ylabel("Rate [counts/s]")
plt.title("Lightcurve")
fig1 = plt.figure()
plt.loglog(f, psd, drawstyle='steps-mid', color='black')
plt.xlabel("Frequency [Hz]")
plt.ylabel("PSD [Leahy normalized]")
plt.title("Power spectrum density")
# Save FITS files with lightcurve and spectrum
pyLCSIM.saveFITSLC("myLC.fits", time, rate)
pyLCSIM.saveFITSPSD("myPSD.fits", f, psd)
plt.show()
.. image:: img/ex3_psd.png
============
Example 4
============
Finally, an example with a fundamental frequency and two harmonics::
import matplotlib.pyplot as plt
import numpy as np
import pyLCSIM
rate_src = 300000.0
rate_bkg = 3000.0
t_exp = 1.0
dt = 0.0001
phase_shift = 0.
nbins = t_exp/dt
# Instantiate a simulation object, this time as coherent
sim = pyLCSIM.Simulation(kind='coherent')
# Run the simulation:
# Fundamental at 500 Hz, 3 harmonics (500, 1000, 1500 Hz)
# with pulsed fractions 10%, 5% and 15% respectively
sim.run(dt, nbins, rate_src, freq=500, nha=3, amp=[0.1, 0.05, 0.15])
# Add Poisson noise to the light curve
sim.poissonRandomize(dt, rate_bkg)
# Get lightcurve and power spectrum as 1-D arrays
time, rate = sim.getLightCurve()
f, psd = sim.getPowerSpectrum()
# Plot the lightcurve and power spectrum
fig0 = plt.figure()
plt.plot(time, rate)
fig1 = plt.figure()
plt.loglog(f, psd, drawstyle='steps-mid', color='black')
# Save FITS files with lightcurve and spectrum
pyLCSIM.saveFITSLC("myLC.fits", time, rate)
pyLCSIM.saveFITSPSD("myPSD.fits", f, psd)
plt.show()
===========
Main module
===========
.. automodule:: pyLCSIM
:members:
=====================
Submodule: psd_models
=====================
.. automodule:: pyLCSIM.psd_models
:members:
=====================
Changelog
=====================
v0.4.0: Added support for astropy.modeling.models functions.
v0.3.0: Added python 3.x compatibility.
v0.2.3: Minor bugfixes.
v0.2.2: Modified rebinning functions (rebin() and logrebin()).
v0.2.1: Added getPSDModel() method to Simulation. Bugfixes.
v0.2.0: Added the possibility to employ user-defined PSD models.
v0.1.2: Added reset method to Simulation. Lightcurve FITS output is now OGIP-compliant.
v0.1.1: Bugfix. Modified Simulation class.
v0.1.0: Initial release.