Ben's Polcal in Stimela/Cult-Cargo

[tmolteno]

A bit of a sprint at Africalim 2024. Requires the 0.2.0 branch which has the casa6 cabs.

#!/usr/bin/env -S stimela run
# Author. Tim tim@elec.ac.nz

opts:
    log:
        dir: logs/log-{config.run.datetime}
        name: log-{info.fqname}
        nest: 2
        symlink: log


_include:
  - (cultcargo):
    - wsclean.yml
    - quartical.yml
    - crystalball.yml
    - tricolour.yml
    - casa:
        - flag.yml
        - listobs.yml
        - calibration.yml
        - split.yml
        - setjy.yml
        - polcal.yml
        - bandpass.yml


moon-polcal:
    info: "A recipe for Lunar Polcal by Ben Hugo"

    inputs:
        ms:
            dtype: MS
            default: moon_uhf.ms
            info: "MeerKAT measurement set"
        cal-ms:
            dtype: MS
            default: moon_uhf_calibrators.ms
            info: "MeerKAT measurement set of calibrators"
        reference-ant:
            dtype: str
            default: 'm002'
            info: "reference antenna"
        image-fov:
            dtype: float
            default: 1.0
            info: "image field of view, in degrees"
        image-pixel-res:
            dtype: float
            default: 3
            info: "pixel size, in arcsec"
        fluxcal:
            dtype: str
            default: 'J1939-6342'
        phasecal:
            dtype: str
            default: J1733-1304
        polcal:
            dtype: str
            default: J1733-1304

    assign:
        image-size: = 2*recipe.image-fov * 3600.0 / recipe.image-pixel-res
        dirs:
            temp: "{config.run.env.HOME}/tmp"    # temp files go here
            base: .                              # base project directory -- directory of recipe by default, but see below
        dir-out: '{recipe.dirs.base}/out'        # output products go here
        log.dir: '{recipe.dir-out}/logs/log-{config.run.datetime}'   # put logs into output dir

    steps:
        listobs-moon:
            info: Look at the moon fields
            cab: casa.listobs
            params:
                ms: =recipe.ms

        listobs-cal:
            info: Look at the moon fields
            cab: casa.listobs
            params:
                ms: =recipe.cal-ms

        #
        # We will set and transfer the flux scale and the spectral (bandpass response) of the system using J1939 here
        # We solve for the amplitudes (mainly) and the leakages (non-orthogonality of the feeds) using this unpolarized source
        #
        # We remove the delay (K Jones - frequency slope), and average gain (constant offset) before doing the (normalized) bandpass.
        #
        # To simplify things for now we will use a single component spectral model for this unresolved source
        # setjy(vis='moon_uhf_calibrators.ms', standard='Stevens-         
        #     Reynolds 2016', field='J1939-6342',usescratch=True)
        #
        flux-density-cal:
            info: Absolute flux density calibration
            cab: casa.setjy
            params:
                ms: =recipe.cal-ms
                standard: 'Stevens-Reynolds 2016'
                field: 'J1939-6342'
                usescratch: true

        # gaincal(vis='moon_uhf_calibrators.ms', caltable='bp.K', field='J1939-6342', gaintype='K', solint='inf', combine='',refant='m002')
        #

        #
        # Delay and Bandpass Calibration
        #
        delay-cal:
            info: Delay Calibration
            cab: casa.gaincal
            params:
                ms: =recipe.cal-ms
                caltable: 'bp.K'
                field: =recipe.fluxcal
                refant: =recipe.reference-ant
                solnorm: true
                gaintype: 'K'
                combine: ''
                solint: 'inf'

        # gaincal(vis='moon_uhf_calibrators.ms', caltable='bp.Gp', field='J1939-6342', gaintype='G', calmode='p', solint='inf', combine='',refant='m002', gaintable=['bp.K'])

        gain-cal-p:
            info: Initial Gain Calibration
            cab: casa.gaincal
            params:
                ms: =recipe.cal-ms
                caltable: 'bp.Gp'
                field: =recipe.fluxcal
                solint: 'inf'
                refant: =recipe.reference-ant
                minsnr: 2.0
                gaintype: 'G'
                gaintable: [bp.K]

        # gaincal(vis='moon_uhf_calibrators.ms', caltable='bp.Ga', field='J1939-6342', gaintype='G', calmode='a', solint='inf', combine='',refant='m002', gaintable=['bp.K', 'bp.Gp'])

        gain-cal-a:
            info: Initial Gain Calibration
            cab: casa.gaincal
            params:
                ms: =recipe.cal-ms
                caltable: 'bp.Ga'
                field: =recipe.fluxcal
                solint: 'inf'
                refant: =recipe.reference-ant
                minsnr: 2.0
                calmode: a
                gaintype: 'G'
                gaintable: [bp.K, bp.Gp]

        
        # bandpass(vis='moon_uhf_calibrators.ms', caltable='bp.B', field='J1939-6342', solint='inf', combine='scan',refant='m002', gaintable=['bp.K', 'bp.Gp', 'bp.Ga'])

        bandpass-cal:
            info: Bandpass Calibration
            cab: casa.bandpass
            params:
                ms: =recipe.cal-ms
                caltable: 'bp.B'
                field: =recipe.fluxcal
                solint: 'inf'
                combine: 'scan'
                refant: =recipe.reference-ant
                gaintable: [bp.K, bp.Gp, bp.Ga]

        #
        # Transfer of gain calibration to the target
        #
        # applycal(vis='moon_uhf_calibrators.ms', field='J1939-6342', gaintable=['bp.K','bp.Gp','bp.Ga','bp.B']) 

        transfer-gain-amplitude:
            info: Apply calibration to the amplitude and phase calibrators
            cab: casa.applycal
            params:
                ms: =recipe.cal-ms
                field: =recipe.fluxcal
                gaintable: [bp.K, bp.Gp, bp.Ga, bp.B]
                
        # polcal(vis='moon_uhf_calibrators.ms', field='J1939-6342', caltable='bp.Df', refant='m002', gaintable=['bp.K','bp.Gp','bp.Ga','bp.B'], combine='', poltype='Df', solint='inf')

        polcal:
            info: Do the polarization cal
            cab: casa.polcal
            params:
                ms: =recipe.cal-ms
                field: =recipe.fluxcal
                poltype: 'Df'
                caltable: bp.Df
                gaintable: [bp.K, bp.Gp, bp.Ga, bp.B]
                refant: =recipe.reference-ant
                solint: 'inf'
        
        # applycal(vis='moon_uhf_calibrators.ms', field='J1939-6342', gaintable=['bp.K','bp.Gp','bp.Ga','bp.B','bp.Df'])

        apply-polcal:
            info: Apply calibration to the amplitude and phase calibrators
            cab: casa.applycal
            params:
                ms: =recipe.cal-ms
                field: =recipe.fluxcal
                gaintable: [bp.K, bp.Gp, bp.Ga, bp.B, bp.Df]

        # # bootstrap the secondary and make some plots!
        # applycal(vis='moon_uhf_calibrators.ms', field='J1733-1304', gaintable=['bp.K','bp.Gp','bp.Ga','bp.B','bp.Df'])

        bootstrap-secondary:
            info: Apply calibration to the target
            cab: casa.applycal
            params:
                ms: =recipe.cal-ms
                field: =recipe.phasecal
                gaintable: [bp.K, bp.Gp, bp.Ga, bp.B, bp.Df]

        # gaincal(vis='moon_uhf_calibrators.ms', caltable='sec.Gp', field='J1733-1304', gaintype='G', calmode='p', solint='inf', combine='',refant='m002', gaintable=['bp.K', 'bp.Ga', 'bp.B', 'bp.Df'])
        secondary-gaincal:
            info: Bandpass Calibration
            cab: casa.gaincal
            params:
                ms: =recipe.cal-ms
                caltable: 'sec.Gp'
                field: =recipe.phasecal
                gaintype: 'G'
                calmode: 'p'
                solint: 'inf'
                refant: =recipe.reference-ant
                gaintable: [bp.K, bp.Ga, bp.B, bp.Df]

        # # please note solnorm=True is set here because we do not have a model for the secondary. In a full workflow you should derive a bootstrapped model for the secondary before running an amplitude gain!!!
        # gaincal(vis='moon_uhf_calibrators.ms', caltable='sec.T', field='J1733-1304', gaintype='T', calmode='ap', solnorm=True, solint='inf', combine='',refant='m002', gaintable=['bp.K', 'sec.Gp', 'bp.Ga', 'bp.B', 'bp.Df'])

        secondary-tcal:
            info: Bandpass Calibration
            cab: casa.gaincal
            params:
                ms: =recipe.cal-ms
                caltable: 'sec.T'
                field: =recipe.phasecal
                gaintype: 'T'
                calmode: 'ap'
                solnorm: true
                solint: 'inf'
                refant: =recipe.reference-ant
                gaintable: [bp.K, sec.Gp, bp.Ga, bp.B, bp.Df]

                
        # run applycal again and check the phase v amp now with the command above
        # applycal(vis='moon_uhf_calibrators.ms', field='J1733-1304', gaintable=['bp.K','sec.Gp','bp.Ga','bp.B','sec.T','bp.Df'])
        reapply-cal:
            info: Apply calibration to the target
            cab: casa.applycal
            params:
                ms: =recipe.cal-ms
                field: =recipe.phasecal
                gaintable: [bp.K, sec.Gp, bp.Ga, bp.B, sec.T, bp.Df]

        # let's apply these phase refined solutions and frequency-dependent amplitude solutions from our primary onto our target, the Moon and see what things look like!
        # applycal(vis='moon_uhf.ms', field='Moon', gaintable=['bp.K','sec.Gp','bp.Ga','bp.B','sec.T','bp.Df'])

        refined-phase-apply:
            info: Apply calibration to the target
            cab: casa.applycal
            params:
                ms: =recipe.ms
                field: 'Moon'
                gaintable: [bp.K, sec.Gp, bp.Ga, bp.B, sec.T, bp.Df]

        first-image:
            info: "make a dirty image"
            cab: wsclean
            params:
                ms: =recipe.ms
                prefix: 'first-image'
                size: =recipe.image-size
                column: CORRECTED_DATA
                auto-threshold: 3
                scale: '{recipe.image-pixel-res}arcsec'
                mgain: 0.9
                niter: 10000
                pol: "IQUV"
                temp-dir: =recipe.dirs.temp

                
        # applycal(vis='moon_uhf_calibrators.ms', field='J1331+3030', gaintable=['bp.K','bp.Gp','bp.Ga','bp.B','sec.T','bp.Df'])

        refined-phase-reapply:
            info: Apply calibration to the target
            cab: casa.applycal
            params:
                ms: =recipe.cal-ms
                field: =recipe.polcal
                gaintable: [bp.K, bp.Gp, bp.Ga, bp.B, sec.T, bp.Df]

        # let's refine the phase -- obviously this source is quite polarized (10% at L-band dropping to <2% at the bottom of UHF at 600MHz).
        # MAKE SURE YOU CALIBRATE ONLY THE PHASE - otherwise you will absorb the source's unmodelled polarization into your transfer 
        # solutions and bias your entire target field as a result!
        # gaincal(vis='moon_uhf_calibrators.ms', caltable='pol.Gp', field='J1331+3030', gaintype='G', calmode='p', solint='inf', combine='',refant='m002', gaintable=['bp.K', 'bp.Ga', 'bp.B', 'bp.Df'])
        
        refine-phase-2:
            info: Phase Calibration
            cab: casa.gaincal
            params:
                ms: =recipe.cal-ms
                caltable: 'pol.Gp'
                field: =recipe.polcal
                gaintype: 'G'
                calmode: 'p'
                solint: 'inf'
                refant: =recipe.reference-ant
                gaintable: [bp.K, bp.Ga, bp.B, bp.Df]

        # applycal(vis='moon_uhf_calibrators.ms', field='J1331+3030', gaintable=['bp.K','pol.Gp','bp.Ga','bp.B','sec.T','bp.Df'])
        refined-phase-reapply-:
            info: Apply calibration to the target
            cab: casa.applycal
            params:
                ms: =recipe.cal-ms
                field: =recipe.polcal
                gaintable: [bp.K, pol.Gp, bp.Ga, bp.B, sec.T, bp.Df]

        # There is fractional V Stokes being picked up here -- the goal is to remove that
        # underneath the hood we take arctan2(V/U) here, assuming V~=0 Jy
        # we can also solve for a crosshand delay (mode KCROSS in gaincal). This is necessary at lower ranges in UHF especially due to SNR
        # concerns but this frequency dependent slope is stable at few ps level so we will just absorb it into crosshand phase (Xf)
        # solutions here

        # polcal(vis='moon_uhf_calibrators.ms', field='J1331+3030', caltable='pol.Xf', refant='m002', gaintable=['bp.K','pol.Gp','bp.Ga','bp.B'], combine='scan', poltype='Xf', solint='inf,20MHz')
        # TODO solint causes error in casa6
        polcal-2:
            info: Do the polarization cal
            cab: casa.polcal
            params:
                ms: =recipe.cal-ms
                field: =recipe.polcal
                caltable: pol.Xf
                gaintable: [bp.K, pol.Gp, bp.Ga, bp.B]
                refant: =recipe.reference-ant
                poltype: Xf
                solint: 'inf'
                combine: scan
                

        # apply and inspect the net result, time variability should be mostly gone (modulo ionosphere) and real dominated as discussed
        # applycal(vis='moon_uhf_calibrators.ms', field='J1331+3030', gaintable=['bp.K','pol.Gp','bp.Ga','bp.B','sec.T','bp.Df','pol.Xf'], parang=True)
                
        apply-polcal-2:
            info: Apply calibration to the amplitude and phase calibrators
            cab: casa.applycal
            params:
                ms: =recipe.cal-ms
                field: =recipe.polcal
                gaintable: [bp.K, pol.Gp, bp.Ga, bp.B, sec.T, bp.Df, pol.Xf]
                parang: true

        # all that is left is transfer of the full chain onto the Moon and verifying that the V goes down on the Lunar maps!
        # applycal(vis='moon_uhf.ms', field='Moon', gaintable=['bp.K','sec.Gp','bp.Ga','bp.B','sec.T','bp.Df','pol.Xf'], parang=True)

        transfer-polcal-2:
            info: Apply calibration to the amplitude and phase calibrators
            cab: casa.applycal
            params:
                ms: =recipe.ms
                field: 'Moon'
                gaintable: [bp.K, sec.Gp, bp.Ga, bp.B, sec.T, bp.Df, pol.Xf]
                parang: true

        last-image:
            info: "make a final image"
            cab: wsclean
            params:
                ms: =recipe.ms
                prefix: 'last-image'
                size: =recipe.image-size
                column: CORRECTED_DATA
                auto-threshold: 3
                scale: '{recipe.image-pixel-res}arcsec'
                mgain: 0.9
                niter: 10000
                pol: "IQUV"
                temp-dir: =recipe.dirs.temp