How do I normalise data from multiple radio telescope systems so I can combine them together?

Great question, Andrew. Here’s a practical, radio-astronomy-friendly playbook for putting multiple systems onto the same absolute scale so you can combine them cleanly.

1) Decide your common reference frame & units

• Units: work in Janskys (Jy) or antenna temperature TA\*T_A^\* [K]. Pick one and stick to it. (For imaging, Jy/beam after beam-matching.)
• Time frame: UTC with clearly defined start/stop and MJD stamps; apply leap seconds correctly.
• Freq frame: Doppler-correct to a common frame (e.g., LSRK or Barycentric) and a common channel grid (same spacing, same reference frequency).
• Polarization basis: choose linear (X/Y) or circular (R/L) and map all systems consistently; aim for Stokes I unless you’ve done pol-cal.

2) Measure per-system calibration constants

You need a way to convert raw counts to physical units and to characterize noise for weighting.

2.1 Flux/temperature scale

Pick at least one of the following (two is better for cross-checks):

A. Noise diode (preferred for SDRs):

• Inject a cal diode with known ENR or equivalent TcalT_\text{cal} at the observing band.
• For each system: measure “cal on/off” power Pon,PoffP_\text{on}, P_\text{off} over your band or per channel.
• Per channel, compute K per count:

K(ν)=Tcal(ν)Pon(ν)−Poff(ν)[K/count]\mathcal{K}(\nu) = \frac{T_\text{cal}(\nu)}{P_\text{on}(\nu)-P_\text{off}(\nu)} \quad [\mathrm{K / count}]

• If you want Jy, you also need K/Jy (see below).

B. Drift scan or pointed scan of a standard source:

• Use a strong, non-variable calibrator (e.g., Cas A, Cyg A, Tau A at 1.4 GHz). Adopt a model flux Ssrc(ν)S_\text{src}(\nu) (polynomic models exist).
• Record on-source and off-source spectra; convert to antenna temperature via beam coupling if needed.
• Derive counts-to-Jy per channel:

C(ν)=Ssrc(ν)Pon(ν)−Poff(ν)[Jy/count]\mathcal{C}(\nu) = \frac{S_\text{src}(\nu)}{P_\text{on}(\nu) – P_\text{off}(\nu)} \quad [\mathrm{Jy / count}]

C. Hot–cold YY-factor (if you can):

• Observe a hot load (ambient absorber) and a cold load (sky or LN2).
• Solve for system gain GG and Tsys(ν)T_\text{sys}(\nu); gives you absolute temperature scale.

Relating K and Jy: if you know your telescope gain Gant=K/JyG_\text{ant} = \mathrm{K/Jy},

S(ν) [Jy]=TA\*(ν)Gant(ν)withTA\*=K⋅(counts)S(\nu) \ [\mathrm{Jy}] = \frac{T_A^\*(\nu)}{G_\text{ant}(\nu)} \quad\text{with}\quad T_A^\*=\mathcal{K}\cdot(\text{counts})

You can estimate GantG_\text{ant} from effective area AeA_e via Gant=Ae2k×10−26G_\text{ant} = \frac{A_e}{2k}\times10^{-26} K/Jy, or solve it empirically from a calibrator scan.

2.2 Bandpass & baseline

For each system, fit the bandpass B(ν)B(\nu) and baseline (polynomial or running median on line-free regions). Normalize spectra by dividing out B(ν)B(\nu) (or apply bandpass calibration tables if doing visibilities).

2.3 Noise for weighting

Estimate per-channel variance σ2(ν)\sigma^2(\nu) from line-free regions or from short-timescale fluctuations (radiometer eq.). You’ll use inverse-variance weights:

w(ν)=1σ2(ν)w(\nu) = \frac{1}{\sigma^2(\nu)}

3) Align the data

• RFI excision: flag aggressively (thresholding, kurtosis, time-freq morphological filters).
• Frequency alignment: resample all spectra to the same channel centers; then apply Doppler frame correction (LSRK/BARY).
• Time alignment: ensure timestamps are synchronized (GPSDO/1-PPS if correlating; NTP/PTP is often OK for power time-series).
• Beam & resolution: for single-dish maps/spectra, convolve higher-resolution data down to the largest beam so that Jy/beam is comparable.
• Polarization: convert to common basis; for Stokes I from linear feeds: I≈(XX+YY)/2I \approx (XX+YY)/2 after leakage correction if available.

