Predict currently assumes local (as opposed to global) sky models

[Joshuaalbert]

@sjperkins I notice that according to this: https://sourceforge.net/p/wsclean/wiki/ComponentList/ we see that major/minor axes of ellipsoids are in arcseconds, however in the plane of sky the direction cosines are dimensionless so something must happen before predict makes sense. Do you know if a conversion is performed on the ellipsoidal parameters, and if so, can you point me to it?

[sjperkins]

This conversion from radec to lm depends on the phase centre, which we usually derive from the MeasurementSet. For e.g. in crystalball:

• https://github.com/caracal-pipeline/crystalball/blob/9235ed1f65a3739842cbd77b098f7f79879c4ad2/crystalball/crystalball.py#L218

The gaussian major/minor are converted to radians in the wsclean file loader:

codex-africanus/africanus/model/wsclean/file_model.py

Lines 73 to 74 in d7d7b6c

	"MajorAxis": arcsec2rad,
	"MinorAxis": arcsec2rad,

[Joshuaalbert]

Hi @sjperkins thanks for pointing it out. However, you realise from the equations that the major/minor axes can't actually be angular quantities right, since they define source distributions in the plane of the sky? For narrow fields of view it wouldn't matter but for all-sky or really wide fields it would matter where on the sky you place the source. Here's an example with Cas A with different phase tracking centres. See the colorbar scale changes, but more subtly the shape changes slightly.

Here is 1) directly pointing at Cas A

and 2) Pointing 20 degrees north of Cas A

I'll post below how I resolve this

Edit: axes labels are accidentally swapped.

[Joshuaalbert]

Instead what I do to resolve this issue is treat the major and minor axes rightfully as angular quantities and use spherical displacement to compute the major and minor axes on the plane of the sky. When you do this, you get almost no difference no matter where you place the source in the plane of the sky. See the colorbars don't move below.

Edit: axes labels are accidentally swapped.

[Joshuaalbert]

And to be clear the above is the gaussian components of the LOFAR Cas-a sky model from https://www.aanda.org/articles/aa/full_html/2020/03/aa36844-19/F2.html. In 1) I point directly at Cas A phase_tracking = ac.SkyCoord("-00h36m28.234s", "58d50m46.396s", frame='icrs'). In 2) I point 20 degrees norther phase_tracking = ac.SkyCoord("-00h36m28.234s", "78d50m46.396s", frame='icrs'). The time of observation is time = at.Time('2021-01-01T00:00:00', scale='utc') (in case you're taking into account aberration -- but that's not noticable here).

