Current state of CLAR simulations
Current scripts reside under CVS in LOFAR/Timba/WH/contrib/AGW.
Creating an MS
We simulate an observation with a VLAx10 (VLA positions multiplied by 10, so it's a big VLA in outer space) composed of CLAR dishes. This contains 27 antennas, so 351 baselines. Full UV coverage requires an 8 hour observation. We use 32 channels of 19 MHz width, from 800MHz to 1408Mhz.
First you need to run ska_sim.g to create an empty MS. See inside for details. Otherwise, you can just copy an MS from my home on birch (.../contrib/AGW). TEST_CLAR_27-480.MS is the smaller one, using 60-second integrations (480 timeslots total). TEST_CLAR_27-4800.MS uses 6-second integrations. I recommend using the smaller MS unless you really want to measure performance under heavy load.
Fitting derived quantities
It is quite expensive (and non-parallelizable) to use AIPS++ to compute the miriad az/el points required for simulating the CLAR beam. We therefore have a separate tree that precomputes the beam parameters and fits them with a polynomial. For a further speedup, it also fits the UVWs. It is necessary to run this script once before doing any further simulations.
Run clar_fit_dq.py and execute the only job in it. You can look at bookmarks to see the fitted parameters. This produces a MEP table called CLAR_DQ_27-480.mep. This MEP table is then used by the predict and solve trees. Note that the 480 timeslot table is still good for the 4800 timeslot MS.
Simulating an observation
clar_fast_predict.py fills the DATA column of the MS with simulated data. It is controlled by a number of parameters at the top of the script.
Imaging
make_image.g will image the DATA column and save the result in clar_data.img. See inside for details.
Calibrating a simulated observation
clar_fast_solve.py attempts to calibrate the simulated observation, solving for source fluxes, spectral indices, and the beam width parameter. Make sure you use the same source model (point source or extended) as you did for the predict step.
The calibration residuals are written to the CORRECTED_DATA column of the MS.
At the moment we can successfully solve for point sources and beam width. For extended sources, the solution does not converge so well. This remains to be investigated.
Simulating error effects
clar_predict_diff.py can be used to simulate two sets of visibilities with different instrumental effects, and write the result to the DATA column of the MS. Imaging this gives you an idea of the image-plane errors introduced by lack of calibration for a specific effect. See discussion on results below.
Initial results
Simulating extended sources
I have added some classes to model elliptical gaussian sources. These can be activated by setting the source_model variable at the top of clar_fast_predict.py. The following is a dirty image with a few extended sources in it. Note the source in the lower-left of the image is a combination of a point source, and a gaussian under it.
The very large fractional bandwidth we're using means that extended sources fade away sharply as we go from 800MHz to 1408Mhz, due to the effective baselines (in wavelengths) getting longer. To compensate for this (so that we still see something at the high frequencies), I have assigned them large positive spectral indices so that their intrinsic fluxes go up with frequency.
For comparison, here's a dirty image of a point source-only model.
Investigating effects of model errors
Beam squishing
Our CLAR model incorporates the beam squishing effect, which makes the station beam "squish" at lower elevations. At zenith the beam is circular; towards the horizon it becomes more and more of an ellipse with a vertical major axis. This is due to the dish being at an oblique angle relative to the pointing direction, decreasing its apparent "vertical" size.
The consequence of this is that each source will have a time-variable gain, in principle different for every station (due to the source being seen at different elevations by stations far enough apart). The clar_predict_diff.py scripts lets us see how this affect our data. We use it to compute the difference between a "true" model of 10 point sources (see dirty map above) that takes beam squishing into account, and a "simplified" model which ues the same circular beam regardless of elevation. Here's a dirty map of the difference at 800MHz:
This map has a number of interesting features:
- Sources are distorted to a different degree, with a characteristic 6-spoke pattern. This pattern corresponds to the instantaneous PSF of the VLA configuration. Basically, as the source goes up and down in elevation, it is observed with a different voltage gain so its apparent flux changes. The spokes are produced when the gain is higher.
- As the sky rotates, the beam "squishes" and rotates w.r.t. to the sky. This means that different sources experience the above "gain distortion" at different times during the 8 hours of observation. This is visible if you compare the spoke pattern of the lower-left and lower-right group of sources: they are not parallel.
- The observation takes 8 hours, so the beam rotates through only part of the sky. The three top-center sources experience relatively little distortion.
- Since the beam becomes smaller at higher frequency, we would expect the distortion to be more pronounced. Indeed, here's a dirty map at 1408MHz. Now the top source has also acquired spokes:
- The peak residual flux is about 0.02Jy (all sources are 1Jy); the spokes are about 0.01Jy.
- I'm still at a loss to explain some other features of this map, so feel free to enlighten me.
Beam averaging
Image-plane-based simulators such as the newsimulator of AIPS++ construct the entire model in the image plane. This makes it impossible (or at least very difficult) to simulate instrumental effects that vary both in direction and from station to station, which MeqTrees can easily model in the uv domain.
Beam squishing is a good example of such an effect. Since each station sees the phase center at a slightly different elevation, each station's beam will be squished to a subtly different extent. In the image plane, however, we can only apply one beam pattern. To repeat our CLAR simulation, the AIPS++ simulator would have to use something that is effectively an "average" beam for all stations and for a given period of time.
Since our longest baseline is only ~30km, the difference is very subtle, but it should be there. I've used clar_predict_diff.py to calculate the difference between a "true" model with individual per-station beams, and an "averaged" model using the same beam for all station (no avergaing in time though). Here's a map of the results at 1408MHz:
The error level is 10-5 Jy, which doesn't sound like much, but remember the target dynamic range for SKA!
A histogram of the differences shows a symmetrical distribution between -10-5 and 10-5 Jy
- The maximum "error" made by the averaged beam should be when the tracking direction is parallel to a long baseline. This probably corresponds to the dips and spokes found around the point sources.