tjgalvin / flint

BSD 3-Clause "New" or "Revised" License
7 stars 1 forks source link

Devlop test for 3C286 data #63

Open tjgalvin opened 7 months ago

tjgalvin commented 7 months ago

Late last year data were taken on a number of known polarised sources as a way of verifying we understood how to correct measurement sets for tools outside of yandasoft and the larger askap pipeline. These data have been useful in verifying the fixms tool and the general continuum calibration and imaging workflow.

In terms of full polarisation, arrakis is starting to show some interesting and difficult to understand results. Using these 3C286 data and components of arrakis it has been highlighted that something is going wrong in the derivation (or application) of the smoothed bandpass.

A proper test should be incorporated into the general pytest suite to verify that the 3C286 data produces a correct stokes I/Q/U spectrum.

The bandpass SBID: 51998 The SBIDs of the science fields: 51997 and 52087

Since 3C286 is bright a minute or two should be sufficent. Bandpass calibration might require some thinking about since the MS are larger.

We would need to store in an accessible way:

AlecThomson commented 7 months ago

For later use, here are the Perley+Bulter reference values:

In Stokes I they have 4-paramter model, fit to log(freq/GHz). The model values can be obtained by doing something like

PB_scales = {
    0: 1.2481,
    1: -0.4507,
    2: -0.1798,
    3:  0.0357,
}
nu_g_log = np.log10(nu_g) # nu_g are frequencies in GHz
pb_model_log = PB_scales[0] + \
            PB_scales[1] * nu_g_log + \
            PB_scales[2] * nu_g_log**2 + \
            PB_scales[3] * nu_g_log**3
pb_model_i = 10**(pb_model_log) # These are the model Stokes I fluxes in Jy

For polarisation, we have reference values at a few frequencies in terms of polarisation fraction and angle. These values are:

pb_pol = {
    "freq": np.array([1.050, 1.450, 1.640]), # reference frequency
    "p": np.array([0.086, 0.096, 0.099]), # fractional polarisation
    "a": np.array([33, 33, 33]), # polarisation angle in deg
}

The angles are conveniently all 33deg.

We can assume 3C286 has a RM of 0, so we can do a 'self-calibration' of the ionosphere by assuming any apparent RM is due to the ionosphere. We can obtain the RM using RM-Tools:

mDict_45, aDict_45 = run_rmsynth(
    [
        freqs, # frequencies in Hz
        stokes_i, # Stokes I flux
        stokes_q, # Stokes Q flux
        stokes_u, # Stokes U flux
        np.ones_like(stokes_i) * 1e-3, # dummy noise
        np.ones_like(stokes_q) * 1e-3, # dummy noise
        np.ones_like(stokes_u) * 1e-3, # dummy noise
    ],
    showPlots=False,
    phiMax_radm2=100 # Ionospheric RM is usually between 1-10rad/m/m
)
measured_intrinsic_angle = mDict_45["polAngle0Fit_deg"]

We should aim to be within a few degrees of 33deg, I think.

tjgalvin commented 7 months ago

Once we have done that correction then I am guessing we can also include some checks to make sure the Q and U is also correct against a known set of spectra, even if it is our own.

Once this basic test has been set up then we can expend upon this,