Observation simulations

Simulating NenuFAR observations are permitted thanks to the HpxSimu class which inherits from the HpxSky class, that enables sky coordinates handling using HEALPix representation.

>>> from nenupy.simulation import HpxSimu

In order to keep control of the parameter units, it is recommanded to use Quantity, Time and TimeDelta objects:

>>> from astropy.time import Time, TimeDelta
>>> import astropy.units as u

Besides HpxSky, simulations make use of HpxGSM class that wraps the PyGSM module (see Oliveira-Costa et. al., 2008) while using all functionalities of HpxSky, namely finding out the visible portion of the sky at a given time. This is intended to build a skymodel at NenuFAR frequencies (10 - 80 MHz). HpxSimu computes the NenuFAR beam corresponding to the desired array configuration, taking into account instrumental effects, and multiply it with the GSM skymodel.

Note

Simulation computation time depends on several parameters. In order to optimize it, here are some recommended configurations:

  • SST simulation: resolution~1deg, dt~10min

  • BST simumation: resolution~0.1deg, dt~1min

Warning

Although most intrumental effects are taken into consideration, some still remains to be understood (such as the pointing shift) and some are difficult to take into account (e.g., individual component misbehaviour). Besides, the GSM is not perfect at these low-frequencies (especially regarding point-source apparent size which can be wider than NenuFAR’s beam). Therefore, these simulations should always be interpreted in light of the aforementioned and underlying caveats.

Transit observation

The following demonstrates how to simulate a transit observation (i.e., NenuFAR is pointing towards the same local coordinates for the whole observation duration) occuring at transit_time:

>>> transit_time = Time('2020-04-01 12:00:00')

from two distinct requirements:

From local transit coordinates

Transit observation simulations can imitate NenuFAR array pointing towards a local sky coordinates and therefore needs to be provided with an AltAz object. It can be easily produced with the ho_coord() function:

>>> from nenupy.astro import ho_coord
>>> transit_altaz = ho_coord(
        az=180 * u.deg,
        alt=90 * u.deg,
        time=transit_time
    )

transit_altaz then corresponds to the local zenith.

Setting up NenuFAR array with only one Mini-Array (namely 00), to simulate a SST obervation at 40 MHz with a sky resolution of 1deg and other keywords passed to HpxSimu:

>>> simu = HpxSimu(
        freq=40*u.MHz,
        resolution=1*u.deg,
        ma=0,
        polar='NW'
    )

The call of azel_transit() allows for the simulation to start. dt and duration keywords stand respectively for simulation time step and total duration of the simulation (t0=transit_time being the mid-time of the simulation).

>>> result = simu.azel_transit(
        acoord=transit_altaz,
        t0=transit_time,
        dt=TimeDelta(15*60, format='sec'),
        duration=TimeDelta(24*3600, format='sec'),
    )

simulates a transit observation of NenuFAR containing a single Mini-Array pointing towards the local zenith for a whole day. Simulation is here stored in the result variable which is a SData object having some properties helpful for analysing and plotting purposes:

>>> import matplotlib.pyplot as plt
>>> plt.plot(
        result.datetime,
        result.db,
        label='Transit simulation'
    )
>>> plt.legend()
>>> plt.title('One Mini-Array')
>>> plt.ylabel('dB')
>>> plt.xlabel('Time since {}'.format(result.time[0].isot))
_images/sst_simulation.png

From equatorial transit coordinates

Rather than using AltAz, it might be of interest to simulate a transit of a given astrophysical source (i.e., in equatorial coordinates with ICRS) at a particular time (here transit_time), which can be eased by eq_coord():

>>> from nenupy.astro import eq_coord
>>> cyga_radec = eq_coord(
        ra=299.868*u.deg,
        dec=40.733*u.deg
    )

cyga_radec is an instance of ICRS for the Cygnus A source position.

Simulation is initialized as previously shown, but for the Mini-Arrays parameter ma set to a list of 56 components (rather than one) implying the full (as of April 2020) NenuFAR array. The beam is significantly smaller in this configuration compared to the single Mini-Array beam, therefore a smaller resolution is also preferred:

>>> simu = HpxSimu(
        freq=40*u.MHz,
        resolution=0.5*u.deg,
        ma=np.arange(56),
        polar='NW'
    )

