A Python script to calculate fiber aperture fraction to the total flux for Prime Focus Spectrograph (PFS) on Subaru Telescope

README.md

pfs_fiber_fraction

Compute the fiber aperture factor (encircled energy fraction) for PFS (Prime Focus Spectrograph). This is a Python port of the gsGeometricThroughput function from gsetc.c. See the documentation in the spt_ExposureTimeCalculator repository.

Dependencies

• numpy
• scipy

Usage

Basic execution

uv run pfs_fiber_fraction.py

Outputs the aperture factor for both an extended source and a point source using the default parameters (seeing=0.80", r_eff=0.30", wavelength=800nm, fieldang=0.45deg, decent=0.10").

Options

Option	Description	Default
--seeing	Seeing FWHM at 800nm [arcsec]	0.80
--wavelength	Wavelength [nm]	800.0
--r_eff	Effective radius of the source [arcsec] (exponential profile)	0.30
--decent	Fiber astrometric offset [arcsec]	0.10
--fieldang	Field angle [deg]	0.45
--mag	Total magnitude of the source (enables fiber magnitude output)	None
--config	Path to spectrograph config file	Built-in values

Examples

# Seeing 0.6" at wavelength 600nm
uv run pfs_fiber_fraction.py --seeing 0.6 --wavelength 600

# Compute fiber magnitude for a 22nd mag source
uv run pfs_fiber_fraction.py --mag 22.0

# Point source only (r_eff=0)
uv run pfs_fiber_fraction.py --r_eff 0.0

As a library

from pfs_fiber_fraction import default_spectrograph_config, geometric_throughput

spectro = default_spectrograph_config()
ee = geometric_throughput(spectro, seeing_fwhm_800=0.8, wavelength=800.0,
                          r_eff=0.0, decent=0.1, fieldang=0.45)

pfs_fiber_fraction.py

#!/usr/bin/env python3
# /// script
# dependencies = [
#   "scipy",
#   "numpy",
# ]
# ///

"""Compute the fiber aperture factor for PFS.

This script extracts the gsGeometricThroughput calculation from gsetc.c
and computes the fraction of light from an astronomical source that enters
a fiber, given the source size, seeing, and instrument parameters.
"""

import argparse
import math

import numpy as np
from scipy.integrate import quad
from scipy.special import j0, j1

ARCSEC_PER_URAD = 0.206264806247097  # 1 microradian in arcseconds
RAT_HL_SL_EXP = 1.67834  # half-light to scale radius ratio for exponential profile


def default_spectrograph_config():
    """Return default PFS spectrograph parameters.

    Values from config/PFS.20240714.dat (OPTICS, SPOT, FIBER lines).
    To update, see the corresponding lines in the config file:
      OPTICS <D_outer> <centobs> <rfov> <EFL[0..4]>
      SPOT <rms_spot[0..4]>
      FIBER <fiber_ent_rad>
    """
    return {
        # Telescope outer diameter [m]
        "D_outer": 8.3,
        # Central obscuration ratio
        "centobs": 0.22,
        # Field of view radius [deg]
        "rfov": 0.65,
        # Effective focal length [m] (center to edge)
        "EFL": [
            23.39,
            23.92,
            24.48,
            25.07,
            25.68,
        ],
        # RMS spot size per axis [um] (center to edge)
        "rms_spot": [
            7.1,
            8.4,
            9.7,
            11.2,
            13.2,
        ],
        # Fiber entrance radius [um]
        "fiber_ent_rad": 63.5,
    }


def read_spectrograph_config(filename):
    """Read spectrograph config file and return parameters needed for throughput calculation.

    Returns a dict with keys: D_outer, centobs, rfov, EFL, rms_spot, fiber_ent_rad.
    """
    params = {}
    with open(filename) as f:
        for line in f:
            line = line.strip()
            if not line or line.startswith("#"):
                continue
            if line.startswith("OPTICS"):
                parts = line.split()[1:]
                params["D_outer"] = float(parts[0])
                params["centobs"] = float(parts[1])
                params["rfov"] = float(parts[2])
                params["EFL"] = [float(x) for x in parts[3:8]]
            elif line.startswith("SPOT"):
                parts = line.split()[1:]
                params["rms_spot"] = [float(x) for x in parts[:5]]
            elif line.startswith("FIBER"):
                params["fiber_ent_rad"] = float(line.split()[1])
    return params


def geometric_throughput(spectro, seeing_fwhm_800, wavelength, r_eff, decent, fieldang):
    """Compute the fiber aperture factor (encircled energy fraction).

    This is a direct port of gsGeometricThroughput from gsetc.c.

    Parameters
    ----------
    spectro : dict
        Spectrograph parameters from read_spectrograph_config.
    seeing_fwhm_800 : float
        Seeing FWHM at 800 nm in arcsec.
    wavelength : float
        Wavelength in nm.
    r_eff : float
        Effective radius of the source in arcsec (0 for point source).
    decent : float
        Fiber astrometric offset in arcsec.
    fieldang : float
        Field angle in degrees.

    Returns
    -------
    float
        Fraction of light entering the fiber (0 to 1).
    """
    # Interpolate RMS spot size and EFL at the given field angle
    rfov = spectro["rfov"]
    # 5 equally spaced sample points from center (0) to edge (rfov)
    xp = [i * rfov / 4 for i in range(5)]
    efl = float(np.interp(fieldang, xp, spectro["EFL"]))
    sigma = float(np.interp(fieldang, xp, spectro["rms_spot"]))

    # Convert RMS spot size from microns to arcsec
    sigma *= ARCSEC_PER_URAD / efl

    # Fiber radius in arcsec
    R = spectro["fiber_ent_rad"] / efl * ARCSEC_PER_URAD

    # Seeing MTF scale
    uscale = 0.465 / seeing_fwhm_800 * (wavelength / 800.0) ** 0.2

    # Galaxy scale length
    rs = r_eff / RAT_HL_SL_EXP

    # Diffraction angle
    D_outer = spectro["D_outer"]
    centobs = spectro["centobs"]
    theta_D = 0.001 * wavelength / D_outer / (1.0 - centobs) * ARCSEC_PER_URAD

    def integrand(u):
        k = 2.0 * math.pi * u
        # Telescope PSF in Fourier space (Kolmogorov seeing + Gaussian spot + diffraction)
        Gtilde = (
            math.exp(-k * k * sigma * sigma / 2.0 - (u / uscale) ** (5.0 / 3.0))
            * j0(k * decent)
            * math.exp(-4.0 / math.pi * u * theta_D)
        )
        # Galaxy profile (Fourier transform of exponential disk)
        ftilde = (1.0 + (k * rs) ** 2) ** (-1.5)
        return R * 2.0 * math.pi * j1(2.0 * math.pi * u * R) * ftilde * Gtilde

    EE, _ = quad(integrand, 0, np.inf)

    return float(EE)


def main():
    parser = argparse.ArgumentParser(
        description="Compute fiber aperture factor for PFS"
    )
    parser.add_argument(
        "--config",
        default=None,
        help="Path to spectrograph config file (default: use built-in PFS parameters)",
    )
    parser.add_argument(
        "--r_eff",
        type=float,
        default=0.3,
        help="Effective radius of source in arcsec (default: 0.3)",
    )
    parser.add_argument(
        "--seeing",
        type=float,
        default=0.80,
        help="Seeing FWHM at 800nm in arcsec (default: 0.80)",
    )
    parser.add_argument(
        "--fieldang",
        type=float,
        default=0.45,
        help="Field angle in degrees (default: 0.45)",
    )
    parser.add_argument(
        "--decent",
        type=float,
        default=0.10,
        help="Fiber astrometric offset in arcsec (default: 0.10)",
    )
    parser.add_argument(
        "--wavelength",
        type=float,
        default=800.0,
        help="Wavelength in nm (default: 800.0)",
    )
    parser.add_argument(
        "--mag",
        type=float,
        default=None,
        help="Total magnitude of the source (optional)",
    )
    args = parser.parse_args()

    if args.config is not None:
        spectro = read_spectrograph_config(args.config)
    else:
        spectro = default_spectrograph_config()

    ee_ext = geometric_throughput(
        spectro, args.seeing, args.wavelength, args.r_eff, args.decent, args.fieldang
    )
    ee_pt = geometric_throughput(
        spectro, args.seeing, args.wavelength, 0.0, args.decent, args.fieldang
    )

    print(
        f'Fiber aperture factor [@{args.wavelength:.0f}nm, r_eff={args.r_eff:.2f}"(exp)] = {ee_ext:10.8f}'
    )
    print(
        f"Fiber aperture factor [@{args.wavelength:.0f}nm,     point source] = {ee_pt:10.8f}"
    )

    if args.mag is not None:
        mag_fiber_ext = args.mag - 2.5 * math.log10(ee_ext)
        mag_fiber_pt = args.mag - 2.5 * math.log10(ee_pt)
        print(f"\nTotal magnitude = {args.mag:.4f}")
        print(f'Fiber magnitude [r_eff={args.r_eff:.2f}"(exp)] = {mag_fiber_ext:.4f}')
        print(f"Fiber magnitude [    point source] = {mag_fiber_pt:.4f}")


if __name__ == "__main__":
    main()