Compare calc11 with DSA2000

[Joshuaalbert]

@aelanman how about I provide some obs parameters and you produ e your delays and I produce mine then we compare.

[aelanman]

Sounds good!

[aelanman]

@Joshuaalbert which parameters would you like to try?

[Joshuaalbert]

Here you go. You should only need astropy installed to create the standard setup.

import itertools
import numpy as np
from astropy import time as at, units as au, coordinates as ac
from astropy.coordinates import offset_by

def prepare_standard_test():
    np.random.seed(42)

    obstime = at.Time("2024-01-01T00:00:00", scale='utc')
    array_location = ac.EarthLocation.of_site('vla')

    phase_centers = []

    for body in ['sun', 'jupiter', 'moon', 'neptune', 'mars']:
        body_position_bcrs, body_velocity_bcrs = ac.get_body_barycentric_posvel(
            body=body,
            time=obstime
        )  # [T, N]
        frame = ac.CartesianRepresentation(body_position_bcrs)
        source = ac.ICRS().realize_frame(frame)

        for displacement_deg in [0., 0.25, 0.5, 1., 2., 4., 8.]:
            # displace by 1 deg in some random direction
            new_ra, new_dec = offset_by(
                lon=source.ra,
                lat=source.dec,
                posang=np.random.uniform(0., 2 * np.pi) * au.rad,
                distance=displacement_deg * au.deg
            )
            phase_centers.append(ac.ICRS(ra=new_ra, dec=new_dec))

    antennas = []

    for b_east in [0, 3, 10, 100]:
        for b_north in [0, 3, 10, 100]:
            az = np.arctan2(b_east, b_north) * au.rad
            dist = np.sqrt(b_east ** 2 + b_north ** 2)
            antenna = ac.AltAz(
                az=az,
                alt=0 * au.deg,
                distance=dist * au.km,
                location=array_location,
                obstime=obstime
            )
            antennas.append(antenna.transform_to(ac.ITRS(obstime=obstime, location=array_location)))

    antennas = ac.concatenate(antennas).earth_location

    antenna_1, antenna_2 = np.asarray(list(itertools.combinations(range(len(antennas)), 2))).T

    return obstime, array_location, antennas, antenna_1, antenna_2, phase_centers

For each phase center produce the UVW coordinates ordered by the antenna1 and antenna2, and save to file. Here's mine.

def test_standard_test_dsa2000():
    obstime, array_location, antennas, antenna_1, antenna_2, phase_centers = prepare_standard_test()

    results = {}

    for i, phase_center in enumerate(phase_centers):
        far_field_engine = FarFieldDelayEngine(
            antennas=antennas,
            phase_center=phase_center,
            start_time=obstime,
            end_time=obstime,
            verbose=True
        )
        times = jnp.repeat(far_field_engine.time_to_jnp(obstime), len(antenna_1))
        uvw = far_field_engine.compute_uvw_jax(
            times=times,
            antenna_1=antenna_1,
            antenna_2=antenna_2
        )  # [N, 3] in meters
        results[str(i)] = uvw.tolist()

    with open('dsa2000_results.json', 'w') as f:
        json.dump(results, f, indent=2)

[Joshuaalbert]

To be sure, I'm using

astropy==6.1.2
astropy-iers-data==0.2024.8.5.0.32.23

[aelanman]

Hi @Joshuaalbert -- Thanks for the code. Here's what I have so far (note, I have astropy==6.1.0 and astropy-iers-data==0.2024.5.6.0.29.28).

Here's a snippet of the code I ran. The json file of result is attached. I noticed that I get significantly different ITRS antenna positions from pyuvdata's ECEF_from_ENU function than I did from your code (transforming AltAz to ITRS), but it's irrelevant to the current comparison.

CALC returns UVWs for baselines formed between the antennas and the geocenter, which I convert to baseline UVWs by subtracting.

import json
from pycalc11.interface import Calc

obstime, array_location, antennas, antenna_1, antenna_2, phase_centers = prepare_standard_test()
station_coords = antennas