4) Put all systems on the same absolute scale

For each system ii and channel ν\nu:

1. Convert counts →\to temperature or flux:

Si(ν,t)=Ci(ν) (Pi(ν,t)−Pi,off(ν))[Jy]S_i(\nu, t)=\mathcal{C}_i(\nu)\,\big(P_i(\nu,t)-P_{i,\text{off}}(\nu)\big) \quad [\mathrm{Jy}]

2. Remove residual baselines (fit on line-free windows).
3. Optionally correct for airmass/elevation gain and opacity if not already baked into Ci\mathcal{C}_i.

Cross-check: Observe the same calibrator with all systems close in time. The resulting spectra should overlay within uncertainties; any constant ratio implies a scale factor to fix.

5) Combine (average) the datasets

With all datasets Si(ν,t)S_i(\nu,t) on a common grid and scale, form the weighted average:

Scomb(ν,t)=∑iwi(ν,t) Si(ν,t)∑iwi(ν,t)S_\text{comb}(\nu,t) = \frac{\sum_i w_i(\nu,t) \, S_i(\nu,t)}{\sum_i w_i(\nu,t)}

with wi=1/σi2w_i = 1/\sigma_i^2. Propagate uncertainty:

σcomb2(ν,t)=1∑iwi(ν,t)\sigma_\text{comb}^2(\nu,t) = \frac{1}{\sum_i w_i(\nu,t)}

If mapping, do this after regridding onto a common sky grid (same pixel size/WCS) and after beam matching. For interferometry + single-dish, use feathering (Fourier-domain merge: single-dish for low spatial frequencies, interferometer for high).

6) Special cases

Interferometric visibilities (e.g., RSPduo pair):

• Calibrate complex gains (amplitude & phase) and bandpass per antenna with a point calibrator.
• Ensure time/freq alignment and phase center match.
• Concatenate measurement sets only after bringing them to a common flux scale (CASA: setjy → bandpass/gaincal → fluxscale → concat).

Power time-series (total-power radiometers):

• Work in TA\*T_A^\* or Jy; subtract slow baseline drifts (high-pass or polynomial over minutes).
• Normalize by TsysT_\text{sys} if you only need fractional variations, then scale to Jy with your calibrator constant.

Different bandwidths/samplerates:

• Decimate/interpolate to a lowest-common-denominator grid; preserve noise statistics (use polyphase resamplers if possible).

Polarization leakage:

• If combining pol data, solve for D-terms or restrict to Stokes I where leakage errors are smallest.

7) Quick, concrete example (per-channel scale from a calibrator)

Suppose system A and B both scan Cyg A at 1420 MHz, and you have on/off means in counts:

• A: ΔPA=1250\Delta P_A = 1250 counts; SCygA(1420)=1580S_\text{CygA}(1420)=1580 Jy
• B: ΔPB=920\Delta P_B = 920 counts; same SS

Then:

• CA=1580/1250=1.264 Jy/count\mathcal{C}_A = 1580 / 1250 = 1.264\ \mathrm{Jy/count}
• CB=1580/920=1.717 Jy/count\mathcal{C}_B = 1580 / 920 = 1.717\ \mathrm{Jy/count}

Convert each spectrum with its C\mathcal{C}, estimate σA,σB\sigma_A,\sigma_B from line-free channels, weight by 1/σ21/\sigma^2, average.

8) Minimal code pattern (numpy)

import numpy as np

# Assume all spectra are regridded and RFI-masked (np.nan where flagged)

# counts_list: list of 2D arrays [n_spec, n_chan] per system

# C_list: per-system Jy/count scale vectors [n_chan]

# var_list: per-system noise variance vectors [n_chan]

