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:
| kind | units | handling |
|---|---|---|
:pseudo | beam pixel channel count adu lambda klambda | registered dimensionless by the extension |
:nounits | _, "" | Unitful.NoUnits |
:unsupported | WU FU fu cy deg_2 sq_deg | clear error + suggested replacement |
:dimension | rad sr | dimensionless 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; directionJ2000/ICRS/B1950/APP/GALACTIC/ECLIPTIC/AZEL/AZELGEO/HADEC/ITRF; frequencyTOPO/GEO/BARY/LSRK/LSRD/GALACTO. - Velocity-frame constants copied verbatim from casacore
MeasTable.cc; every conversion cross-checked againstcasatools.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.J2000is treated asICRS(~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.ymlgains adocsjob (julia-actions/julia-docdeploy) — builds the Documenter site and deploys it (stable + dev + PR previews) via theDOCUMENTER_KEYsecret thatTagBot.ymlalready references — plus a coverage upload on thetestjob (julia-actions/julia-processcoverage+codecov/codecov-action)..github/workflows/CompatHelper.ymladded (the canonical Julia[compat]-bumping bot, root +docs/)..github/dependabot.ymltrimmed to thegithub-actionsecosystem only — CompatHelper now owns the Julia dependency updates.- README badge rows:
docs-stable/docs-devon the first, CI / codecov / MIT-license on the second. docs/make.jldeploydocsguarded onCI, withversions = ["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 theMEASINFOkeyword 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 theVarRefColcode column) via the column's own whole-column fast path.query/groupby/update!/delete!loaded the columns aWHERE/GROUP BY/SET-RHS touches as lazyColumns and indexedcol[i]per row — for anIncrementalStMancolumn (most MAIN metadata:TIME,FIELD_ID,SCAN_NUMBER, …) that re-walks the bucket index on every access. A new_load_colmaterialises 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 (aWHEREusually 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:
getcellon aDATAcell 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/dowcdate/ctime/cmonth/cdow/ctod(formatted strings)hms(rad)/dms(rad)— sexagesimal stringsnormangle(rad)— wrap to(−π, π]angdist/angdistx— great-circle distance,angdist(lon1, lat1, lon2, lat2)orangdist([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.Quantitycolumn (scalar or array-cell) isustripped to its first element's unit and getsQuantumUnitsstamped. A new_ms_ustring(inUnitfulExt) 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
Measurecolumn —MEpoch{R}→mjd*86400s,MDirection{R}→[lon,lat]rad,MPosition{R}→[x,y,z]m,MFrequency{R}→ Hz,MRadialVelocity{R}→ m/s — getsMEASINFO(type+ fixedReffrom the first row's frame) +QuantumUnits. - New
units = Dict(col => "Hz" | ["rad","rad"])kwarg onwrite_tablestampsQuantumUnitswithout aMEASINFO. An explicitmeasures=/units=entry for a column overrides the auto-detection (wins silently). _casatype_ofgains anAnyfallback with an actionable message (import Unitful/ check the eltype) instead of aMethodError.edit'saddcolumn!(t, name, data)does the same detection.copyms/write_msunchanged — 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 casacoreMCRadialVelocity. - Verified against
casatools.measures(): BARY / LSRD / GALACTO to < 1 mm/s; GEO / TOPO to ~0.25 m/s (the SOFAepv00vs casacore-ephemeris residual, the same one the frequency test carries asrtol = 2e-9). SOFA-onlyLSRK → R → LSRKround-trips to 1e-12. _MEAS_ENUM[:radialvelocity]added so a bare-codeVarRefColcolumn 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 casacoreMCDoppler's hub =RATIO(F = ν/ν₀) route formulas;BETA=RELATIVISTIC=v/c,Z=OPTICAL.doppler/frequency/radialvelocity/restfrequencymirror casacoreMFrequency::to{Doppler,Rest}/fromDoppler/MRadialVelocity::{to,from}Doppler— the frame of afromDopplerresult defaults toLSRK(casacore's choice).:dopplerMEASINFO read + write (_MEAS_FRAMES/_MEAS_ENUM/_measure_column_spec) — completeness; no real MS stores one.- Verified against
casatools.measures()tortol = 1e-12(casatools reports every dopplerm0asvalue·cin "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 axisshiftfreq(d, ν)—νa Hz number, anMFrequency(frame label kept), or a vector of either (the Doppler factor is computed once). Unlike casacore it converts a non-BETAdtoBETAfirst.- The two new
radialvelocity/frequencymethods compose the Phase 72 pieces so the commonMFrequency↔MRadialVelocitystep is one call; they broadcast, so no dedicatedVectormethod is needed for the bridges (onlyshiftfreqhas one, for the factor-once path). _beta_factorextracted and reused byfrequency/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)measurenow reads atype = "uvw"column asMuvW(wasMPosition) and atype = "baseline"column asMBaseline(was an error). Both frame sets mirrorMDirection::Types.MBaselineconversion = convert the unit direction via the existingMDirectioncode, rescale by the original length — matches casacoreMCBaselinefor every route (pure rotation, or itsadjust/readjustfor the aberration routes).MuvWconversion adds the phase-centre pole rotation (MCuvw::toPole/fromPole,R = RotMatrix(Euler(-π/2+lat, 2, -lon, 3))), so it also needsframe.direction.- Writing is automatic for an
MBaseline/MuvWcolumn (kind:baseline/:uvw,["m","m","m"]);copymsround-trips the keyword. - Non-goal (unchanged):
MEarthMagnetic, the IGRF field model,MVuvw(baseline, direction)construction, aVarRefColuvw/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 asMDirection{SUN}(the stored(lon, lat)is a placeholder); a per-rowVarRefColcode ≥ 32 decodes through_BODY_ENUM(40 => SUN,41 => MOON, …).measconvertcomputes the body's geocentric astrometric direction (light-time iterated), routed through the existing_icrs_to_dirchain for the target frame. The Moon's topocentric parallax is applied whenframe.positionis set (~1° for the Moon, ≲30″ planets).- Accuracy is the
plan94/moon98floor: ~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(noplan94),COMET/MeasCometephemeris 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_DIRin ITRF,[lon, lat]);mscal.uvw_j2000()([u, v, w]m — theUVWcolumn in J2000);mscal.delay()(geometric delay, s).1/2picksANTENNA1/ANTENNA2; no suffix uses antenna 0. - The tokenizer already lexed
mscal.ha1as one identifier (Phase 63); aTQLMScalAST node + a precompute hook in_tql_cols(coversquery+groupby) and_vtq_prepare!(VirtualTaQLColumn). Angle functions memo by(antenna, field, TIME). - Cross-checked against real
derivedmscal:ha1to < 1″,pa1/last1to < 1″ (pa1uses 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 aqueryclosure / ajoinonstring / on aGroupedTable.
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-handRX..YLoutputs are not supported.- Input basis comes from
POLARIZATION.CORR_TYPErow 1 (one setup per MS, as casacore's UDF does). The result is a(nOut, nchan)matrix; a(ncorr,)WEIGHT/SIGMAcell 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 (noLinearAlgebra);rescaleapplies casacore's 0.5 fudge factor to codes 5–12. - Threading mirrors
mscal.*: aTQLStokesnode, a"::stokes::…"sentinel ref split by_stokes_splitin_tql_cols(coversquery+groupby) and_vtq_prepare!. - Cross-checked against real
derivedmscalmscal.stokestortol = 1e-5. - Non-goals:
RX..YLinput frames; thePtotal/Plinear/ … pseudo outputs; per-DATA_DESC_IDCORR_TYPE;mscal.stokesin aqueryclosure /joinonstring.
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')")specis 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.baselineadditionally:L & R/L && R(baseline between two antenna sets —&drops autocorrelations,&&keeps them;L &repeatsL), and a whole-spec!negation.spwmaps the row'sDATA_DESC_IDthroughDATA_DESCRIPTION.SPECTRAL_WINDOW_ID.- A hand-written comma-list parser (not a port of casacore's per-type yacc grammars). Threads like
mscal.*— aTQLMSSelnode, a"::mssel::<fn>::<spec>"sentinel split by_mssel_splitin_tql_cols/_vtq_prepare!. - Cross-checked against real
derivedmscal(row counts match). - Non-goals:
mscal.time/mscal.uvdist; channel sub-selection onspw(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 oft0~t1ranges (or>t0/<t1). Each endpoint is an ISO (2024-05-24T10:00:00) orYYYY/MM/DD[/HH:MM:SS]datetime (parsed by the Phase-69_tql_parse_datetime), or a bare number = MJD days. Compared against the MAINTIMEcolumn (UTC seconds).mscal.uvdist('spec')— a comma-list ofa~branges (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 spwSPECTRAL_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_oneinsrc/taql/mscal.jl; threading unchanged (the"::mssel::<fn>::<spec>"sentinel from Phase 80). uvdistrow counts cross-checked against realderivedmscal;mscal.timeis 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 attenuationmscal.pbcorr(valexpr, 'spec' [, dir])/mscal.pbatten(valexpr, 'spec' [, dir])— pure parser sugar desugaring tovalexpr / mscal.pbresponse('spec', dir)/valexpr * mscal.pbresponse(...)(no new AST node, no new_mscal_columnsbranch —TQLArith's existing elementwise_bcastalready broadcasts the division/ multiplication over an array cell likeDATAagainst the scalar response). Usable anywhere an expression is, including anupdate!SET RHS to primary-beam-correct a column in place using the same per-rowTIME/ANTENNA1/FIELD_ID/POINTING.DIRECTIONgeometry asmscal.pbresponse.- Exact inverses of each other (up to storage precision) —
pbattenthenpbcorrround-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 themscal.*family (not a realderivedmscalUDF): theGaussianBeam/AiryBeampower response towarddir(defaultFIELD.PHASE_DIR, same direction-argument mini-language asmscal.azel1()etc.) as seen through ANTENNA1's actual pointing (POINTING.DIRECTION, matched by antenna + nearest-pastTIME) rather than its nominal position — the attenuation from a pointing/tracking error, computed automatically from the row'sTIME/ANTENNA1/FIELD_IDgeometry (needs aPOINTINGsubtable +import SOFA).- Both directions are compared in
AZEL(matching_cache's existing frame), so no extra conversion beyondPOINTING.DIRECTION's ownAZELGEO→AZELstep. - No cross-check oracle (a MeasurementSets-only extension); verified by a controlled, time-aligned synthetic
POINTINGfixture (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-liteEllipticalGaussianBeam(hpbw_major, hpbw_minor, pa, reffreq)— position-angle-rotated Gaussian power pattern (pafrom north through east, theMDirectionconvention); needs a 2-D(dlon, dlat)offset (a scalarθerrors clearly — ambiguous for a non-circular beam).SquintBeam(base, squint)offsets anyPrimaryBeam's centre by a fixed(dlon, dlat)— feed/pointing squint; composes with any base beam includingEllipticalGaussianBeam.pointing_offset(pointing::MDirection, target::MDirection) -> (dlon, dlat)— the small-angle tangent-plane offset (both directions in the same frame).- Every
PrimaryBeamnow 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 (noPrimaryBeamobject in TaQL) usable inquery/groupbyWHERE and computedselect.
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—PrimaryBeamabstract 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 viaSpecialFunctions.besselj1, already a dependency);PolynomialBeam(coeffs, maxrad, reffreq)(CASAPBMath1DPolyformpb = 1 + Σcₖ·(ν[GHz]·θ[arcmin])^(2k),0beyondmaxrad— 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_fluxare 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 toomscal.baseline/mscal.feed'sL & Rgrammar gainsL &&&— casacoreMSAntennaParse::AutoCorrOnly(self-correlations only;&=cross-only,&&=cross+auto, unchanged)._mssel_baseline_predchecks&&&before&&/&(a substring of both).mscal.baselinealso accepts a bare physical baseline-length range/bound with no&('100~500m'/'<200m'/'>1km', unitmdefault /km), computed fromANTENNA.POSITION— casacore'sblengthlist(distinct frommscal.uvdist, which is the per-row, frequency-dependentuvw).- Cross-checked against real
derivedmscal/tableCommandwhere 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
libmeasUDF 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;mjdMJD days for app/azel/hadec/itrf,x,y,zITRF 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 wrappingmeasconvert/_lst; args are ordinary expressions (columns, arithmetic). Needsimport 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),stecin TECU. Positive RM ⇔ field toward the observer.faraday_rotation(rm, freq)=rm · (c/freq)²(freqa number orMFrequency);derotate_angle(χ, rm, freq)=χ − Δχ.const RM_IONOSPHERE = 2.631e-6. All exported. Pure helpers are core;rotation_measuredefers the SOFA requirement toemm_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_uparsetry inUnitfulExt) now lexes as a quantity literal (casacoresimexpr unit), socol > 3 km/BETWEEN 1.4 GHz AND 1.5 GHzwork alongside the adjacent1.4GHzform. 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 klambdamscal.timenow implements casacore'sMSTimeParsegrammar: a single time (|TIME − t0| ≤ EXPOSURE/2),t0~t1, edge-inclusive[t0~t1], bufferedN[t0~t1],t0+dur,>t0/<t1. Each time isY/[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 ydatetimes still parse.mscal.uvdist('<expr>:P%')widens the range by ±P percent (casacoreuvwdistexpr COLON FNUMBER PERCENT); a bare value now needs a:P%to be a range.- Hand-computed (no
derivedmscalUDF 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 whenNUM_POLY > 0(or the cell is(2, n+1)withn > 1):dir = c[:,1] + Σ c[:,k]·dtᵏ,dt = epoch − FIELD.TIME(s) — casacoreMSFieldColumns::interpolateDirMeas. Withoutepoch, ordt ≈ 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 optionalDiskLong/DiskLatcolumns, great-circle (SLERP) interpolated between the bracketing rows (casacoreMeasComet::getDisk).Ephemerisgainsdisklon/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.sublatemm_lineofsight(dir, height, pos, epoch)/EarthMagneticMachine— port of casacoremeasures/Measures/EarthMagneticMachine. Intersects the line of sight todirwith a sphereheightm above the observer's geocentric radius, samples the bundled IGRF-14 field there (_earthfield_itrf), and projects onto the line of sight.dirmay be in any direction frame (rotated to ITRF viaepoch+pos). Real method inext/SOFAExt.jl; exportedEarthMagneticMachine,emm_lineofsight.- CASA cross-check:
test/measures_fixture.pyre-derives the same geometry withme+ numpy andme.earthmagneticat 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 theIGRFmodel frame.Measureunion,reftype,show, exported. - New
src/measures/earthfield.jl—earthfield(pos, epoch)evaluates the IGRF-14 geomagnetic field._earthfield_itrfis a direct port of casacore'sEarthField::calcFieldspherical-harmonic synthesis;_igrf_ghlinearly 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 — noSOFA. measconvert(viaext/SOFAExt.jl):MEarthMagnetic{A}rotates between direction frames like anMBaseline(rotate the unit direction, keep the length);MEarthMagnetic{IGRF}evaluates the model atframe.position/frame.epochthen rotates to the target.- Read (
measure(t, col)) + write (write_tableauto-stampsMEASINFOtype=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).