srcs = ac.SkyCoord(ra=1.5, dec=-30.4, unit='deg', frame='icrs')
ci = Calc(
    station_names=[f"VLA_{ii}" for ii in range(16)],        # Each antenna needs a unique name
    station_coords=station_coords,
    start_time=obstime - at.TimeDelta(30, format='sec'),
    duration_min=5.,
    source_coords=ac.concatenate(phase_centers),
    check_sites=False,
    dry_atm=False,
    wet_atm=False,
)

ci.run_driver()

uvws = ci.uvw_interp(obstime)

# Form baseline UVWs from antenna_1 and 2
bl_uvws = uvws[0, antenna_1, :, :] - uvws[0, antenna_2, :, :]

# Write to file
results = {str(si) : bl_uvws[:, si, :].to_value('m').tolist() for si in range(len(phase_centers))}
with open('pycalc_results.json', 'w') as f:
    json.dump(results, f, indent=2)

pycalc_results.json

Let's see if we're at least in the same ballpark and then track down finer differences as needed.

[Joshuaalbert]

I'm getting some erratic differences, that suggest there may have been a problem with your calculation of antenna coordinates.

def test_compare_to_calc():
    with open('dsa2000_results.json', 'r') as f:
        results_dsa = json.load(f)
    with open('pycalc_results.json', 'r') as f:
        results_calc = json.load(f)

    obstime, array_location, antennas, antenna_1, antenna_2, phase_centers = prepare_standard_test()

    antennas_gcrs = antennas.get_gcrs(obstime)
    antennas_gcrs = antennas_gcrs.cartesian.xyz.T # [N, 3]

    baselines = antennas_gcrs[antenna_2, :] - antennas_gcrs[antenna_1, :]
    baselines_norm = np.linalg.norm(baselines, axis=1).to('m').value

    for key, phase_center in enumerate(phase_centers):
        alt = phase_center.transform_to(ac.AltAz(obstime=obstime, location=array_location)).alt.to('deg')

        value_calc = np.asarray(results_calc[str(key)])
        value_dsa = np.asarray(results_dsa[str(key)])
        value_dsa_norm = np.linalg.norm(value_dsa, axis=1)
        value_calc_norm = np.linalg.norm(value_calc, axis=1)
        for i1, i2, b_norm, dsa_norm, calc_norm in zip(antenna_1, antenna_2, baselines_norm, value_dsa_norm, value_calc_norm):
            if np.abs(dsa_norm - calc_norm) > 10:
                print(f"alt: {alt}, phase center: {phase_center.ra} {phase_center.dec} i1: {i1} i2: {i2} Baseline: {b_norm}, DSA: {dsa_norm}, Calc: {calc_norm}")

        diff = (value_dsa - value_calc)
        mean_diff = np.mean(diff)
        std_diff = np.std(diff)
        print(f"alt: {alt}, phase center: {phase_center.ra} {phase_center.dec} Mean diff: {mean_diff}, Std diff: {std_diff}")

Firstly, there are numerous situations where calc11 seems to have baseline norms much larger than expected.

alt: -1.7089649902803516 deg, phase center: 101.23523170067013 deg 23.05625429855563 deg i1: 0 i2: 1 Baseline: 3000.0000000001514, DSA: 2999.995567159286, Calc: 5130.00832333731
alt: -1.7089649902803516 deg, phase center: 101.23523170067013 deg 23.05625429855563 deg i1: 0 i2: 2 Baseline: 9999.999999999716, DSA: 9999.985230146925, Calc: 19215.574769375286
alt: -1.7089649902803516 deg, phase center: 101.23523170067013 deg 23.05625429855563 deg i1: 0 i2: 3 Baseline: 100000.00000000015, DSA: 99999.85310914909, Calc: 175525.9270168726
alt: -1.7089649902803516 deg, phase center: 101.23523170067013 deg 23.05625429855563 deg i1: 0 i2: 4 Baseline: 3000.0000000001137, DSA: 2999.995586034651, Calc: 8249.578161041158
...