The method radec_transit() is used to compute the simulation on a total duration of 30 minutes with a time step dt of 1 min:

>>> result = simu.radec_transit(
        acoord=cyga_radec,
        t0=transit_time,
        dt=TimeDelta(60, format='sec'),
        duration=TimeDelta(30*60, format='sec'),
    )

Simulation result is a SData object and can be shown along with a vertical line indicating the requested transit time of Cyg A.

>>> import matplotlib.pyplot as plt
>>> plt.plot(
        result.datetime,
        result.db,
        label='Transit simulation'
    )
>>> plt.axvline(
        transit_time.datetime,
        color='black',
        linestyle=':'
    )
>>> plt.legend()
>>> plt.title('56 Mini-Arrays')
>>> plt.ylabel('dB')
>>> plt.xlabel('Time since {}'.format(result.time[0].isot))
_images/simulation_radectransit.png

Tracking observation

Simulating a tracking observation (i.e., fixed pointing in equatorial coordinates) is fairly similar to the aforementioned simulation of a transit observation (see above for a more detailed description regarding the initialization).

Once the HpxSimu instance is properly set up, the radec_tracking() method allows for simulating a tracking. Here, Cygnus A cyga_radec is followed since 2020-04-01 07:00:00 and for 30 minutes with a time step of 10 sec:

>>> result = simu.radec_tracking(
        acoord=cyga_radec,
        t0=Time('2020-04-01 07:00:00'),
        dt=TimeDelta(10, format='sec'),
        duration=TimeDelta(30*60, format='sec'),
    )

Again, result is a SData object, and the simulation profile can be shown:

>>> plt.plot(
        result.datetime,
        result.db,
        label='Tracking simulation'
    )
>>> plt.legend()
>>> plt.title('56 Mini-Arrays')
>>> plt.ylabel('dB')
>>> plt.xlabel('Time since {}'.format(result.time[0].isot))
_images/simulation_radectracking.png

Note

Irregularities in the time profile of the Cygnus A tracking is due to the discrete grid of analog-pointable positions (see NenuFAR website Figure 9).

Simulating BST observations

There is a need in simulating real observations, taking into account actual instrument pointing, in order to assess the correct behaviour of NenuFAR. If NenuFAR BST observation files are available, there is the possiblity to simulate with HpxSimu the contained observations as they occurred, since every instrument and observation configuration metadata are included within the BST FITS files.

A BST observation is loaded using the BST_Data class and, in the following, instanciated as the bst variable:

>>> from nenupy.beamlet import BST_Data
>>> bst = BST_Data('20200219_175900_BST.fits')
>>> data = bst.select(
        freqrange=60*u.MHz
    )

The file '20200219_175900_BST.fits' contains a Taurus A tracking observation. The select() method allows for data selection. Here, in order to simulate a time-profile, a single frequency is selected (the BST_Data object is updated while performing any data selection). The data selected are stored in the data variable, a SData object as well.

The classmethod from_bst from HpxSimu takes as input a BST_Data object (here bst) and the current selection are read in order to perform the simulation:

>>> result = HpxSimu.from_bst(
        bstdata=bst,
        resolution=0.2*u.deg,
        dt=TimeDelta(60, format='sec')
    )

In order to compare simulation result with real BST data, the simulation needs to be scaled as it is not calibrated, this can be done, in a quick manner by:

>>> scale = np.median(data.db) / np.median(result.db)

And here is the result:

>>> plt.plot(
        data.datetime,
        data.db,
        label='Observation'
    )
>>> plt.plot(
        result.datetime,
        result.db * scale,
        label='Simulation'
    )
>>> plt.ylabel('dB')
>>> plt.xlabel('Time since {}'.format(result.time[0].isot))
>>> plt.legend()
>>> plt.axvline(
        (dd.time[0] + TimeDelta(6*60, format='sec')).datetime,
        color='black',
        linestyle=':'
    )
>>> plt.axvline(
        (dd.time[0] + TimeDelta(12*60, format='sec')).datetime,
        color='black',
        linestyle=':'
    )
_images/simulation_frombst.png

Note

The simulation is of course much smoother than the observation which is affected by noise and also by a few tenths of dB variation due to analog pointing switchs, highlighted in the figure by the vertical dotted black lines (see also Identified features under study).