Projection

def project_polarizations_to_network(  # noqa: PLR0913
    polarizations: Mapping[str, GWpyTimeSeries],
    detector_names: Sequence[DetectorSpec],
    *,
    right_ascension: float,
    declination: float,
    polarization_angle: float,
    earth_rotation: bool = True,
) -> dict[str, GWpyTimeSeries]:
    """Project tensor plus/cross strains onto detectors using detector geometry.

    Built-in and custom detector codes are resolved through the LAL cached
    detector registry. Polarizations are interpolated in time with cubic
    splines (see user guide for caveats at edges).

    Args:
        polarizations: Mapping containing ``plus`` and ``cross`` GWpy time series
            on a common grid.
        detector_names: Sequence of IFO codes (e.g. ``H1``, ``L1``, ``V1``) or
            :class:`~gwmock_signal.detector.CustomDetector` instances, or a mix
            of both.
        right_ascension: Source right ascension in radians.
        declination: Source declination in radians.
        polarization_angle: Polarization angle psi in radians (tensor modes).
        earth_rotation: If ``True``, evaluate antenna patterns at time-dependent
            GPS times (recommended for longer signals). If ``False``, use a single
            reference time at the segment midpoint for patterns and delays.

    Returns:
            Mapping from each detector name to the projected strain as a GWpy time
            series (same length and sample rate as the inputs).

    Raises:
        TypeError: If ``polarizations`` is not a mapping of GWpy series as required.
        ValueError: If keys are missing, time grids disagree, or a detector name
            is not recognized.
    """
    hp, hc = _validate_polarizations(polarizations)
    normalized_names = [d if isinstance(d, str) else d.name for d in detector_names]
    if len(set(normalized_names)) != len(normalized_names):
        raise ValueError("detector_names must not contain duplicates.")
    detectors = _make_detectors(list(detector_names))

    time_array = cast(np.ndarray, hp.times.to_value())
    reference_time = float(0.5 * (time_array[0] + time_array[-1]))
    time_array_wrt_reference = time_array - reference_time

    minimum_number_of_data_points = 4
    interp_kind = "cubic" if len(time_array_wrt_reference) >= minimum_number_of_data_points else "linear"

    hp_func = interp1d(
        time_array_wrt_reference,
        hp.to_value(),
        kind=interp_kind,
        bounds_error=False,
        fill_value=0.0,
    )
    hc_func = interp1d(
        time_array_wrt_reference,
        hc.to_value(),
        kind=interp_kind,
        bounds_error=False,
        fill_value=0.0,
    )

    strains: dict[str, GWpyTimeSeries] = {}

    cosdec = np.cos(declination)
    sindec = np.sin(declination)
    cospsi = np.cos(polarization_angle)
    sinpsi = np.sin(polarization_angle)
    gmst_array = _gmst_accurate_array(time_array)
    gha_array = gmst_array - right_ascension
    cosgha = np.cos(gha_array)
    singha = np.sin(gha_array)

    for name, prefix in detectors:
        if earth_rotation:
            response, location = _reconstructed_geometry(prefix)

            # Vectorized time delay: time_delay = -location · prop_dir / c
            prop_dir = np.stack([cosdec * cosgha, -cosdec * singha, np.full(len(time_array), sindec)], axis=-1)
            time_delays = -np.dot(prop_dir, location) / constants.c.value

            shifted_times = time_array_wrt_reference - time_delays

            # Vectorized antenna pattern (same math as _antenna_pattern_lal, batched)
            gmst_antenna = _gmst_accurate_array(time_array + time_delays)
            gha_a = gmst_antenna - right_ascension
            cosgha_a = np.cos(gha_a)
            singha_a = np.sin(gha_a)

            # Shape (N, 3) — polarization basis vectors
            x_vec = np.stack(
                [
                    -cospsi * singha_a - sinpsi * cosgha_a * sindec,
                    -cospsi * cosgha_a + sinpsi * singha_a * sindec,
                    np.full(len(time_array), sinpsi * cosdec),
                ],
                axis=-1,
            )
            y_vec = np.stack(
                [
                    sinpsi * singha_a - cospsi * cosgha_a * sindec,
                    sinpsi * cosgha_a + cospsi * singha_a * sindec,
                    np.full(len(time_array), cospsi * cosdec),
                ],
                axis=-1,
            )

            # dx[n] = response @ x_vec[n], using row-vector form: x_vec @ response.T
            dx = x_vec @ response.T
            dy = y_vec @ response.T
            fp_vals = np.sum(x_vec * dx - y_vec * dy, axis=-1)
            fc_vals = np.sum(x_vec * dy + y_vec * dx, axis=-1)
        else:
            time_delay = _time_delay_from_earth_center_lal(
                prefix,
                right_ascension=right_ascension,
                declination=declination,
                t_gps=reference_time,
            )
            fp_vals, fc_vals = _antenna_pattern_lal(
                prefix,
                right_ascension=right_ascension,
                declination=declination,
                polarization_angle=polarization_angle,
                t_gps=reference_time,
            )
            shifted_times = time_array_wrt_reference - time_delay

        hp_shifted = hp_func(shifted_times)
        hc_shifted = hc_func(shifted_times)
        response = fp_vals * hp_shifted + fc_vals * hc_shifted

        strains[name] = GWpyTimeSeries(
            response,
            t0=float(time_array[0]),
            sample_rate=hp.sample_rate,
            name=name,
        )

    return strains