Astrophysics Data Guide

from astropy.io import fits
import numpy as np

def inspect_fits(filepath: str) -> dict:
    """
    Inspect the structure of a FITS file.
    Returns information about each HDU (Header/Data Unit).
    """
    with fits.open(filepath) as hdul:
        info = []
        for i, hdu in enumerate(hdul):
            entry = {
                "index": i,
                "name": hdu.name,
                "type": type(hdu).__name__,
            }
            if hdu.data is not None:
                entry["shape"] = hdu.data.shape
                entry["dtype"] = str(hdu.data.dtype)
            if hasattr(hdu, "columns") and hdu.columns is not None:
                entry["columns"] = [c.name for c in hdu.columns]
            info.append(entry)
        return {"filename": filepath, "n_hdus": len(hdul), "hdus": info}

def read_fits_image(filepath: str, hdu_index: int = 0) -> tuple:
    """Read a FITS image and its WCS (World Coordinate System)."""
    from astropy.wcs import WCS

    with fits.open(filepath) as hdul:
        data = hdul[hdu_index].data
        header = hdul[hdu_index].header
        wcs = WCS(header)

    return data, wcs, header

from astropy.table import Table

def read_fits_catalog(filepath: str, hdu: int = 1) -> Table:
    """Read a FITS binary table extension as an Astropy Table."""
    catalog = Table.read(filepath, hdu=hdu)
    print(f"Catalog: {len(catalog)} objects, {len(catalog.columns)} columns")
    print(f"Columns: {catalog.colnames}")
    return catalog

from astropy.coordinates import SkyCoord, EarthLocation, AltAz
from astropy.time import Time
import astropy.units as u

def coordinate_transforms(ra_deg: float, dec_deg: float) -> dict:
    """
    Transform between astronomical coordinate systems.
    ra_deg, dec_deg: right ascension and declination in degrees (ICRS/J2000)
    """
    coord = SkyCoord(ra=ra_deg * u.degree, dec=dec_deg * u.degree, frame="icrs")

    return {
        "icrs": {
            "ra": coord.ra.to_string(unit=u.hourangle, precision=2),
            "dec": coord.dec.to_string(unit=u.degree, precision=2),
        },
        "galactic": {
            "l": round(coord.galactic.l.degree, 4),
            "b": round(coord.galactic.b.degree, 4),
        },
        "ecliptic": {
            "lon": round(coord.geocentricmeanecliptic.lon.degree, 4),
            "lat": round(coord.geocentricmeanecliptic.lat.degree, 4),
        },
    }

def compute_altaz(ra_deg: float, dec_deg: float,
                   obs_time: str, location: tuple) -> dict:
    """
    Compute altitude and azimuth for a target from a given location and time.
    location: (latitude_deg, longitude_deg, elevation_m)
    """
    target = SkyCoord(ra=ra_deg * u.degree, dec=dec_deg * u.degree)
    time = Time(obs_time)
    loc = EarthLocation(
        lat=location[0] * u.degree,
        lon=location[1] * u.degree,
        height=location[2] * u.m,
    )
    altaz_frame = AltAz(obstime=time, location=loc)
    altaz = target.transform_to(altaz_frame)

    return {
        "altitude_deg": round(altaz.alt.degree, 2),
        "azimuth_deg": round(altaz.az.degree, 2),
        "airmass": round(altaz.secz.value, 3) if altaz.alt.degree > 0 else None,
        "is_observable": altaz.alt.degree > 10,
    }

from photutils.aperture import CircularAperture, CircularAnnulus
from photutils.aperture import aperture_photometry