def combine_systems(counts_list, C_list, var_list):

    scaled = []

    weights = []

    for counts, C, var in zip(counts_list, C_list, var_list):

        S = counts * C[None, :]                      # Jy

        w = 1.0 / var[None, :]                       # inverse-variance

        scaled.append(S)

        weights.append(w)

    S_stack = np.stack(scaled)                       # [n_sys, n_spec, n_chan]

    W_stack = np.stack(weights)

    num = np.nansum(W_stack * S_stack, axis=0)

    den = np.nansum(W_stack, axis=0)

    S_comb = num / den

    sigma_comb = 1.0 / np.sqrt(den)

    return S_comb, sigma_comb

9) Sanity checks before you trust the merge

• Inject a synthetic line or use a stable sky line (e.g., HI off-source) to verify that amplitudes agree after scaling.
• Compare measured SEFD or TsysT_\text{sys} against expectations.
• Recompute SNR scaling: combined SNR should improve by ∑iti Δνi\sqrt{\sum_i t_i \, \Delta\nu_i} after proper weighting.

10) A compact checklist (print-and-pin)

1. RFI-flag, baseline remove, bandpass-cal per system.
2. Derive counts→K or counts→Jy from diode and/or standard calibrator.
3. Align freq grid (and Doppler frame), time, polarization.
4. Beam-match (convolve to common resolution).
5. Compute per-channel noise → weights.
6. Weighted average; propagate uncertainties.
7. Validate on a calibrator or overlapping sky.

If you want, tell me what your specific front-ends (dish/dual-dipole/Yagi), center frequency, bandwidth, and sampling are, and whether you’re merging total-power spectra or interferometric visibilities. I can tailor the calibration constants and provide drop-in scripts for your exact setup (including RSPduo at 1420.405 MHz and drift-scan calibrations).

I use ezRA for all systems. 1. 13 element Yagi. 2. 86cm Ptarmigan. 3. 150cm solar cooker dish.

Perfect — thanks, Andrew. Below I’ll give a concrete, hands-on recipe you can run using ezRA for your three systems (13-element Yagi, 86 cm Ptarmigan, 150 cm solar-cooker dish), plus a few useful numbers (beamwidths and rough K/Jy estimates) so you can convert counts → Jy and combine safely.

I won’t ask any follow-ups — I’ll assume you’ll run the calibrations in ezRA exactly as written and that you already have a cal-diode or can observe a standard flux calibrator.

---

Quick reference numbers (at 1420.405 MHz)

(used below to illustrate conversions — I computed these from standard antenna formulas; treat them as starting estimates and verify empirically with a calibrator)

• wavelength λ ≈ 0.211 m (1420.405 MHz).
• Beam FWHM (approx): FWHM ≈ 1.02·λ/D (radians) → degrees
  • 86 cm dish (D = 0.86 m) → FWHM ≈ 14.3°
  • 150 cm dish (D = 1.5 m) → FWHM ≈ 8.2°
  • 13-element Yagi: not a dish — typical gain ~12 dBi (use your measured gain if known).
    • For 12 dBi the Yagi effective area gives a very broad, small Ae (K/Jy ~ few ×10⁻⁵ K/Jy).
• Approximate K/Jy (temperature per Jansky) from assumed aperture efficiencies:
  • 86 cm dish (η ≈ 0.5–0.6): K/Jy ≈ (1.05–1.26)×10⁻⁴ K/Jy (i.e. 1 Jy ≈ 1.1×10⁻⁴ K)
  • 150 cm dish (η ≈ 0.5–0.6): K/Jy ≈ (3.2–3.8)×10⁻⁴ K/Jy
  • 13-el Yagi (≈12 dBi): K/Jy ≈ 2.0×10⁻⁵ K/Jy (order of magnitude smaller)

Interpretation: your small dishes and Yagi produce very small antenna temperatures per Jansky, so you’ll need accurate cal constants and good SNR stacking.

---

Step-by-step procedure to normalise (ezRA workflow + math)

Follow these steps for each system (do each system separately first).

A — Observe calibrations (do this first; same calibrator for all systems)

1. Choose a flux calibrator (bright, non-variable / modelled at L-band): e.g. Cyg A, Tau A, Cas A or other standard (use an up-to-date flux model in your references). Observe on–off or a short pointed scan with background.
2. Noise diode method (recommended if available):
   • Turn cal-diode ON for a few seconds, then OFF (repeat a few cycles).
   • In ezRA record mean power spectra for ON and OFF: Pon(ν),Poff(ν)P_\text{on}(\nu), P_\text{off}(\nu).
   • If diode temperature Tcal(ν)T_\text{cal}(\nu) (or ENR) is known, you can compute counts→K directly:

