NIRCam coronagraphy coordinate transformation inconsistencies

[JarronL]

All of the NIRCam coronagraphic apertures appear to not be reversible when transforming between sci/det pixel coordinates and idl/tel angular coordinates, then back.

For instance:

ap = nrc_siaf['NRCA5_MASK335R']
xsci, ysci = (160., 160.)
# Convert to IDL coord frame
xidl, yidl = ap.sci_to_idl(xsci, ysci)
# Convert back to sci pixels (should be the same!)
xsci2, ysci2 = ap.idl_to_sci(xidl, yidl)
print(xsci2, ysci2)
# Returns (160.000490 159.549532), but expected (160.0, 160.0)

So, we're seeing an inconsistency of order -0.45 pixels along the y-axis. I tried different values across the entire aperture and get a similar result.

This appears to be primarily an issue with the coronagraphic apertures. I tried this on a couple direct imaging subarray (NRCA5_SUB320, NRCA5_SUB160P), and everything looked fine, with pixel consistency on the order 1E-4.

Here's a little script to test aperture difference across their field of view:

def test_ap_reverse(ap):
    xsci_arr = np.arange(1., ap.XSciSize+1)
    ysci_arr = np.arange(1., ap.YSciSize+1)
    xg, yg = np.meshgrid(xsci_arr, ysci_arr)
    xidl_arr, yidl_arr = ap.sci_to_idl(xg, yg)
    xg2, yg2 = ap.idl_to_sci(xidl_arr, yidl_arr)
    xdiff = xg - xg2
    ydiff = yg - yg2

    print(f"xdiff_mean: {np.mean(xdiff):.6f}, xdiff_min: {np.min(xdiff):.6f}, xdiff_max: {np.max(xdiff):.6f}")
    print(f"ydiff_mean: {np.mean(ydiff):.6f}, ydiff_min: {np.min(ydiff):.6f}, ydiff_max: {np.max(ydiff):.6f}")

And here are some results:

apname: NRCA5_TAMASK335R
xdiff_mean: -0.003125, xdiff_min: -0.003496, xdiff_max: -0.002652
ydiff_mean: 0.450492, ydiff_min: 0.450422, ydiff_max: 0.450539

apname: NRCA5_MASK335R
xdiff_mean: -0.000443, xdiff_min: -0.003339, xdiff_max: 0.001951
ydiff_mean: 0.450349, ydiff_min: 0.449878, ydiff_max: 0.450602

apname: NRCA5_FULL_MASK335R
xdiff_mean: -0.000106, xdiff_min: -0.013535, xdiff_max: 0.026609
ydiff_mean: 0.448994, ydiff_min: 0.444262, ydiff_max: 0.455327

apname: NRCA2_MASK335R
xdiff_mean: 0.003562, xdiff_min: 0.002166, xdiff_max: 0.004319
ydiff_mean: 0.700701, ydiff_min: 0.699110, ydiff_max: 0.702016

apname: NRCA2_FULL_MASK335R
xdiff_mean: 0.000496, xdiff_min: -0.014397, xdiff_max: 0.011275
ydiff_mean: 0.699327, ydiff_min: 0.690276, ydiff_max: 0.719769


# Direct imaging aperture
apname: NRCA5_FULL
xdiff_mean: -0.000180, xdiff_min: -0.042135, xdiff_max: 0.040512
ydiff_mean: 0.000038, ydiff_min: -0.001396, ydiff_max: 0.001195

apname: NRCA5_SUB160P
xdiff_mean: 0.000223, xdiff_min: -0.000471, xdiff_max: 0.000770
ydiff_mean: -0.000329, ydiff_min: -0.001342, ydiff_max: 0.000363

So, I suspect the transformation in one direction (pix->ang, or ang->pix) is "correct" but the other direction might no have been correctly updated during a revision to the distortion correction.

[mgennaro]

This has to do with the fact that we force the Sci2IdlY to NOT be 0 for coron apertures, to circumvent the problems with the COM discontinuity. I believe the corresponding CRDS distortion files correctly include that discontinuity explicitly as a step function (@bhilbert4 please confirm), but in SIAF we elected to use a trick that would make a polynomial solution work across the COM with focus towards the area behind it (more weight given to stars behind the COM when we determine the full-frame solution). I am not sure whether there is something with the inverse coeff, that is not fully accounting for that offset. I believe this doesn't impact operations, but indeed it seems to be impacting pysiaf.

[mgennaro]

To elaborate from yesterday's late night ramble: my suspicion is that the problem @JarronL has found is due to the fact that our transform now includes an offset of the origin and thus it is not linear but affine (i.e. for Xsci=0, Ysci=0 we get Xidl, Yidl != (0,0), but rather (Xidl, Yidl) == (0, Ycom_offset). Probably, somewhere in how @juliengirard derived the inverse coefficients, there is an assumption of linearity (and NOT affinity). So you have a direct transform that effectively goes Yidl = A*Ysci + B, where A is the linear part of the transform (from the polynomial coeff) and B is the Yoffset (the 0th term in the coeff). While I suspect the inverse was derived only for the A part, the equivalent of a "minus B" offset is missing in the inverse transform. See here: https://en.wikipedia.org/wiki/Affine_transformation. Do I know how to fix this? Probably not.

[bhilbert4]

The updated coron distortion files have not yet been delivered to CRDS. Last I heard, there was still some testing being done to make sure they work correctly. So it could be that @juliengirard has better/correct Sci2IdlY and inverse coefficients that just need to be delivered to pysiaf?

[JarronL]

@bhilbert4 It shouldn't really matter if the distortion files are old or new. It should still be reversible (to a relatively high precision) in the implementation.

[bhilbert4]

Of course. I was trying to say that we (probably) have new coefficients in hand that work correctly, in which case we can deliver these as new distortion files to CRDS, as well as to pysiaf (since pysiaf does not know anything about distortion files). At that point, things should round-trip correctly.

[mgennaro]

I suspect that new distortion coefficients will NOT fix this problem. My bet is that the problem is in the assumptions on how the inverse coefficients are calculated in the scripts that Armin/Julien tweaked ad-hoc for handling the COM. It's almost as if one should add the offset Jarron is seeing into the (0,0) coefficient for the Y direction in the Idl2Sci coefficients.