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.