K(ν)=Tcal(ν)Pon(ν)−Poff(ν)[K/count]\mathcal{K}(\nu) = \frac{T_\text{cal}(\nu)}{P_\text{on}(\nu)-P_\text{off}(\nu)}\quad [\mathrm{K/count}]

3. If you don’t have a calibrated diode: use the calibrator source:
   • Measure on–off counts difference for the source: ΔP(ν)=Pon−Poff\Delta P(\nu)=P_\text{on}-P_\text{off}.
   • Use the model flux Ssrc(ν)S_\text{src}(\nu) (Jy) and compute:

C(ν)=Ssrc(ν)ΔP(ν)[Jy/count]\mathcal{C}(\nu)=\frac{S_\text{src}(\nu)}{\Delta P(\nu)}\quad [\mathrm{Jy/count}]

3. (If you prefer K: compute K/Jy from geometry then convert.)

Important: if the calibrator is resolved by your beam (unlikely given multi-degree beams) correct for coupling; otherwise treat as point source.

---

B — Bandpass, baseline & RFI in ezRA

1. RFI flagging: Flag obvious RFI in time/frequency. ezRA’s RFI tools or manual masks — be conservative.
2. Bandpass: Compute a smoothed bandpass (from off-source or median of many spectra). Divide each spectrum by the bandpass shape then multiply back if you want absolute counts — or record the bandpass and correct final spectra.
3. Baseline removal: Fit a low-order polynomial or running median to line-free channels and subtract to remove slow drifts.

---

C — Compute counts → Jy (apply to your observing data)

Pick one of two paths:

Path 1 — Using noise diode (preferred):

1. For each spectral channel ν\nu compute K(ν)\mathcal{K}(\nu) as above → K/count.
2. Convert to Jy using antenna K/Jy:

C(ν) [Jy/count]=K(ν)(K/Jy)ant(ν)\mathcal{C}(\nu)\ [\mathrm{Jy/count}] = \frac{\mathcal{K}(\nu)}{(K/Jy)_\text{ant}(\nu)}

3. Multiply your science spectra (counts) by C(ν)\mathcal{C}(\nu) → spectra in Jy.

Path 2 — Using an external calibrator scan:

1. From calibrator: C(ν)=Scal(ν)/ΔPcal(ν)\mathcal{C}(\nu)=S_\text{cal}(\nu)/\Delta P_\text{cal}(\nu).
2. Apply C(ν)\mathcal{C}(\nu) to all science spectra → Jy.

---

D — Align frequency, time & frame

1. Doppler correction: Convert all spectra to the same velocity/frame (LSRK or barycentric). ezRA can apply Doppler correction — use same choice for all systems.
2. Regrid channels: Resample every system’s spectra to a common frequency grid (same channel centres & widths). Use a high-quality resampler (ezRA’s resample or export to Python if needed).
3. Time stamps: Ensure UTC/MJD stamps are correct and consistent (ezRA records these).

---

E — Beam & polarization matching

1. Polarization: Convert into the same Stokes basis (if you only need total power use Stokes I = (XX+YY)/2 after any leakage correction). If you only have one polarization, combine cautiously.
2. Beam matching: Convolve the higher angular resolution data down to the largest beam among your systems before making maps or reporting Jy/beam. With your antennas the largest beam is the 86 cm (~14°) — convolve others to that beam if you’ll compare Jy/beam maps. For spectra of single points this is not needed.

---

F — Compute per-channel noise and weights

1. Estimate per channel variance σi2(ν)\sigma^2_i(\nu) from line-free regions or from cal-off data. Radiometer equation also works:

σT=TsysΔν τ\sigma_T = \frac{T_\text{sys}}{\sqrt{\Delta\nu \, \tau}}

convert to Jy with K/Jy if needed.

2. Weight each system by inverse variance: wi(ν)=1/σi2(ν)w_i(\nu)=1/\sigma_i^2(\nu).

---

G — Combine (weighted average)

