Addressed referee 2, report 1, comments 4, 7

import os

import glob

import numpy as np

import matplotlib.pyplot as plt

from spectral_cube import SpectralCube

from regions import CircleSkyRegion

import astropy.units as u

import warnings

# Suppress verbose astropy WCS warnings

warnings.filterwarnings('ignore')

def compare_tp_feather_spectralcube(source, species, ap_radius_arcsec=40.0, data_dir="assets/data",

                                    output_dir="stg/figures"):

    os.makedirs(output_dir, exist_ok=True)

    feather_files = glob.glob(f"{data_dir}/{source}_{species}_*feather_pbcor.fits")

    tp_files = glob.glob(f"{data_dir}/*TP*{source}_{species}*.fits")

    if not feather_files or not tp_files:

        print(f"[Error] Missing files for {source} - {species}.")

        return

    feather_path, tp_path = feather_files[0], tp_files[0]

    print(f"\n--- Processing {source} ({species}) ---")

    # 1. Load cubes and set spectral axis to GHz

    cube_f = SpectralCube.read(feather_path).with_spectral_unit(u.GHz)

    cube_tp = SpectralCube.read(tp_path).with_spectral_unit(u.GHz)

    # 2. Get center and define aperture

    wcs_cel = cube_f.wcs.celestial

    ny, nx = cube_f.shape[1:]

    center_coord = wcs_cel.pixel_to_world((nx - 1) / 2., (ny - 1) / 2.)

    aperture = CircleSkyRegion(center=center_coord, radius=ap_radius_arcsec * u.arcsec)

    # 3. Extract the subcube FIRST (drastically reduces data size)

    subcube_f = cube_f.subcube_from_regions([aperture])

    subcube_tp = cube_tp.subcube_from_regions([aperture])

    # Enable huge ops just in case the subcube is still large

    subcube_f.allow_huge_operations = True

    subcube_tp.allow_huge_operations = True

    # 4. Convert only the small subcube to Kelvin

    subcube_f_k = subcube_f.to(u.K)

    subcube_tp_k = subcube_tp.to(u.K)

    # 5. Calculate spatial mean

    spec_f_k = subcube_f_k.mean(axis=(1, 2)).value

    spec_tp_k = subcube_tp_k.mean(axis=(1, 2)).value

    freq_f = cube_f.spectral_axis.value

    freq_tp = cube_tp.spectral_axis.value

    # Interpolate TP onto ALMA frequency grid for clean residuals

    sort_tp = np.argsort(freq_tp)

    spec_tp_k_interp = np.interp(freq_f, freq_tp[sort_tp], spec_tp_k[sort_tp])

    diff_k = spec_f_k - spec_tp_k_interp

    # --- Plotting ---

    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6), sharex=True, gridspec_kw={"height_ratios": [3, 1]})

    ax1.plot(freq_f, spec_tp_k_interp, label="Single Dish (TP)", color="#e41a1c", lw=2.0, linestyle="--")

    ax1.plot(freq_f, spec_f_k, label="Feathered (ALMA + TP)", color="#377eb8", lw=1.5, alpha=0.9)

    ax1.set_ylabel("Brightness Temperature [K]")

    ax1.set_title(f"Feathering Verification: {source} — {species} (Aperture $R = {ap_radius_arcsec}''$)")

    ax1.legend(frameon=True, loc="upper right")

    ax1.grid(True, linestyle=":", alpha=0.6)

    ax2.plot(freq_f, diff_k, color="black", lw=1.2, label="Residual (Feather - TP)")

    ax2.axhline(0, color="gray", linestyle="--", lw=0.8)

    ax2.set_xlabel("Frequency [GHz]")

    ax2.set_ylabel("Residual [K]")

    ax2.legend(frameon=True, loc="upper right")

    ax2.grid(True, linestyle=":", alpha=0.6)

    plt.tight_layout()

    out_png = os.path.join(output_dir, f"tp_vs_feather_{source}_{species}_SpectralCube.png")

    plt.savefig(out_png, dpi=300)

    plt.close()

    print(f"  [Success] Plot saved to: {out_png}")

if __name__ == "__main__":

    sources = ["G332.96", "G28", "G14"]

    species_list = ["CH3OH", "N2H+"]

    for src in sources:

        for sp in species_list:

            compare_tp_feather_spectralcube(src, sp, ap_radius_arcsec=40.0)

 No newline at end of file