sc-depol

Introduction

This is not fully-fledged Python script. It is a simple shim availed to create a depolarization ratio map.

The depolarization ratio is the ratio of fractional polarisation at a shorter wavelength, to that at a longer wavelength. The fractional polarisation, \(P\), is defined using Stokes \(\mathit{I, Q, U}\) terms as:

\[ \begin{equation*} P = \frac{\sqrt{Q^{2} + U^{2}}}{\mathit{I}} \end{equation*} \]

Given fractional polarization at different frequencies, \(\nu_{max}\) and \(\nu_{min}\), the depolarization ratio is given as:

\[ \begin{equation*} \text{Depolarisation ratio}=\frac{P\text{ at }\nu_{max}}{P\text{ at }\nu_{min}} \end{equation*} \]

Docstring Documentation

Calculate the depolarization ratio

Parameters

(i|q|u)data: np.ndarray Data cubes containing the stokes e_file: str Name of the file containing errors. This is found in the file generated from scrappy's 'sc-los' tool and is a numpy pickle file. mask: str Name of mask to use while calculating depolarisation ratio output: str Name of output depolarisation ratio map

Returns

Maps of depolarization ratio and its error

Source code in src/rmsynthesis/rmmap_x2.py

def depolarisation_ratio(idata, qdata, udata, e_file, mask=None, output="depol-map"):
    """
    Calculate the depolarization ratio

    Parameters
    ----------
    (i|q|u)data: np.ndarray
        Data cubes containing the stokes
    e_file: str
        Name of the file containing errors. This is found in the file generated
        from scrappy's 'sc-los' tool and is a numpy pickle file.
    mask: str
        Name of mask to use while calculating depolarisation ratio
    output: str
        Name of output depolarisation ratio map        

    Returns
    -------
    Maps of depolarization ratio and its error

    """

    def depolzn_err(fp_lo, fer_lo, fp_hi, fer_hi):
        return np.sqrt(np.square(fer_lo/fp_lo) + np.square(fer_hi/fp_hi))

    def read_npz(fname):
        return dict(np.load(fname))

    ihdr = fits.getheader(idata)
    idata = fits.getdata(idata).squeeze()
    qdata = fits.getdata(qdata).squeeze()
    udata = fits.getdata(udata).squeeze()

    sel = [0,-1]

    npz = read_npz(e_file)
    ierr, qerr, uerr = npz["I_err"][sel], npz["Q_err"][sel], npz["U_err"][sel]

    idata, qdata, udata = idata[sel], qdata[sel], udata[sel]

    if mask:
        mask = fits.getdata(mask).squeeze()
        idata = np.ma.masked_equal(idata*mask, 0).filled()
        qdata = np.ma.masked_equal(qdata*mask, 0).filled()
        udata = np.ma.masked_equal(udata*mask, 0).filled()

    fpol = frac_polzn(idata, qdata, udata)



    fperr_lo = frac_polzn_error(idata[0], qdata[0], udata[0], ierr[0], qerr[0], uerr[0])
    fperr_hi = frac_polzn_error(idata[-1], qdata[-1], udata[-1], ierr[-1], qerr[-1], uerr[-1])


    fp_lo, fp_hi = fpol[0], fpol[-1]    

    depoln = fp_lo / fp_hi
    depoln_err = depolzn_err(fp_lo, fperr_lo, fp_hi, fperr_hi)
    depoln_err_ratio = depoln_err/depoln


    depol_hdr = modify_fits_header(ihdr, ctype='depol', unit='unit')

    #depolarization
    fits.writeto(output+".fits", depoln, depol_hdr,
        overwrite=True)

    #depolarization
    fits.writeto(output+"-err.fits", depoln_err, depol_hdr,
        overwrite=True)

    fits.writeto(output+"-err-ratio.fits", depoln_err_ratio, depol_hdr,
        overwrite=True)

    # write out fpol error at the highest freq (shortest wavelength)
    fperr_hdr = modify_fits_header(ihdr, ctype='fpol_err', unit='unit')
    fits.writeto(output+"-fpol-err.fits", fperr_hi, fperr_hdr,
        overwrite=True)


    snitch.info("Depolarization maps ready")
    return