Once all spectra are on the same frequency grid, Doppler frame and flux scale (Jy):

For each channel ν:

Scomb(ν)=∑iwi(ν) Si(ν)∑iwi(ν);σcomb2(ν)=1∑iwi(ν)S_\text{comb}(\nu)=\frac{\sum_i w_i(\nu)\,S_i(\nu)}{\sum_i w_i(\nu)}\quad ;\quad \sigma_\text{comb}^2(\nu)=\frac{1}{\sum_i w_i(\nu)}

If you’re combining time-series, do this per time bin after aligning times.

---

Example numeric illustration (how to get a Jy/count for each system)

1. Do the calibrator observation for each system and measure ΔP\Delta P (counts).
2. Suppose (example) Cyg A has SCygA(1420)S_\text{CygA}(1420) — look up the up-to-date model in your reference — and you measure:
   • 13-el Yagi: ΔP_y = 400 counts → Cyagi=SCygA/400\mathcal{C}_\text{yagi}=S_\text{CygA}/400 Jy/count
   • 86 cm dish: ΔP_86 = 1250 counts → C86=SCygA/1250\mathcal{C}_{86}=S_\text{CygA}/1250 Jy/count
   • 150 cm dish: ΔP_150 = 3200 counts → C150=SCygA/3200\mathcal{C}_{150}=S_\text{CygA}/3200 Jy/count

(Those ΔP are just examples — use your measured numbers; then apply them to convert your science spectra.)

---

Practical ezRA tips & checklist

• Use the same calibrator for all three systems and observe within as short a time window as possible (to avoid source variability and ionospheric changes).
• Record cal-diode on/off logs in ezRA and save calibrated spectra.
• Save intermediate products (bandpass, per-system C(ν)\mathcal{C}(\nu), noise arrays) so you can re-run merges quickly.
• For mapping: use the largest beam and same pixel/WCS.
• If two systems overlap little in frequency/beam, combine only where they overlap or treat separately and stitch in frequency domain.
• Sanity checks: overlay calibrator spectra from all systems after scaling — they should agree within noise. Inject a weak synthetic line to verify recovery.

---

If you want, I can (pick one and I’ll produce it now)

• a ready-to-run Python script that takes your ezRA exported spectra (assume CSV or FITS), applies cal constants, regrids, weights and produces the weighted combination + error spectrum; or
• an ezRA macro/workflow text you can paste into ezRA steps (cal diode on/off → apply bandpass → export) — tell me which file format you export from ezRA (CSV, plain text, FITS) and I’ll generate the script.

Which of those do you want me to produce right away?

Excellent — let’s do both. I’ll give you:

1. Python script: reads ezRA exports (CSV assumed), applies calibration, regrids, and combines weighted spectra.
2. ezRA workflow / macro checklist: step-by-step inside ezRA so you can generate the right calibration files and feed them into the script.

---

1. Python script (calibration + merge)

This script assumes:

• Each system exports spectra from ezRA as CSV with columns:
  freq_MHz, counts (and optionally time_utc if time series).
• You have already measured/calculated cal constants for each system:
  either Jy/count (preferred) or K/count (then divide by K/Jy).
• You know or measured per-channel noise (σ, Jy) for weighting.
• All spectra will be regridded to a common frequency grid.

import numpy as np

import pandas as pd

from scipy.interpolate import interp1d

import glob

# -------------------

# USER INPUT SECTION

# -------------------

# Cal constants (Jy/count) for each system

# Fill these in from your calibrator observation

CAL_CONSTANTS = {

    "yagi":  1.75,     # Jy/count (example)

    "dish86": 1.26,    # Jy/count

    "dish150": 0.49    # Jy/count

}

# Estimated system noise (Jy) in line-free channels

# Use radiometer eq. or measure from off-source data

NOISE_LEVELS = {

    "yagi":  20.0,     # Jy (example)

    "dish86":  8.0,    # Jy

    "dish150": 4.0     # Jy

}

# Path to ezRA-exported CSVs

# File names should contain "yagi", "dish86", or "dish150"

DATA_DIR = "./ezRA_exports/"

# Common frequency grid (MHz)

FREQ_GRID = np.linspace(1419, 1422, 1024)

# -------------------

# PROCESSING