And, then with regard to mean diff and std of diff some are ~cm's but others are much larger. Definitely, the error from above would affect this.

alt: -16.000295314055546 deg, phase center: 258.14092536218783 deg -23.32562920151148 deg Mean diff: 1.0024403872218832, Std diff: 1.5010545936659367
alt: -15.385621889939655 deg, phase center: 258.04189759660414 deg -21.969838967355134 deg Mean diff: 1.0933058952702472, Std diff: 1.6264004712524232
alt: -16.140209723664153 deg, phase center: 255.84210751347544 deg -20.006281256701556 deg Mean diff: 0.989912851472601, Std diff: 1.4516475312872446
alt: -14.402022982582498 deg, phase center: 255.40734734738976 deg -15.983604017405893 deg Mean diff: 1.2726459170162336, Std diff: 1.828275605965074
...

Here's mine btw: dsa2000_results.json

[Joshuaalbert]

I noticed that I get significantly different ITRS antenna positions from pyuvdata's ECEF_from_ENU function than I did from your code (transforming AltAz to ITRS), but it's irrelevant to the current comparison.

Not sure this is the same thing but to convert AltAz to ITRS you need topocentric ITRS. See this ticket where I encountered this. https://github.com/astropy/astropy/issues/16317#issuecomment-2068394464

[aelanman]

Is it just those four cases that are far off? It's possible that something's going wrong with pycalc's handling of atmospheric corrections in those cases.

Here's a json file with the delays on each baseline in μs: calc_delays.json

[aelanman]

For another comparison, here are the uvw and delay results with the atmospheric corrections disabled: calc_delays_noatmo.json pycalc_results_noatmo.json

[Joshuaalbert]

Indeed, with atmosphere turned off the big errors disappear, and we get:

alt: -53.03969179667606 deg, phase center: 199.0411607343037 deg -6.5183238618360235 deg Mean diff: 0.005133341332866969, Std diff: 0.053784907357853226
alt: -52.90016882015525 deg, phase center: 198.9645146243513 deg -6.280209687589366 deg Mean diff: 0.005250530950046044, Std diff: 0.053806254741862077
alt: -53.39143050850489 deg, phase center: 198.54106850487892 deg -6.574524132727084 deg Mean diff: 0.004777458674998493, Std diff: 0.053784950017652974
alt: -54.038236401747284 deg, phase center: 198.45544576295353 deg -7.331915438882114 deg Mean diff: 0.004190159219863669, Std diff: 0.05372104129738027
alt: -51.13048039018375 deg, phase center: 200.70979183749034 deg -5.402176859153181 deg Mean diff: 0.006931350226070136, Std diff: 0.053810741667083425
alt: -49.21373249754742 deg, phase center: 202.37193302071898 deg -4.280990902474353 deg Mean diff: 0.008649890041554266, Std diff: 0.05375199617336047
alt: -45.3573997373658 deg, phase center: 201.88868529881321 deg 0.959793548491772 deg Mean diff: 0.011754301761207788, Std diff: 0.054067486412539424
alt: 40.5645441456519 deg, phase center: 43.40583670184876 deg 15.475912094428624 deg Mean diff: 0.00865157138112473, Std diff: 0.052799119946482564
alt: 40.58450242197394 deg, phase center: 43.252044496367446 deg 15.274639376046006 deg Mean diff: 0.008528766954864844, Std diff: 0.05286841303464343
alt: 40.90578168182581 deg, phase center: 42.90523975253777 deg 15.345150286680525 deg Mean diff: 0.008451425014748064, Std diff: 0.05302867218075489
alt: 40.96992052801057 deg, phase center: 43.54032271664111 deg 16.467518590316228 deg Mean diff: 0.00908435751420347, Std diff: 0.05276341120586875
alt: 41.90037465775388 deg, phase center: 43.01192742221875 deg 17.43991501042299 deg Mean diff: 0.009335328971397376, Std diff: 0.053034343686807
alt: 44.56218127540172 deg, phase center: 39.76393655154233 deg 17.42580498937576 deg Mean diff: 0.008336316025678059, Std diff: 0.0544394381265021
alt: 34.720080072077224 deg, phase center: 51.547863703913094 deg 17.19829327493339 deg Mean diff: 0.010897006954859378, Std diff: 0.048208568221101575
alt: -1.7089649902803516 deg, phase center: 101.23523170067013 deg 23.05625429855563 deg Mean diff: 0.0009366121501465664, Std diff: 0.031512465789100316
alt: -1.825243387562478 deg, phase center: 101.48368206747512 deg 23.157642076553902 deg Mean diff: 0.0008813139517504472, Std diff: 0.03162269084168736
alt: -2.182344546732068 deg, phase center: 101.74674772653613 deg 22.888305144192973 deg Mean diff: 0.000729795869561182, Std diff: 0.03192577418735819
alt: -2.1995347136652814 deg, phase center: 101.06806983295284 deg 22.068240397238263 deg Mean diff: 0.0007414073047633628, Std diff: 0.031880370928632495
alt: -3.487570757482375 deg, phase center: 102.12487448841954 deg 21.233869808952733 deg Mean diff: 0.00017832551176060935, Std diff: 0.032977548698183826
alt: -5.377100968603685 deg, phase center: 105.40505017151173 deg 21.977290336136996 deg Mean diff: -0.000647499025244416, Std diff: 0.034917165391090965
alt: -1.4193298212816532 deg, phase center: 95.84162530439148 deg 16.864626712746325 deg Mean diff: 0.0013413789798532348, Std diff: 0.030922380921618833
alt: 53.225191856549095 deg, phase center: 357.33881848021724 deg -2.512062003289685 deg Mean diff: -0.027679094084581975, Std diff: 0.05131619996153931
alt: 53.12992053416528 deg, phase center: 357.58034882135206 deg -2.577469216659949 deg Mean diff: -0.027547145970373752, Std diff: 0.051384743483953145
alt: 52.84677044360912 deg, phase center: 357.7114884986282 deg -2.8458596305599086 deg Mean diff: -0.027554530201485614, Std diff: 0.05146218834386254
alt: 52.23645153364766 deg, phase center: 357.61183630740857 deg -3.4741775152858105 deg Mean diff: -0.02782200687661093, Std diff: 0.05156009110264898
alt: 53.85846091689707 deg, phase center: 355.38621011276956 deg -2.072245354859094 deg Mean diff: -0.02868860271310538, Std diff: 0.05077028950041659
alt: 53.808195485694355 deg, phase center: 1.1411298385724182 deg -1.2635840686725994 deg Mean diff: -0.024722038614073977, Std diff: 0.05208547963589685
alt: 45.377620916048 deg, phase center: 356.6144894001537 deg -10.47967957817007 deg Mean diff: -0.030004875497425806, Std diff: 0.05289992383723178
alt: -16.999865394016005 deg, phase center: 257.1849045823964 deg -23.80743220393574 deg Mean diff: -0.003319799903335034, Std diff: 0.10915266755172588
alt: -16.817112437033597 deg, phase center: 257.26338326445006 deg -23.56798454591868 deg Mean diff: -0.003499723437833069, Std diff: 0.10988942743295312
alt: -17.466368571748365 deg, phase center: 256.84205252036094 deg -24.19717965052778 deg Mean diff: -0.0026703721158995108, Std diff: 0.10625577949146608
alt: -16.000295314055546 deg, phase center: 258.14092536218783 deg -23.32562920151148 deg Mean diff: -0.005081875688521128, Std diff: 0.117494907406081
alt: -15.385621889939655 deg, phase center: 258.04189759660414 deg -21.969838967355134 deg Mean diff: -0.004945661704209236, Std diff: 0.11654237682834045
alt: -16.140209723664153 deg, phase center: 255.84210751347544 deg -20.006281256701556 deg Mean diff: -0.0009574952564445145, Std diff: 0.0984002836613996
alt: -14.402022982582498 deg, phase center: 255.40734734738976 deg -15.983604017405893 deg Mean diff: 0.0014624907625040528, Std diff: 0.09043430430779925

