← older entries

Phase 59 — UPDATE SET (col, maskcol) = expr

update!'s set key can be a (datacol, maskcol) tuple (TaQL's SET (NAME, MASKNAME) = …). ("D", "M") => "dexpr" writes dexpr to D and a non-finite flag (!isfinite, element-wise) of the result to M; ("D", "M") => ("dexpr", "mexpr") writes an explicit mask expression to M instead. Either name may be a slice / mask target. taql's UPDATE … SET parses (D, M) = dexpr and (D, M) = (dexpr, mexpr). Since TaQL-lite has no masked-array expressions the data and mask are given explicitly — a documented divergence from casacore, which writes expr's own attached mask, so there is no cross-check. New nonfinite(x) / isnonfinite(x) function (element-wise !isfinite) is available in expressions generally.

Phase 60 — masked arrays in the expression engine

A TQLMArray value type (data + a true-means-masked-out Bool mask, casacore's MArray). V[boolexpr] in an expression now yields the whole cell with !boolexpr masked out (not a flat selection); marray(d, m) / arraydata(m) / arraymask(m) build and unpack one. Reductions (sum / mean / min / max / median / variance / stddev / rms / any / all / ntrue / nfalse) skip masked elements, nelements counts the unmasked ones, and arithmetic propagates the mask (union). So mean(V[!FLAG]), gmax(mean(V[FLAG])) (a masked reduction inside a g* aggregate), and SET (D, M) = V[goodcond] (writes V's data and the !goodcond mask to the companion column — the faithful form Phase 59 approximated) all work. A computed select column of masked arrays persists as plain data (arraymask(expr) for a separate mask column).

Phase 61 — SELECT expr AS (val, mask)

query's select ("valname", "maskname") => "expr" pair form (and taql's SELECT expr AS (v, m)) emits two output columns from a masked-array expression — the data and the mask. So query(t; select = [("D", "F") => "marray(DATA, FLAG)"]) reads a column together with its mask column, and taql(t, "SELECT V[V > 0] AS (D, M)") splits a masked selection into data + mask. A non-masked RHS gives a non-finite flag for the mask column (consistent with Phase 59). casacore has no on-disk column↔mask association — the mask is always named, so marray(col, maskcol) is "read a masked column".

Phase 62 — masked g* aggregates (first-class)

A g* aggregate with a masked-array argument now reduces over the group's unmasked elements: scalar g* (gmean / gsum / gmin / gmax / gmedian / gstddev / grms / gany / gall / gproduct) pools every row's unmasked elements into one flat reduction — gmean(V[!FLAG]) is the mean of every unflagged visibility in the group — and per-element gs* reduces per cell position over the rows where that cell is unmasked (gmeans(V[!FLAG]) = the per-cell mean spectrum, flagged cells ignored; an all-masked cell → NaN/0/false). The 14 Phase-52 gs* closures collapse to one generic _perelem_reduce sharing the non-s scalar reducer (so gproducts is a genuine per-cell product, no matmul). NOT now broadcasts, so V[!FLAG] and NOT FLAG negate an array-cell mask elementwise. casacore's own masked-g* semantics are murky (no reliable cross-check) — the pooled-unmasked interpretation is the clean one.

Phase 63 — string L. / R.-qualified non-equi join condition

join's on can now be a string — a TaQL-lite expression with L. / R. table-qualified column references, join(l, r; on = "L.TIME BETWEEN R.T0 AND R.T1"). It is the declarative form of the Phase-56 predicate closure: same nested loop, same O(N×M) cost, same unmatched join-type semantics. Must reference at least one L.<col> and one R.<col> (aggregates rejected); a bare-identifier on string is still the Phase-28 index-lookup join. The tokenizer now lexes one optional .suffix on an identifier. TaQL has no cartesian non-equi join, so no cross-check.

Phase 64 — src/ restructured into subdirectories

Pure refactor, no behaviour change. src/ is now organised into io/ (the aips.jl / lock.jl wire-format codecs), tables/ (the CTDS type system + table model + I/O: typeenum.jl, record.jl, table.jl, column.jl, interface.jl, writer.jl, create.jl, edit.jl, resync.jl), datamanagers/ (unchanged), and taql/ (below); only the module file and the two MS-domain files measurementset.jl / schema.jl stay at the top. MeasurementSets.jl's include list keeps its existing order (dependencies force table files to interleave with the datamanagers/ and taql/ includes), just gaining directory prefixes.

The ~2300-line src/query.jl and src/write_commands.jl become the src/taql/ directory: ast.jl (nodes + visitors + _bcast + masked arrays + indexing), parse.jl (tokenizer + parser + operator tables + pattern→regex), functions.jl (function / aggregate registries + _make_func), query.jl (ORDER BY + select + query(::AbstractTable)), groupby.jl (GroupSlice / GroupedTable / groupby + _geval), join.jl, and commands.jl (update! / delete! / insert! / taql). taql/taql.jl includes the first six; commands.jl stays included last (it builds on edit / create / resync). Test files renamed to taql_query_tests.jl / taql_command_tests.jl. The _select_spec helper also moved from table.jl to taql/query.jl, so the whole query engine — parser, evaluators, verbs, write commands — now lives under src/taql/; create.jl / resync.jl / edit.jl keep only thin GroupedTable / VirtualTaQLColumn dispatch adapters.

Phase 65 — physical units (Unitful weak-dependency extension)

import Unitful, UnitfulAngles, UnitfulAstro loads UnitfulExt, which maps a column's QuantumUnits keyword onto a Unitful unit:

columnunit(t, "CHAN_FREQ")   # u"Hz"
qcolumn(t, "UVW")            # the whole column as `… m` quantities (materialised)

UnitfulAngles supplies casacore's angle vocabulary (arcsec, mas, °), UnitfulAstro the astronomy units (Jy, pc, AU), and the extension registers the dimensionless "pseudo-units" casacore uses that no Julia package provides — beam, pixel, channel, count, adu, lambda, klambda (so Jy/beam parses to Jy beam⁻¹). The string-normalisation layer (_normalize_unit, always loaded) handles casacore's .-as-multiply, %/%%, and the FITS/long-form aliases.