# -------------------

def load_and_calibrate(fname, system_name):

    """Load ezRA CSV, convert counts to Jy, interpolate to freq grid."""

    df = pd.read_csv(fname)

    freqs = df['freq_MHz'].values

    counts = df['counts'].values

    # Convert counts → Jy

    scale = CAL_CONSTANTS[system_name]

    S = counts * scale

    # Interpolate to common grid

    interp = interp1d(freqs, S, bounds_error=False, fill_value=np.nan)

    S_grid = interp(FREQ_GRID)

    return S_grid

def combine_systems(file_map):

    """Weighted combination across systems."""

    spectra = []

    weights = []

    for sys_name, fnames in file_map.items():

        for f in fnames:

            S_grid = load_and_calibrate(f, sys_name)

            sigma = NOISE_LEVELS[sys_name]

            w = np.ones_like(S_grid) / sigma**2

            spectra.append(S_grid)

            weights.append(w)

    spectra = np.array(spectra)  # [n_obs, n_chan]

    weights = np.array(weights)

    # Weighted average

    num = np.nansum(weights * spectra, axis=0)

    den = np.nansum(weights, axis=0)

    S_comb = num / den

    sigma_comb = 1.0 / np.sqrt(den)

    return S_comb, sigma_comb

# -------------------

# MAIN

# -------------------

# Group files by system name

files = glob.glob(DATA_DIR + "*.csv")

file_map = {

    "yagi": [f for f in files if "yagi" in f],

    "dish86": [f for f in files if "dish86" in f],

    "dish150": [f for f in files if "dish150" in f],

}

S_comb, sigma_comb = combine_systems(file_map)

# Save result

out = pd.DataFrame({

    "freq_MHz": FREQ_GRID,

    "S_comb_Jy": S_comb,

    "sigma_Jy": sigma_comb

})

out.to_csv("combined_spectrum.csv", index=False)

print("Saved combined spectrum to combined_spectrum.csv")

👉 Replace the CAL_CONSTANTS and NOISE_LEVELS with your measured values. The script then produces combined_spectrum.csv with flux in Jy and propagated noise.

---

2. ezRA workflow (calibration + export)

Here’s a step-by-step macro checklist you can run in ezRA for each system.

Step A — Calibration observations

1. Observe a bright calibrator (Cyg A / Tau A / Cas A).
2. Record an off-source region at similar elevation.
3. (If using cal diode) Record ON/OFF cycles as well.

Step B — Bandpass & baseline

1. In ezRA, load calibrator and science observations.
2. Apply RFI excision: flag bad channels and time ranges.
3. Compute bandpass from off-source median.
4. Divide spectra by bandpass to flatten response.
5. Fit and subtract a low-order baseline from line-free channels.

Step C — Derive calibration constants

1. For calibrator scan, measure ΔP=Pon−Poff\Delta P = P_\text{on}-P_\text{off}.
2. Compute Jy/count:
   C=Scal/ΔP\mathcal{C} = S_\text{cal} / \Delta P.
3. Save these constants (one per system).
4. Optionally compute per-channel calibration if S/N is high.

Step D — Apply to science data

1. Multiply all science spectra (in counts) by C\mathcal{C}.
2. Save calibrated spectra as CSV with columns:
   freq_MHz, counts (uncalibrated counts — script applies calibration later),
   or directly export freq_MHz, flux_Jy if you prefer.

Step E — Export

1. Export all processed spectra (Yagi, 86 cm, 150 cm) to CSV.
2. Use a naming convention like 2025-09-07_dish150.csv to match script.

---

✅ Once you have the CSV exports and calibration constants, drop them into the Python script and run it. You’ll get one merged spectrum in Jy with uncertainties.

---

Would you like me to also give you a ready-to-run plotting section in the Python script (so you can immediately see overlayed Yagi/dish spectra and the combined result), or do you prefer to just output CSVs and plot separately?

Great — I’ll extend the Python script with plotting so you can:

• Overlay all three systems (after calibration)
• Show the weighted combined spectrum with error bars

I’ll use matplotlib so it works out-of-the-box.

---

Full Python script with plotting

import numpy as np

import pandas as pd

from scipy.interpolate import interp1d