def perform_aperture_photometry(image: np.ndarray,
                                  positions: list[tuple],
                                  aperture_radius: float = 5.0,
                                  annulus_inner: float = 10.0,
                                  annulus_outer: float = 15.0) -> list[dict]:
    """
    Perform aperture photometry with local background subtraction.
    image: 2D numpy array (flux/counts)
    positions: list of (x, y) pixel coordinates of sources
    """
    apertures = CircularAperture(positions, r=aperture_radius)
    annuli = CircularAnnulus(positions, r_in=annulus_inner, r_out=annulus_outer)

    # Measure flux in aperture and annulus
    phot_table = aperture_photometry(image, [apertures, annuli])

    results = []
    for row in phot_table:
        # Background per pixel from annulus
        annulus_area = np.pi * (annulus_outer**2 - annulus_inner**2)
        bkg_per_pixel = row["aperture_sum_1"] / annulus_area

        # Background-subtracted flux
        aperture_area = np.pi * aperture_radius**2
        net_flux = row["aperture_sum_0"] - bkg_per_pixel * aperture_area

        # Instrumental magnitude
        if net_flux > 0:
            inst_mag = -2.5 * np.log10(net_flux)
        else:
            inst_mag = float("nan")

        results.append({
            "x": float(row["xcenter"]),
            "y": float(row["ycenter"]),
            "raw_flux": float(row["aperture_sum_0"]),
            "net_flux": round(float(net_flux), 2),
            "bkg_per_pixel": round(float(bkg_per_pixel), 2),
            "inst_mag": round(inst_mag, 4),
        })

    return results

from photutils.detection import DAOStarFinder
from astropy.stats import sigma_clipped_stats

def detect_sources(image: np.ndarray, fwhm: float = 3.0,
                    threshold_sigma: float = 5.0) -> Table:
    """
    Detect point sources in an astronomical image using DAOFind algorithm.
    """
    mean, median, std = sigma_clipped_stats(image, sigma=3.0)
    daofind = DAOStarFinder(fwhm=fwhm, threshold=threshold_sigma * std)
    sources = daofind(image - median)

    if sources is not None:
        sources.sort("flux", reverse=True)
        print(f"Detected {len(sources)} sources")
    return sources

from specutils import Spectrum1D, SpectralRegion
from specutils.analysis import line_flux, equivalent_width, centroid
import astropy.units as u

def analyze_spectrum(wavelength: np.ndarray,
                      flux: np.ndarray,
                      line_center: float,
                      line_width: float = 10.0) -> dict:
    """
    Analyze an emission or absorption line in a 1D spectrum.
    wavelength: array in Angstroms
    flux: array in erg/s/cm2/Angstrom
    line_center: expected line center in Angstroms
    line_width: width of spectral region to analyze
    """
    spectrum = Spectrum1D(
        spectral_axis=wavelength * u.Angstrom,
        flux=flux * u.Unit("erg / (s cm2 Angstrom)"),
    )

    region = SpectralRegion(
        (line_center - line_width) * u.Angstrom,
        (line_center + line_width) * u.Angstrom,
    )

    measured_flux = line_flux(spectrum, regions=region)
    ew = equivalent_width(spectrum, regions=region)
    center = centroid(spectrum, region)

    # Redshift from line center offset
    rest_wavelength = line_center  # assumed rest frame
    z = (center.value - rest_wavelength) / rest_wavelength

    return {
        "line_flux": f"{measured_flux:.4e}",
        "equivalent_width": f"{ew:.2f}",
        "measured_center_A": round(center.value, 2),
        "redshift": round(z, 6),
        "velocity_km_s": round(z * 299792.458, 1),
    }

from astroquery.vizier import Vizier
from astroquery.simbad import Simbad
from astroquery.sdss import SDSS

def query_simbad(object_name: str) -> dict:
    """Query SIMBAD for basic object information."""
    result = Simbad.query_object(object_name)
    if result is None:
        return {"found": False}

    return {
        "found": True,
        "name": object_name,
        "ra": str(result["RA"][0]),
        "dec": str(result["DEC"][0]),
        "object_type": str(result["OTYPE"][0]),
    }

def cone_search_vizier(ra_deg: float, dec_deg: float,
                        radius_arcmin: float = 1.0,
                        catalog: str = "II/246") -> Table:
    """
    Cone search in a VizieR catalog.
    Default catalog II/246 = 2MASS Point Source Catalog.
    """
    coord = SkyCoord(ra=ra_deg * u.degree, dec=dec_deg * u.degree)
    result = Vizier.query_region(
        coord, radius=radius_arcmin * u.arcmin, catalog=catalog
    )
    return result[0] if result else None