MeasurementSets.UNITS_NO_JULIA_COUNTERPART documents every casacore unit without a third-party Julia implementation and how it is handled:

kindunitshandling
:pseudobeam pixel channel count adu lambda klambdaregistered dimensionless by the extension
:nounits_, ""Unitful.NoUnits
:unsupportedWU FU fu cy deg_2 sq_degclear error + suggested replacement
:dimensionrad srdimensionless here; base dimensions in casacore

Divergence from casacore: angles are SI-dimensionless (dimension(u"rad") == NoDims) — UnitfulAngles gives the unit names and all angle↔angle / angle↔scalar conversions, but not casacore's treatment of angle as a base dimension. DimensionfulAngles.jl is the strict alternative.

Read-side only for now. TaQL unit literals (3km, 10arcsec) and writing QuantumUnits from Unitful-typed columns are follow-ups.

Phase 66 — measures / reference frames (SOFA weak-dependency extension)

A measure-valued column declares its physical quantity and reference frame in a MEASINFO keyword. src/measures/ (always loaded) parses it and reads a cell as a typed value:

measinfo(t, "CHAN_FREQ")           # MeasInfo(:frequency, …, VarRefCol = "MEAS_FREQ_REF", …)
measure(main, "TIME", 1)           # MEpoch{UTC}(60454.4… d)
measure(subtable(ms, "FIELD"), "PHASE_DIR", 1)   # MDirection{J2000}(…, …)

import SOFA loads SOFAExt (pure-Julia SOFA.jl v2, IAU SOFA port), which converts between frames; import EarthOrientation additionally loads EarthOrientationExt for IERS ΔUT1 / polar motion:

fr = MeasFrame(epoch = measure(main, "TIME", 1),
               position = measure(subtable(ms, "ANTENNA"), "POSITION", 1),
               direction = measure(subtable(ms, "FIELD"), "PHASE_DIR", 1))

measconvert(measure(main, "TIME", 1), TAI)                 # epoch:  UTC → TAI
measconvert(MDirection{J2000}(2.0, 0.5), AZEL; frame = fr) # direction
measconvert(MFrequency{TOPO}(100e9), LSRK; frame = fr)     # frequency (radio/relativistic Doppler)
  • Epoch UTC/TAI/TT/TDB/UT1; direction J2000/ICRS/ B1950/APP/GALACTIC/ECLIPTIC/AZEL/AZELGEO/HADEC/ITRF; frequency TOPO/GEO/BARY/LSRK/LSRD/GALACTO.
  • Velocity-frame constants copied verbatim from casacore MeasTable.cc; every conversion cross-checked against casatools.measures() (epoch < 1 µs, direction < 5″, frequency < 1 Hz on 100 GHz).
  • Without EarthOrientation.jl: ΔUT1 = 0, no polar motion (~1″), a one-time warning. J2000 is treated as ICRS (~0.02″ frame bias).

Write path: write_table(dir, name, cols; measures = Dict("D" => (; kind = :direction, ref = "J2000"))) (or the per-row (; kind, varrefcol, tabtypes, tabcodes) form) stamps a MEASINFO keyword; copyms / copytable round-trip it verbatim.

Non-goals: solar-system-body direction frames (SUN/MOON/planets), MeasComet / ephemeris tables, MBaseline / MEarthMagnetic, standalone MDoppler, RefOff application, pulsar-timing-grade precision, TaQL measures functions, in-place MEASINFO edit.

Phase 67 — GitHub / CI wiring for the JuliaAstro repository

The package moved to github.com/JuliaAstro/MeasurementSets.jl. The stale Paul Barrett owner in the README CI badge and the docs/make.jl repo / edit links were fixed, and the git remote re-pointed.

  • .github/workflows/CI.yml gains a docs job (julia-actions/julia-docdeploy) — builds the Documenter site and deploys it (stable + dev + PR previews) via the DOCUMENTER_KEY secret that TagBot.yml already references — plus a coverage upload on the test job (julia-actions/julia-processcoverage + codecov/codecov-action).
  • .github/workflows/CompatHelper.yml added (the canonical Julia [compat]-bumping bot, root + docs/).
  • .github/dependabot.yml trimmed to the github-actions ecosystem only — CompatHelper now owns the Julia dependency updates.
  • README badge rows: docs-stable / docs-dev on the first, CI / codecov / MIT-license on the second.
  • docs/make.jl deploydocs guarded on CI, with versions = ["stable" => "v^", "v#.#", "dev" => "dev"], matching the JuliaAstro convention.

No source or test changes; the suite is unchanged at 2591.

Performance — bulk column reads in the query engine and measure

Two hot paths were re-decoding storage-manager cells one at a time inside a per-row loop:

  • measure(t, col) (whole-column) re-parsed the MEASINFO keyword and re-resolved the frame every row, and read the value cell by cell. Now it parses once and bulk-reads the value column (and the VarRefCol code column) via the column's own whole-column fast path.
  • query / groupby / update! / delete! loaded the columns a WHERE / GROUP BY / SET-RHS touches as lazy Columns and indexed col[i] per row — for an IncrementalStMan column (most MAIN metadata: TIME, FIELD_ID, SCAN_NUMBER, …) that re-walks the bucket index on every access. A new _load_col materialises each referenced scalar column once before the loop (array-valued columns stay lazy, so a predicate over a big cube column still streams). Orders-of-magnitude faster on a real-sized table; bounded extra memory (a WHERE usually references 1–3 scalar columns).

No API or behaviour change.

Performance — tiled reads no longer box every element

reinterpret(T, ::Vector{UInt8}) followed by per-element indexing is an allocating slow path in Julia when sizeof(T) > 1 (each access boxes a ReinterpretArray value). Every TiledStMan read went through it:

  • getcell on a DATA cell allocated ~144 KB (the cell is 2 KB); column(t, "DATA")[:] allocated ~22 KB/cell.