import glob

import matplotlib.pyplot as plt

# -------------------

# USER INPUT SECTION

# -------------------

# Cal constants (Jy/count) for each system

# Fill these in from your calibrator observation

CAL_CONSTANTS = {

    "yagi":  1.75,     # Jy/count (example)

    "dish86": 1.26,    # Jy/count

    "dish150": 0.49    # Jy/count

}

# Estimated system noise (Jy) in line-free channels

# Use radiometer eq. or measure from off-source data

NOISE_LEVELS = {

    "yagi":  20.0,     # Jy (example)

    "dish86":  8.0,    # Jy

    "dish150": 4.0     # Jy

}

# Path to ezRA-exported CSVs

# File names should contain "yagi", "dish86", or "dish150"

DATA_DIR = "./ezRA_exports/"

# Common frequency grid (MHz)

FREQ_GRID = np.linspace(1419, 1422, 1024)

# -------------------

# PROCESSING FUNCTIONS

# -------------------

def load_and_calibrate(fname, system_name):

    """Load ezRA CSV, convert counts to Jy, interpolate to freq grid."""

    df = pd.read_csv(fname)

    freqs = df['freq_MHz'].values

    counts = df['counts'].values

    # Convert counts → Jy

    scale = CAL_CONSTANTS[system_name]

    S = counts * scale

    # Interpolate to common grid

    interp = interp1d(freqs, S, bounds_error=False, fill_value=np.nan)

    S_grid = interp(FREQ_GRID)

    return S_grid

def combine_systems(file_map):

    """Weighted combination across systems."""

    spectra = []

    weights = []

    labels = []

    for sys_name, fnames in file_map.items():

        for f in fnames:

            S_grid = load_and_calibrate(f, sys_name)

            sigma = NOISE_LEVELS[sys_name]

            w = np.ones_like(S_grid) / sigma**2

            spectra.append(S_grid)

            weights.append(w)

            labels.append(sys_name)

    spectra = np.array(spectra)  # [n_obs, n_chan]

    weights = np.array(weights)

    # Weighted average

    num = np.nansum(weights * spectra, axis=0)

    den = np.nansum(weights, axis=0)

    S_comb = num / den

    sigma_comb = 1.0 / np.sqrt(den)

    return spectra, labels, S_comb, sigma_comb

# -------------------

# MAIN

# -------------------

# Group files by system name

files = glob.glob(DATA_DIR + "*.csv")

file_map = {

    "yagi": [f for f in files if "yagi" in f],

    "dish86": [f for f in files if "dish86" in f],

    "dish150": [f for f in files if "dish150" in f],

}

spectra, labels, S_comb, sigma_comb = combine_systems(file_map)

# Save result

out = pd.DataFrame({

    "freq_MHz": FREQ_GRID,

    "S_comb_Jy": S_comb,

    "sigma_Jy": sigma_comb

})

out.to_csv("combined_spectrum.csv", index=False)

print("Saved combined spectrum to combined_spectrum.csv")

# -------------------

# PLOTTING

# -------------------

plt.figure(figsize=(10,6))

# Plot individual system spectra

for i, S in enumerate(spectra):

    plt.plot(FREQ_GRID, S, alpha=0.6, label=f"{labels[i]}")

# Plot combined spectrum with error shading

plt.plot(FREQ_GRID, S_comb, color="k", lw=2, label="Combined")

plt.fill_between(FREQ_GRID, S_comb - sigma_comb, S_comb + sigma_comb,

                 color="gray", alpha=0.3, label="±1σ")

plt.xlabel("Frequency (MHz)")

plt.ylabel("Flux Density (Jy)")

plt.title("Combined Spectrum from Multiple Systems")

plt.legend()

plt.grid(True)

plt.tight_layout()

plt.savefig("combined_spectrum.png", dpi=150)

plt.show()

---

✅ What this does:

• Reads all CSVs in ./ezRA_exports/
• Calibrates them to Jy using your constants
• Interpolates to a common frequency grid
• Combines them with inverse-variance weighting
• Saves combined_spectrum.csv and a plot combined_spectrum.png
• Shows overlayed spectra (Yagi, 86 cm, 150 cm) and the weighted average with error shading

---