[Joshuaalbert]

[aelanman]

Nice! Is this just all Us, Vs, and Ws together? So we're getting ~ 10 cm errors for some coordinates but most are consistent to within cm?

[Joshuaalbert]

Yes the flattened diff. It looks okay. I also don't have any atmospheric effects.

[aelanman]

@Joshuaalbert What do you want to do next?

[Joshuaalbert]

Maybe a call to discuss?On Sept 13, 2024 16:26, Adam Lanman @.***> wrote: @Joshuaalbert What do you want to do next?

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: @.***>

[aelanman]

@Joshuaalbert I sent an email to what I think is the correct email for you (at leidenuniv.nl), but let me know if you have a different preferred address.

[Joshuaalbert]

Sent an invite to your email

[Joshuaalbert]

With DSA using DE430 and pycalc using DE421:

alt: -53.03969179667606 deg, phase center: 199.0411607343037 deg -6.5183238618360235 deg Mean diff: 0.005133341332866969, Std diff: 0.053784907357853226
alt: -52.90016882015525 deg, phase center: 198.9645146243513 deg -6.280209687589366 deg Mean diff: 0.005250530950046044, Std diff: 0.053806254741862077
alt: -53.39143050850489 deg, phase center: 198.54106850487892 deg -6.574524132727084 deg Mean diff: 0.004777458674998493, Std diff: 0.053784950017652974
alt: -54.038236401747284 deg, phase center: 198.45544576295353 deg -7.331915438882114 deg Mean diff: 0.004190159219863669, Std diff: 0.05372104129738027
alt: -51.13048039018375 deg, phase center: 200.70979183749034 deg -5.402176859153181 deg Mean diff: 0.006931350226070136, Std diff: 0.053810741667083425
alt: -49.21373249754742 deg, phase center: 202.37193302071898 deg -4.280990902474353 deg Mean diff: 0.008649890041554266, Std diff: 0.05375199617336047
alt: -45.3573997373658 deg, phase center: 201.88868529881321 deg 0.959793548491772 deg Mean diff: 0.011754301761207788, Std diff: 0.054067486412539424
alt: 40.5645441456519 deg, phase center: 43.40583670184876 deg 15.475912094428624 deg Mean diff: 0.00865157138112473, Std diff: 0.052799119946482564
alt: 40.58450242197394 deg, phase center: 43.252044496367446 deg 15.274639376046006 deg Mean diff: 0.008528766954864844, Std diff: 0.05286841303464343
alt: 40.90578168182581 deg, phase center: 42.90523975253777 deg 15.345150286680525 deg Mean diff: 0.008451425014748064, Std diff: 0.05302867218075489
alt: 40.96992052801057 deg, phase center: 43.54032271664111 deg 16.467518590316228 deg Mean diff: 0.00908435751420347, Std diff: 0.05276341120586875
alt: 41.90037465775388 deg, phase center: 43.01192742221875 deg 17.43991501042299 deg Mean diff: 0.009335328971397376, Std diff: 0.053034343686807
alt: 44.56218127540172 deg, phase center: 39.76393655154233 deg 17.42580498937576 deg Mean diff: 0.008336316025678059, Std diff: 0.0544394381265021
alt: 34.720080072077224 deg, phase center: 51.547863703913094 deg 17.19829327493339 deg Mean diff: 0.010897006954859378, Std diff: 0.048208568221101575
alt: -1.7089649902803516 deg, phase center: 101.23523170067013 deg 23.05625429855563 deg Mean diff: 0.0009366121501465664, Std diff: 0.031512465789100316
alt: -1.825243387562478 deg, phase center: 101.48368206747512 deg 23.157642076553902 deg Mean diff: 0.0008813139517504472, Std diff: 0.03162269084168736
alt: -2.182344546732068 deg, phase center: 101.74674772653613 deg 22.888305144192973 deg Mean diff: 0.000729795869561182, Std diff: 0.03192577418735819
alt: -2.1995347136652814 deg, phase center: 101.06806983295284 deg 22.068240397238263 deg Mean diff: 0.0007414073047633628, Std diff: 0.031880370928632495
alt: -3.487570757482375 deg, phase center: 102.12487448841954 deg 21.233869808952733 deg Mean diff: 0.00017832551176060935, Std diff: 0.032977548698183826
alt: -5.377100968603685 deg, phase center: 105.40505017151173 deg 21.977290336136996 deg Mean diff: -0.000647499025244416, Std diff: 0.034917165391090965
alt: -1.4193298212816532 deg, phase center: 95.84162530439148 deg 16.864626712746325 deg Mean diff: 0.0013413789798532348, Std diff: 0.030922380921618833
alt: 53.225191856549095 deg, phase center: 357.33881848021724 deg -2.512062003289685 deg Mean diff: -0.027679094084581975, Std diff: 0.05131619996153931
alt: 53.12992053416528 deg, phase center: 357.58034882135206 deg -2.577469216659949 deg Mean diff: -0.027547145970373752, Std diff: 0.051384743483953145
alt: 52.84677044360912 deg, phase center: 357.7114884986282 deg -2.8458596305599086 deg Mean diff: -0.027554530201485614, Std diff: 0.05146218834386254
alt: 52.23645153364766 deg, phase center: 357.61183630740857 deg -3.4741775152858105 deg Mean diff: -0.02782200687661093, Std diff: 0.05156009110264898
alt: 53.85846091689707 deg, phase center: 355.38621011276956 deg -2.072245354859094 deg Mean diff: -0.02868860271310538, Std diff: 0.05077028950041659
alt: 53.808195485694355 deg, phase center: 1.1411298385724182 deg -1.2635840686725994 deg Mean diff: -0.024722038614073977, Std diff: 0.05208547963589685
alt: 45.377620916048 deg, phase center: 356.6144894001537 deg -10.47967957817007 deg Mean diff: -0.030004875497425806, Std diff: 0.05289992383723178
alt: -16.999865394016005 deg, phase center: 257.1849045823964 deg -23.80743220393574 deg Mean diff: -0.003319799903335034, Std diff: 0.10915266755172588
alt: -16.817112437033597 deg, phase center: 257.26338326445006 deg -23.56798454591868 deg Mean diff: -0.003499723437833069, Std diff: 0.10988942743295312
alt: -17.466368571748365 deg, phase center: 256.84205252036094 deg -24.19717965052778 deg Mean diff: -0.0026703721158995108, Std diff: 0.10625577949146608
alt: -16.000295314055546 deg, phase center: 258.14092536218783 deg -23.32562920151148 deg Mean diff: -0.005081875688521128, Std diff: 0.117494907406081
alt: -15.385621889939655 deg, phase center: 258.04189759660414 deg -21.969838967355134 deg Mean diff: -0.004945661704209236, Std diff: 0.11654237682834045
alt: -16.140209723664153 deg, phase center: 255.84210751347544 deg -20.006281256701556 deg Mean diff: -0.0009574952564445145, Std diff: 0.0984002836613996
alt: -14.402022982582498 deg, phase center: 255.40734734738976 deg -15.983604017405893 deg Mean diff: 0.0014624907625040528, Std diff: 0.09043430430779925

dsa2000_jpl_results.json

[aelanman]

@Joshuaalbert Here's a conference paper I wrote last year that says more about pycalc11 https://ieeexplore.ieee.org/abstract/document/10464974

This is the main reference for the Consensus model used in CALC. In particular, the Eubanks paper starting on page 60: https://www.researchgate.net/profile/Marshall-Eubanks/publication/330766254_Proceedings_of_the_US_Naval_Observatory_Workshop_on_Relativistic_Models_for_use_in_Space_Geodesy/links/5c5394eea6fdccd6b5d87021/Proceedings-of-the-US-Naval-Observatory-Workshop-on-Relativistic-Models-for-use-in-Space-Geodesy.pdf

[Joshuaalbert]

Thanks @aelanman for the nice reference