read_plane / read_cube_whole / _read_cube_bulk now copy each contiguous on-disk run with unsafe_load over a pinned pointer (_rd_run!, handling Real and Complex endianness). Measured on the fixture MS: DATA getcell 144 KB → ~3 KB/cell, UVW 2.0 → 0.9 KB, DATA[:] bulk 22 → 2.2 KB/cell. The RefTable / partial-rows= materialise path (always per-cell) benefits directly.

Also: the tile-layout computation (sortperm + allocation) is memoised per tileshape on the TiledStMan instead of recomputed on every read_plane; fixed-size index vectors became NTuples; and the SOFA extension uses StaticArrays SVector for its 3-vector / rotation math (StaticArrays added as a weak dependency — it is already a transitive dependency of SOFA.jl, so import SOFA still activates the extension).

No API, return-type, or behaviour change. 2592 tests.

Phase 69 — TaQL-lite quantity literals + date/time & angle functions

The string expression grammar (shared by query / groupby / join / update! / VirtualTaQLColumn CALC) gains two loosely-coupled features.

Quantity literals. 1.4GHz, 10arcsec, 30deg — a number token immediately followed (no space) by a unit. Compared against a column that carries a QuantumUnits keyword, the comparison goes through Unitful: query(spw, "CHAN_FREQ > 1.4GHz"). When an expression uses a quantity literal, every referenced unit-bearing column is loaded with its unit attached, so a bare-number comparison against such a column then raises a DimensionError (casacore's "units do not conform" — a deliberately stricter reading; only bites when bare and unit literals are mixed against one column). A computed select / CALC expression whose result is a dimensionless quantity ("CHAN_FREQ / 1GHz") is stripped to a plain number; a dimensional result errors. An update! SET RHS that evaluates to a quantity is converted to the target column's unit. Needs the Unitful extension (import Unitful, UnitfulAngles, UnitfulAstro); the parser errors clearly without it.

Date/time + angle functions (pure, no weak-dep gating — Dates stdlib + trig). Every date value is an MJD Float64 (days), so _bcast / isless / ORDER BY all keep working:

  • datetime('2020-02-12') / datetime() / mjd(x) / mjd() / mjdtodate(x) / date(x) / time(x)
  • year / month / day / week / weekday / dow
  • cdate / ctime / cmonth / cdow / ctod (formatted strings)
  • hms(rad) / dms(rad) — sexagesimal strings
  • normangle(rad) — wrap to (−π, π]
  • angdist / angdistx — great-circle distance, angdist(lon1, lat1, lon2, lat2) or angdist([lon1, lat1], [lon2, lat2])

Also: general array literals [a, b, c] in any expression position (not just IN); scientific-notation number literals (1.4e9, 1e-9) — previously an "unknown column" error.

Non-goals: measure literals / a TQLMeasure node (casacore core TaQL has none); mscal.azel() / measures-frame functions; the multi-point (2N-element) angdist array form; a BigFloat / Float128 high-precision MJD variant (considered, deferred); the spaced postfix 3 km unit form (no-space 3km only); input angle literals (30d15m as a value — only the hms(rad) function).

2638 tests.

Phase 70 — write QuantumUnits / MEASINFO from typed columns

The write side of Phases 65-66/69, previously asymmetric: a write_table column could only get a QuantumUnits / MEASINFO keyword from an explicit measures = Dict(...) kwarg with hand-written unit strings; a column whose Julia element type was already a Unitful.Quantity or a Measure hit a bare MethodError.

Now such a column is automatically stored as plain numbers on disk with the right keyword stamped, so write_table → readtable → qcolumn / measure round-trips:

  • A Unitful.Quantity column (scalar or array-cell) is ustripped to its first element's unit and gets QuantumUnits stamped. A new _ms_ustring (in UnitfulExt) is the inverse of _ms_uparse — a curated reverse-alias map + a round-trip check; a unit with no casacore spelling errors clearly (symmetric with the read side).
  • A Measure column — MEpoch{R} → mjd*86400 s, MDirection{R} → [lon,lat] rad, MPosition{R} → [x,y,z] m, MFrequency{R} → Hz, MRadialVelocity{R} → m/s — gets MEASINFO (type + fixed Ref from the first row's frame) + QuantumUnits.
  • New units = Dict(col => "Hz" | ["rad","rad"]) kwarg on write_table stamps QuantumUnits without a MEASINFO. An explicit measures= / units= entry for a column overrides the auto-detection (wins silently).
  • _casatype_of gains an Any fallback with an actionable message (import Unitful / check the eltype) instead of a MethodError.
  • edit's addcolumn!(t, name, data) does the same detection.
  • copyms / write_ms unchanged — keyword records are already copied verbatim.

Non-goals: deriving a per-row VarRefCol from a mixed-frame Measure column (fixed Ref only); non-canonical stored units for a Measure; _ms_ustring for the full casacore unit vocabulary (curated

  • verified, else a clear error); reading a plain column back as a

Quantity automatically (column still returns plain numbers — use qcolumn / measure).

2677 tests.

Phase 71 — MRadialVelocity reference-frame conversion

Phase 66 shipped measconvert for epoch / direction / frequency and left MRadialVelocity (m/s, SOURCE.SYSVEL) as a clear-error stub. Phase 71 implements it, reusing the frequency path's velocity machinery in ext/SOFAExt.jl (_v_earth_bary / _v_obs_geo / _n_hat / the _VEL_LSRK/_VEL_LSRD/_VEL_LSRGAL constants):

measconvert(MRadialVelocity{LSRK}(2e4), BARY; frame = fr)   # needs frame.direction
  • Frames LSRK / LSRD / BARY / GEO / TOPO / GALACTO, hub = BARY, relativistic velocity addition (β_out = (β_in ∓ g)/(1 ∓ β_in·g), g = V·n̂/c) — matches casacore MCRadialVelocity.
  • Verified against casatools.measures(): BARY / LSRD / GALACTO to < 1 mm/s; GEO / TOPO to ~0.25 m/s (the SOFA epv00 vs casacore-ephemeris residual, the same one the frequency test carries as rtol = 2e-9). SOFA-only LSRK → R → LSRK round-trips to 1e-12.
  • _MEAS_ENUM[:radialvelocity] added so a bare-code VarRefCol column decodes.

Non-goals: LGROUP / CMB velocity frames (no singletons — a measconvert targeting one errors clearly; no real MS SYSVEL uses them); a standalone MDoppler measure or an MFrequency ↔ MRadialVelocity bridge given a rest frequency (a distinct feature).

2699 tests.

Phase 72 — MDoppler measure + frequency ↔ velocity bridge

The Phase 71 non-goal: the Doppler-shift value — its convention (RADIO / OPTICAL / RATIO / BETA / GAMMA), orthogonal to the reference frame. Pure algebra (c is the only constant, exact SI) → core, no import SOFA.

d = MDoppler{RADIO}(0.01)
measconvert(d, OPTICAL)                        # convention -> convention

ν₀ = 1.42040575e9                              # HI rest frequency
doppler(MFrequency{LSRK}(1.4e9), ν₀)           # -> MDoppler{BETA}
frequency(d, ν₀)                               # MDoppler + rest -> MFrequency{LSRK}
restfrequency(MFrequency{LSRK}(1.4e9), d)      # -> MFrequency{REST}
radialvelocity(d)                              # -> MRadialVelocity{LSRK}  (c·β)
doppler(MRadialVelocity{LSRK}(3e5))            # -> MDoppler{BETA}         (β = v/c)
  • measconvert(::MDoppler, C2) uses casacore MCDoppler's hub = RATIO (F = ν/ν₀) route formulas; BETA = RELATIVISTIC = v/c, Z = OPTICAL.
  • doppler / frequency / radialvelocity / restfrequency mirror casacore MFrequency::to{Doppler,Rest} / fromDoppler / MRadialVelocity::{to,from}Doppler — the frame of a fromDoppler result defaults to LSRK (casacore's choice).
  • :doppler MEASINFO read + write (_MEAS_FRAMES / _MEAS_ENUM / _measure_column_spec) — completeness; no real MS stores one.
  • Verified against casatools.measures() to rtol = 1e-12 (casatools reports every doppler m0 as value·c in "m/s" — the fixture divides it back out).

Non-goals: MDoppler reference-frame handling (frame-agnostic — measconvert the frequency first); MDoppler::shiftFrequency (an array helper); a VarRefCol :doppler write path.

2725 tests.

Phase 73 — shared constants file

The physical / astronomical / calendar constants that were repeated across src/ and the extensions (speed of light — a bare literal in Phase 72 and SOFA.LIGHTSPEED in the SOFA extension; seconds- and milliseconds-per-day; the MJD reference epoch; one arcsecond; the AU; the casacore LSR-motion velocity vectors) are now defined once in src/constants.jl — MeasurementSets.C_LIGHT, SEC_PER_DAY, MSEC_PER_DAY, MJD_JD_OFFSET, MJD_EPOCH, ARCSEC, AU_METRES, VEL_LSRK / VEL_LSRD / VEL_LSRGAL. src/measures/, src/taql/, ext/SOFAExt.jl and ext/EarthOrientationExt.jl reference these (the extensions keep short local aliases). Every value is an exact SI / IAU definition and was verified equal to the SOFA.jl constant it replaced — behaviour-preserving, test count unchanged. Format-internal magic numbers (Dysco bit widths, CTDS type codes, container header offsets) stay with their subsystems.

2725 tests.

Phase 74 — shiftfreq + one-step frequency ↔ velocity bridges

The Phase 72 non-goals — casacore MDoppler::shiftFrequency and the spectral-axis array forms:

shiftfreq(d, νs)                              # νs .* √((1−β)/(1+β))
radialvelocity(f::MFrequency, ν₀)            # = radialvelocity(doppler(f, ν₀))
frequency(v::MRadialVelocity, ν₀)            # = frequency(doppler(v), ν₀)

radialvelocity.(measure(spw, "CHAN_FREQ", 1), ν₀)   # a whole velocity axis
  • shiftfreq(d, ν) — ν a Hz number, an MFrequency (frame label kept), or a vector of either (the Doppler factor is computed once). Unlike casacore it converts a non-BETA d to BETA first.
  • The two new radialvelocity / frequency methods compose the Phase 72 pieces so the common MFrequency ↔ MRadialVelocity step is one call; they broadcast, so no dedicated Vector method is needed for the bridges (only shiftfreq has one, for the factor-once path).
  • _beta_factor extracted and reused by frequency / restfrequency.

2735 tests.

Phase 75 — MBaseline + MuvW vector measures

The UVW column and antenna-to-antenna baselines as typed measures. Both are 3-vectors (metres) in a direction frame (ITRF default); measconvert rotates them.

measure(main, "UVW", 1)                        # -> MuvW{ITRF}(u, v, w m)
measconvert(MBaseline{ITRF}(x, y, z), J2000; frame = fr)
measconvert(measure(main, "UVW", 1), J2000; frame = fr)   # needs fr.direction

write_table(dir, "T", ["W" => [MuvW{J2000}(10.0, 20.0, 30.0)]]; nrow = 1)
  • measure now reads a type = "uvw" column as MuvW (was MPosition) and a type = "baseline" column as MBaseline (was an error). Both frame sets mirror MDirection::Types.
  • MBaseline conversion = convert the unit direction via the existing MDirection code, rescale by the original length — matches casacore MCBaseline for every route (pure rotation, or its adjust/readjust for the aberration routes).
  • MuvW conversion adds the phase-centre pole rotation (MCuvw::toPole / fromPole, R = RotMatrix(Euler(-π/2+lat, 2, -lon, 3))), so it also needs frame.direction.
  • Writing is automatic for an MBaseline / MuvW column (kind :baseline / :uvw, ["m","m","m"]); copyms round-trips the keyword.
  • Non-goal (unchanged): MEarthMagnetic, the IGRF field model, MVuvw(baseline, direction) construction, a VarRefCol uvw/baseline write path.

Phase 76 — solar-system-body direction reference frames

A FIELD.PHASE_DIR column can name a solar-system body as its frame (a moving target). Nine body singleton frames — SUN MOON MERCURY VENUS MARS JUPITER SATURN URANUS NEPTUNE — resolved via SOFA.plan94 / moon98 (no new dependency).

measure(fld, "PHASE_DIR", 1)                          # -> MDirection{SUN}(0.0, 0.0)
measconvert(MDirection{SUN}(0, 0), AZEL; frame = fr)  # the Sun's az/el
measconvert(MDirection{JUPITER}(0, 0), J2000; frame = fr)
  • measure() reads a body-frame column as MDirection{SUN} (the stored (lon, lat) is a placeholder); a per-row VarRefCol code ≥ 32 decodes through _BODY_ENUM (40 => SUN, 41 => MOON, …).
  • measconvert computes the body's geocentric astrometric direction (light-time iterated), routed through the existing _icrs_to_dir chain for the target frame. The Moon's topocentric parallax is applied when frame.position is set (~1° for the Moon, ≲30″ planets).
  • Accuracy is the plan94 / moon98 floor: ~arcsec for Sun / Moon / Venus / Mercury, ~20-40″ for Mars, ~arcmin for Jupiter / Saturn — fine for a pointing / scheduling reference, not astrometry.
  • Body frames are source-only: measconvert(d, SUN) errors.
  • Non-goals: PLUTO (no plan94), COMET / MeasComet ephemeris tables, a stored pointing offset in a body-frame cell, solar light deflection (< 2″, below the ephemeris floor).

Phase 77 — mscal.* derived-MS TaQL functions

The astronomy-value subset of casacore's derivedmscal UDFs in the TaQL-lite grammar — per-MAIN-row values from TIME + the ANTENNA / FIELD subtables via the Measures engine (needs import SOFA).

query(main, "mscal.el1() > 0.35")                       # elevation cut
groupby(main, "FIELD_ID"; select = ["az" => "gmean(mscal.az1())"])
query(main, ""; select = ["w" => "mscal.uvw_j2000()", "d" => "mscal.delay()"])
  • Functions: mscal.ha1() / ha2() / ha() (hour angle, rad); mscal.hadec1() ([ha, dec]); mscal.azel1() ([az, el]); mscal.az1() / el1() (scalar); mscal.pa1() (parallactic angle); mscal.last1() (local apparent sidereal time, rad — casacore returns a raw MVEpoch day count here; we return the angle); mscal.itrf() (PHASE_DIR in ITRF, [lon, lat]); mscal.uvw_j2000() ([u, v, w] m — the UVW column in J2000); mscal.delay() (geometric delay, s). 1 / 2 picks ANTENNA1 / ANTENNA2; no suffix uses antenna 0.
  • The tokenizer already lexed mscal.ha1 as one identifier (Phase 63); a TQLMScal AST node + a precompute hook in _tql_cols (covers query + groupby) and _vtq_prepare! (VirtualTaQLColumn). Angle functions memo by (antenna, field, TIME).
  • Cross-checked against real derivedmscal: ha1 to < 1″, pa1 / last1 to < 1″ (pa1 uses casacore's HADEC-pole, not the J2000 pole).
  • Non-goals: the CASA-MSSelection selection functions (mscal.baseline / mscal.spw / …), the Observatories-table array centre, mscal.* in a query closure / a join on string / on a GroupedTable.

Phase 78 — mscal.stokes() polarization conversion

mscal.stokes(col [, 'types'] [, rescale]) — a TaQL-lite function that converts a DATA / FLAG / WEIGHT array cell between correlation bases, a port of casacore's StokesConverter.

query(main, "any(abs(mscal.stokes(DATA, 'I')) > 5.0)")
groupby(main, "FIELD_ID"; select = ["p" => "gmean(mean(abs(mscal.stokes(DATA, 'I'))))"])
  • types (default 'IQUV') is an alias (IQUV / STOKES, CIRC / CIRCULAR, LIN / LINEAR) or a comma-list of Stokes / circular / linear names ('I', 'I,V', 'XX,YY'). Mixed-hand RX..YL outputs are not supported.
  • Input basis comes from POLARIZATION.CORR_TYPE row 1 (one setup per MS, as casacore's UDF does). The result is a (nOut, nchan) matrix; a (ncorr,) WEIGHT / SIGMA cell yields a (nOut,) vector.
  • Per-cell: a matrix multiply for Complex data (out[o,ch] = Σⱼ conv[o,j]·in[j,ch]), an any-of-contributing test for Bool flags, the weight-propagation formula for Float weights. The 6 basis-change 4×4 matrices are hardcoded (no LinearAlgebra); rescale applies casacore's 0.5 fudge factor to codes 5–12.
  • Threading mirrors mscal.*: a TQLStokes node, a "::stokes::…" sentinel ref split by _stokes_split in _tql_cols (covers query + groupby) and _vtq_prepare!.
  • Cross-checked against real derivedmscal mscal.stokes to rtol = 1e-5.
  • Non-goals: RX..YL input frames; the Ptotal / Plinear / … pseudo outputs; per-DATA_DESC_ID CORR_TYPE; mscal.stokes in a query closure / join on string.

Phase 79 — mscal.uvw_j2000() rotation memo

The ITRF→J2000 uvw transform is a linear map (pole rotation + baseline rotation) that depends only on the frame — i.e. on (antenna, field, TIME). Phase 77 evaluated it with one measconvert(::MuvW) per MAIN row; Phase 79 memoizes the 3×3 (as its three result columns, three measconverts per distinct key) and applies it to each row's stored UVW. On a real MAIN this is one matrix per integration-baseline pair instead of one per row. No behaviour change — cross-checked against a direct per-row measconvert to rtol = 1e-9.

Phase 80 — mscal.<sel>() MSSelection-lite selection functions

The derivedmscal row-selection UDFs in the TaQL-lite grammar — mscal.baseline / mscal.field / mscal.spw / mscal.scan / mscal.state / mscal.array / mscal.obs, each mscal.<sel>('spec') returning a per-MAIN-row Bool.

query(main, "mscal.baseline('ea01 & *') AND mscal.field('3C286')")
query(main, "mscal.spw('0~3') AND NOT mscal.scan('1,2')")
  • spec is a comma-list of terms; a row passes if it matches ANY. Terms: N, N~M (inclusive range), >N / <N / >=N / <=N, an exact name, a glob (* ? [...]), or /regex/ — names matched against the type's NAME column (ANTENNA.NAME, FIELD.NAME, SPECTRAL_WINDOW.NAME, STATE.OBS_MODE). A !-prefixed term is subtracted (a spec of only !-terms selects everything else).
  • mscal.baseline additionally: L & R / L && R (baseline between two antenna sets — & drops autocorrelations, && keeps them; L & repeats L), and a whole-spec ! negation. spw maps the row's DATA_DESC_ID through DATA_DESCRIPTION.SPECTRAL_WINDOW_ID.
  • A hand-written comma-list parser (not a port of casacore's per-type yacc grammars). Threads like mscal.* — a TQLMSSel node, a "::mssel::<fn>::<spec>" sentinel split by _mssel_split in _tql_cols / _vtq_prepare!.
  • Cross-checked against real derivedmscal (row counts match).
  • Non-goals: mscal.time / mscal.uvdist; channel sub-selection on spw (0:5~20); mscal.corr / mscal.feed.

Phase 81 — mscal.time() / mscal.uvdist() selection

The last two derivedmscal selection UDFs.

query(main, "mscal.time('2024/05/24/10:00:00~2024/05/24/11:00:00')")
query(main, "mscal.uvdist('20~200klambda') AND NOT mscal.uvdist('<50m')")
  • mscal.time('spec') — a comma-list of t0~t1 ranges (or >t0 / <t1). Each endpoint is an ISO (2024-05-24T10:00:00) or YYYY/MM/DD[/HH:MM:SS] datetime (parsed by the Phase-69 _tql_parse_datetime), or a bare number = MJD days. Compared against the MAIN TIME column (UTC seconds).
  • mscal.uvdist('spec') — a comma-list of a~b ranges (or >a / <b) with a trailing unit: m (default) / km / lambda / klambda / mlambda. Uses the 2-D uv-distance √(u²+v²) (casacore's fast path); wavelength units scale per row by the row's spw SPECTRAL_WINDOW.REF_FREQUENCY. A spec must not mix distance and wavelength units; a bare single value (no range/bound) is an error.
  • Both extend _mssel_one in src/taql/mscal.jl; threading unchanged (the "::mssel::<fn>::<spec>" sentinel from Phase 80).
  • uvdist row counts cross-checked against real derivedmscal; mscal.time is hand-computed only (casacore's time-string defaults make an exact edge match fiddly).
  • Non-goals: casacore's full time grammar (* wildcards, + durations, [...] edge buffers, MS-derived field defaults); the :P% percent-tolerance on a uvdist value.

Phase 102 — mscal.pbcorr() / mscal.pbatten(): primary-beam write path

update!(ms; set = ["DATA" => "mscal.pbcorr(DATA, 'gaussian:0.008727')"])   # true flux
update!(ms; set = ["DATA" => "mscal.pbatten(DATA, 'airy:25.0:8.0e9')"])   # simulate attenuation
  • mscal.pbcorr(valexpr, 'spec' [, dir]) / mscal.pbatten(valexpr, 'spec' [, dir]) — pure parser sugar desugaring to valexpr / mscal.pbresponse('spec', dir) / valexpr * mscal.pbresponse(...) (no new AST node, no new _mscal_columns branch — TQLArith's existing elementwise _bcast already broadcasts the division/ multiplication over an array cell like DATA against the scalar response). Usable anywhere an expression is, including an update! SET RHS to primary-beam-correct a column in place using the same per-row TIME/ANTENNA1/FIELD_ID/POINTING.DIRECTION geometry as mscal.pbresponse.
  • Exact inverses of each other (up to storage precision) — pbatten then pbcorr round-trips a value.

Phase 101 — mscal.pbresponse(): primary beam ↔ mscal.* integration

query(main, "mscal.pbresponse('gaussian:0.008727') < 0.5")   # tracking-error cut
groupby(main, "ANTENNA1"; select = ["a"=>:ANTENNA1, "m"=>"gmean(mscal.pbresponse('airy:25.0:8e9'))"])
  • mscal.pbresponse('gaussian:HPBW' | 'airy:D:FREQ[:BLOCKAGE]' [, dir]) — a MeasurementSets extension to the mscal.* family (not a real derivedmscal UDF): the GaussianBeam / AiryBeam power response toward dir (default FIELD.PHASE_DIR, same direction-argument mini-language as mscal.azel1() etc.) as seen through ANTENNA1's actual pointing (POINTING.DIRECTION, matched by antenna + nearest-past TIME) rather than its nominal position — the attenuation from a pointing/tracking error, computed automatically from the row's TIME/ANTENNA1/FIELD_ID geometry (needs a POINTING subtable + import SOFA).
  • Both directions are compared in AZEL (matching _cache's existing frame), so no extra conversion beyond POINTING.DIRECTION's own AZELGEO → AZEL step.
  • No cross-check oracle (a MeasurementSets-only extension); verified by a controlled, time-aligned synthetic POINTING fixture (offset 0 → response ≈ 1; a known offset → the closed-form Gaussian/Airy value).

Phase 100 — elliptical / squinted primary beams + TaQL-lite

θ = pointing_offset(pointing, source)                          # (dlon, dlat), same frame
attenuate(EllipticalGaussianBeam(hmaj, hmin, pa, freq), flux, θ)
power_response(SquintBeam(base_beam, (dlon0, dlat0)), θ)
query(cat, "pbairy(OFFSET, 25.0, 1.4e9) > 0.5")                 # TaQL-lite
  • EllipticalGaussianBeam(hpbw_major, hpbw_minor, pa, reffreq) — position-angle-rotated Gaussian power pattern (pa from north through east, the MDirection convention); needs a 2-D (dlon, dlat) offset (a scalar θ errors clearly — ambiguous for a non-circular beam). SquintBeam(base, squint) offsets any PrimaryBeam's centre by a fixed (dlon, dlat) — feed/pointing squint; composes with any base beam including EllipticalGaussianBeam.
  • pointing_offset(pointing::MDirection, target::MDirection) -> (dlon, dlat) — the small-angle tangent-plane offset (both directions in the same frame).
  • Every PrimaryBeam now accepts either a scalar θ or a (dlon, dlat) pair everywhere (power_response, voltage_response, attenuate, correct_flux) — a circularly symmetric beam falls back to the pair's magnitude.
  • TaQL-lite: pbgaussian(θ, hpbw), pbairy(θ, diameter, freq[, blockage]), pbellipse(dlon, dlat, hpbw_major, hpbw_minor, pa) — pure-numeric wrappers (no PrimaryBeam object in TaQL) usable in query/groupby WHERE and computed select.

Phase 99 — analytic primary-beam models

pb = GaussianBeam(1.4e9; diameter = 25.0)          # HPBW = 1.02λ/D
θ = angular_separation(pointing, source)            # both MDirection, same frame
correct_flux(pb, apparent_flux, θ)                  # -> true flux
power_response(AiryBeam(25.0; blockage=2.5), θ, 1.4e9)
  • New src/beam/beam.jl — PrimaryBeam abstract type; GaussianBeam(freq; diameter, k=1.02) (HPBW power pattern, freq scaling); AiryBeam(diameter; blockage=0.0) (uniformly illuminated circular aperture, optional central obstruction — the standard two-term closed form via SpecialFunctions.besselj1, already a dependency); PolynomialBeam(coeffs, maxrad, reffreq) (CASA PBMath1DPoly form pb = 1 + Σcₖ·(ν[GHz]·θ[arcmin])^(2k), 0 beyond maxrad — no coefficient table is bundled, none is vendored on this machine; supply your own).
  • power_response(beam, θ, freq=reffreq(beam)) is the one method each subtype implements; voltage_response, attenuate, correct_flux are generic over it. angular_separation(d1::MDirection, d2::MDirection) (great-circle, both directions in the same frame) — the θ input.
  • Standalone Julia feature, not a TaQL-lite integration and not a casacore port (Gaussian/Airy are textbook optics; no CASA oracle) — verified by half-power-point, Airy-null, and dish-scaling sanity checks, not a cross-check.

Phase 98 — mscal.baseline &&& + physical baseline-length selection

query(main, "mscal.baseline('ea01 &&&')")     # self-correlations only
query(main, "mscal.baseline('100~500m')")     # physical antenna-pair distance
query(main, "mscal.feed('0 &&&')")            # &&& works on mscal.feed too
  • mscal.baseline / mscal.feed's L & R grammar gains L &&& — casacore MSAntennaParse::AutoCorrOnly (self-correlations only; &=cross-only, &&=cross+auto, unchanged). _mssel_baseline_pred checks &&& before &&/& (a substring of both).
  • mscal.baseline also accepts a bare physical baseline-length range/bound with no & ('100~500m' / '<200m' / '>1km', unit m default / km), computed from ANTENNA.POSITION — casacore's blengthlist (distinct from mscal.uvdist, which is the per-row, frequency-dependent uvw).
  • Cross-checked against real derivedmscal/tableCommand where registered.

Phase 97 — meas.* measure conversions in TaQL-lite

query(t, "meas.galactic(RA, DEC)[2] > 0")                       # b > 0
query(t, "meas.azel(RA, DEC, TIME/86400, X, Y, Z)[2] > 0.3")    # elevation
query(t, "meas.epoch('TAI', TIME/86400)")                       # UTC → TAI MJD
query(t, "meas.last(TIME/86400, X, Y, Z)")                      # local apparent sidereal time
  • A subset of casacore's libmeas UDF library: meas.<frame>(['SRC', ]lon, lat[, mjd[, x, y, z]]) converts a direction (<frame> = j2000/b1950/app/galactic/ecliptic/ azel/hadec/itrf/icrs; optional string-literal source frame, default J2000; mjd MJD days for app/azel/hadec/itrf, x,y,z ITRF m also for azel/hadec/itrf) → [lon, lat] rad. meas.epoch(scale, mjd) converts an epoch's time scale; meas.last/meas.lst(mjd, x, y, z) gives the local apparent sidereal time (rad).
  • Plain TQLFuncs wrapping measconvert / _lst; args are ordinary expressions (columns, arithmetic). Needs import SOFA.

Phase 96 — ionospheric Faraday rotation

m = EarthMagneticMachine(350e3, observatory("VLA"), MEpoch{UTC}(mjd))
rm = rotation_measure(m, MDirection{J2000}(ra, dec); stec = 12.0)   # rad/m², STEC in TECU
Δχ = faraday_rotation(rm, 1.4e9)                                     # RM·λ²
χ_true = derotate_angle(χ_obs, rm, MFrequency{TOPO}(1.4e9))
  • rotation_measure(dir, epoch, pos; stec, height=350e3) / rotation_measure(m::EarthMagneticMachine, dir; stec) — thin-shell ionospheric RM (rad/m²): RM_IONOSPHERE · stec · B∥, B∥ the line-of-sight field along the propagation direction at the shell pierce point (emm_lineofsight), stec in TECU. Positive RM ⇔ field toward the observer.
  • faraday_rotation(rm, freq) = rm · (c/freq)² (freq a number or MFrequency); derotate_angle(χ, rm, freq) = χ − Δχ.
  • const RM_IONOSPHERE = 2.631e-6. All exported. Pure helpers are core; rotation_measure defers the SOFA requirement to emm_lineofsight.

Phase 95 — spaced unit literals + observatory() in TaQL-lite

query(spw, "CHAN_FREQ > 1.4 GHz")                        # spaced postfix unit
query(ant, "sqrt(sum((POSITION - observatory('VLA'))**2)) < 1e5")
  • A number followed by a space then a bare identifier that is not a column and is a known unit (_tql_known_unit — a common-unit set in core, a full _ms_uparse try in UnitfulExt) now lexes as a quantity literal (casacore simexpr unit), so col > 3 km / BETWEEN 1.4 GHz AND 1.5 GHz work alongside the adjacent 1.4GHz form. A trailing non-unit ident stays an unknown-column error.
  • observatory('NAME') — a TaQL-lite function returning a telescope's ITRF [x, y, z] (m) from the bundled Observatories table; unknown name errors. Composes with array arithmetic / indexing.

Phase 94 — full MSSelection time grammar + uvdist :P%

query(main, "mscal.time('2024/05/24/09:00:00~11:00:00')")   # t1 inherits t0's date
query(main, "mscal.time('[09:00:00~11:00:00]')")            # edge-inclusive
query(main, "mscal.time('09:00:00 + 02:00:00')")            # t0 .. t0+2h
query(main, "mscal.time('2024/05/24/10:08:00')")            # ± EXPOSURE/2
query(main, "mscal.uvdist('100klambda:10%')")               # 90–110 klambda
  • mscal.time now implements casacore's MSTimeParse grammar: a single time (|TIME − t0| ≤ EXPOSURE/2), t0~t1, edge-inclusive [t0~t1], buffered N[t0~t1], t0+dur, >t0 / <t1. Each time is Y/[M/[D/]][h:[m:[s]]] with * wildcards; a missing component defaults to the first MAIN-row TIME (a ~ range's upper bound inherits from the lower). ISO / d U y datetimes still parse.
  • mscal.uvdist('<expr>:P%') widens the range by ±P percent (casacore uvwdistexpr COLON FNUMBER PERCENT); a bare value now needs a :P% to be a range.
  • Hand-computed (no derivedmscal UDF for the exotic time forms).

Phase 93 — polynomial PHASE_DIR + ephemeris sub-Earth point

  • measure(fld, "PHASE_DIR", row; epoch) now evaluates a FIELD direction as a time polynomial when NUM_POLY > 0 (or the cell is (2, n+1) with n > 1): dir = c[:,1] + Σ c[:,k]·dtᵏ, dt = epoch − FIELD.TIME (s) — casacore MSFieldColumns::interpolateDirMeas. Without epoch, or dt ≈ 0, the 0-order term (unchanged). mscal.* memoises a polynomial field per (field, TIME) like an ephemeris field.
  • ephemeris_diskpos(e, mjd) -> (lon, lat) — the sub-observer point on a body's surface from an ephemeris table's optional DiskLong / DiskLat columns, great-circle (SLERP) interpolated between the bracketing rows (casacore MeasComet::getDisk). Ephemeris gains disklon / disklat. Exported.

Phase 92 — EarthMagneticMachine (line-of-sight field)

m = EarthMagneticMachine(350e3, observatory("VLA"), MEpoch{UTC}(mjd))
r = m(MDirection{J2000}(ra, dec))
r.losfield        # nT parallel to the line of sight (× slant TEC × 2.63e-13 → RM)
r.field, r.subpoint, r.sublon, r.sublat
  • emm_lineofsight(dir, height, pos, epoch) / EarthMagneticMachine — port of casacore measures/Measures/EarthMagneticMachine. Intersects the line of sight to dir with a sphere height m above the observer's geocentric radius, samples the bundled IGRF-14 field there (_earthfield_itrf), and projects onto the line of sight. dir may be in any direction frame (rotated to ITRF via epoch + pos). Real method in ext/SOFAExt.jl; exported EarthMagneticMachine, emm_lineofsight.
  • CASA cross-check: test/measures_fixture.py re-derives the same geometry with me + numpy and me.earthmagnetic at the pierce point — the geometry matches exactly, the field to ~3% (IGRF-12 vs -14).

Phase 91 — MEarthMagnetic measure + IGRF-14 model

earthfield(MPosition{ITRF}(x, y, z), MEpoch{UTC}(mjd))   # -> MEarthMagnetic{ITRF}, nT
measconvert(MEarthMagnetic{IGRF}(0, 0, 1e-6), J2000; frame)
  • New MEarthMagnetic{R} measure (a 3-vector in nano-tesla, in a direction-family frame) and the IGRF model frame. Measure union, reftype, show, exported.
  • New src/measures/earthfield.jl — earthfield(pos, epoch) evaluates the IGRF-14 geomagnetic field. _earthfield_itrf is a direct port of casacore's EarthField::calcField spherical-harmonic synthesis; _igrf_gh linearly interpolates the bundled 5-year Gauss coefficients (src/measures/igrf14_data.jl — IAGA / NOAA IGRF-14, degree 13, epochs 1900.0–2025.0 plus the 2025–2030 secular variation). Pure arithmetic — no SOFA.
  • measconvert (via ext/SOFAExt.jl): MEarthMagnetic{A} rotates between direction frames like an MBaseline (rotate the unit direction, keep the length); MEarthMagnetic{IGRF} evaluates the model at frame.position / frame.epoch then rotates to the target.
  • Read (measure(t, col)) + write (write_table auto-stamps MEASINFO type=earthmagnetic, QuantumUnits=["nT","nT","nT"]) round-trip.
  • CASA me.earthmagnetic('igrf') cross-check (test/measures_fixture.py): casacore ships IGRF-12, so components agree to ~100–200 nT (model generation) — the frame rotation is exact (magnitude ITRF ≡ J2000).

newer entries →