In the first plots (#303 (comment)) I use the method you use:

major_tangent = major.to(au.rad).value * au.dimensionless_unscaled
minor_tangent = minor.to(au.rad).value * au.dimensionless_unscaled

In the second I use a more nuanced approach:

s1_ra, s1_dec = offset_by(
    lon=source_directions.ra, lat=source_directions.dec,
    posang=theta, distance=major / 2.
)
s1 = ac.ICRS(s1_ra, s1_dec)
lmn1 = icrs_to_lmn(s1, time, phase_tracking)
major_tangent = 2. * np.linalg.norm(lmn1[:, :2] - lmn[:, :2], axis=-1) * au.dimensionless_unscaled
s2_ra, s2_dec = offset_by(
    lon=source_directions.ra, lat=source_directions.dec,
    posang=theta + au.Quantity(90, 'deg'), distance=minor / 2.
)
s2 = ac.ICRS(s2_ra, s2_dec)
lmn2 = icrs_to_lmn(s2, time, phase_tracking)
minor_tangent = 2. * np.linalg.norm(lmn2[:, :2] - lmn[:, :2], axis=-1) * au.dimensionless_unscaled

# position angle a bit more involved.

[landmanbester]

I think the assumption is that components coming out of wsclean are already in projected radians (for lack of a better unit, open to suggestions) since they just correspond to clean components defined on a tangent plane. If you want to simulate a wide field of view where the sources are actually defined on the celestial sphere you would have to take this projection effect into account. I suspect we have never really needed to worry about this but it would be a useful feature to support. Keep in mind that there is also an implicit assumption that the source is small enough so that we can approximate any DDE (including the w-term) to be constant across it

[Joshuaalbert]

Hi @landmanbester!

Yes, I know it's a clean component list, and typically these important standard calibrator sources are observed/imaged in the phase centre. For predicting for a a different phase centre should take into account this projection. However, I do not see that in any code. I have a simple correction for it that I use. I'm working on calibrating an all-sky LWA snapshot as a test bed for DSA calibration software. Going through a lot of other calibration software to understand conventions etc.

Re: constant w over source. I use a first-order approx I use to do slightly better:

F[I(l,m) * (C + (l - l0) * A + (m - m0) * B)] 
= F[I(l,m)] * (C - l0 * A - m0 * B) + A * i / (2pi) * d/du F[I(l,m)] + B * i / (2pi) * d/dv F[I(l,m)]

Using auto-diff it's easy to compute the terms. Works for any differentiable source description where we know the fourier transform.

[landmanbester]

For predicting for a a different phase centre should take into account this projection. However, I do not see that in any code. I have a simple correction for it that I use.

Indeed, we don't currently account for this in any of our software as far as I am aware. I think you can also predict to the reference phase center and rephase in visibility space but that does sound a bit cumbersome. I'm not sure if an ellipse on the sphere always stays an ellipse on the tangent plane though. Is that the case? If you can share your simple correction we could have a crack at putting that in.

Re: constant w over source. I use a first-order approx I use to do slightly better:

Yeah, I've had some vague thoughts along these same lines. It's probably time to give some serious thought to a jax implementation of the RIME which would make a lot of this much easier and also provide out of the box GPU support. If you can render your source to a pixelated grid on a tangent plane you can probably do pretty well with the wgridder if it's just the w-term you need to account for (have a look at the functions in ducc0.wgridder.experimental for versions that support arbitrary phase centers).

[Joshuaalbert]

I'm not sure if an ellipse on the sphere always stays an ellipse on the tangent plane though. Is that the case?

Nope, it does not. But you can approximate it as such.

If you can render your source to a pixelated grid on a tangent plane...

Yes, this does always seem to be the best way to generate model visibilities, however the clean component list is used for many applications, for it's compactness and lack of need for wgridding.

If you can share your simple correction...

It's pretty simple:

def get_constraint_points(posang, distance):
    s_ra, s_dec = offset_by(
        lon=source_directions.ra, lat=source_directions.dec,
        posang=posang, distance=distance
    )
    s = ac.ICRS(s_ra, s_dec)
    lmn = icrs_to_lmn(s, time, phase_tracking)
    return lmn[:, 0], lmn[:, 1]

# Offset by theta and a distance of half-major axis ==> half-major in tangent
l1, m1 = get_constraint_points(theta, major / 2.)
# Offset by theta + 90 and a distance of half-minor axis ==> half-minor in tangent
l2, m2 = get_constraint_points(theta + au.Quantity(90, 'deg'), minor / 2.)

# Approximation: Assume the remain orthogonal. They don't, but acceptable error.
major_tangent = 2. * np.sqrt((l1 - l0) ** 2 + (m1 - m0) ** 2)
minor_tangent = 2. * np.sqrt((l2 - l0) ** 2 + (m2 - m0) ** 2)

# Approximation: Assume orientations stays the same. It doesn't, but acceptable error.
theta_tangent = theta

This is the effect if the correction is not applied, when predicting from clean components:

[Joshuaalbert]

I think you can also predict to the reference phase center and rephase in visibility space but that does sound a bit cumbersome

For this you'd need to know the phase centre of the original observation from which clean component list was made wouldn't you?

[landmanbester]

For this you'd need to know the phase centre of the original observation from which clean component list was made wouldn't you?

Yeah you do. Isn't it stored in the component list?

[Joshuaalbert]

Nope

[o-smirnov]

Yeah it's the old distinction. In LOFAR days we used to call them global and local sky models.

A clean component list is an LSM by definition, and makes no sense without the original phase centre.

@Joshuaalbert is obviously in need of some GSM support, so let's discuss how best to implement that.

[landmanbester]

Yeah it's the old distinction. In LOFAR days we used to call them global and local sky models.

And are these a thing of the past ;-)?

Given that parametrizations are not preserved (eg. Gaussians do not remain Gaussian) under the projection from the sphere to a tangent plane this will have to be approximate. Unless we use a more appropriate coordinate system but then the analytic transforms will take a different form. @o-smirnov do you know of any other GSM predict implementations?