from __future__ import annotations

import numbers
from typing import Any, Dict, Final, Literal, Union, cast, get_args

import lal
import lalsimulation as lalsim
import numpy as np

from .eob.hamiltonian.Ham_align_a6_apm_AP15_DP23_gaugeL_Tay_C import (
    Ham_align_a6_apm_AP15_DP23_gaugeL_Tay_C as Ham_aligned_opt,
)
from .eob.hamiltonian.Ham_AvgS2precess_simple_cython_PA_AD import (
    Ham_AvgS2precess_simple_cython_PA_AD as Ham_prec_pa_cy,
)
from .eob.waveform.waveform import SEOBNRv5RRForce
from .eob.waveform.waveform_ecc import SEOBNRv5RRForceEcc
from .models import SEOBNRv5EHM, SEOBNRv5HM
from .models.model import Model

#: Supported approximants
SupportedApproximants = Literal["SEOBNRv5HM", "SEOBNRv5PHM", "SEOBNRv5EHM"]


def generate_prec_hpc_opt(
    q: float,
    chi1: np.ndarray,
    chi2: np.ndarray,
    omega_start: float,
    omega_ref: float | None = None,
    settings: Dict[Any, Any] | None = None,
    debug: bool = False,
) -> tuple[np.array, np.array] | tuple[np.array, np.array, Model]:
    """
    Generate the GW wave polarizations for precessing model in an optimised way.
    In particular, do not compute the inertial frame modes, instead write the
    polarizations directly summing over the co-precessing frame modes.

    Args:
        q (float): Mass ratio >=1
        chi1 (np.ndarray): Dimensionless spin of the primary
        chi2 (np.ndarray): Dimensionless spin of the secondary
        omega_start (float): Starting orbital frequency, in geometric units
        omega_ref (float, optional): Reference orbital frequency, in geometric
            units. Defaults to None, in which case equals omega_start
        settings (Dict[Any,Any], optional): Any additional settings to pass to
            model. Defaults to None
        debug (bool, optional): Run in debug mode and return the model object.
            Defaults to False

    Returns:
        Tuple: Either (time, polarization) or (time, polarizations, model)
            if debug=True
    """

    if settings is None:
        settings = {"polarizations_from_coprec": True}
    else:
        settings.update(polarizations_from_coprec=True)

    if q < 1.0:
        raise ValueError("mass-ratio has to be positive and with convention q>=1")
    if omega_start < 0:
        raise ValueError("omega_start has to be positive")

    if np.linalg.norm(chi1) > 1 or np.linalg.norm(chi2) > 1:
        raise ValueError("chi1 and chi2 have to respect Kerr limit (|chi|<=1)")

    RR_f = SEOBNRv5RRForce()
    model = SEOBNRv5HM.SEOBNRv5PHM_opt(
        q,
        *chi1,
        *chi2,
        omega_start,
        Ham_prec_pa_cy,
        RR_f,
        omega_ref=omega_ref,
        settings=settings,
    )

    model()
    if debug:
        return model.t, model.hpc, model
    else:
        return model.t, model.hpc


def _check_spins(
    chi1: Union[
        float,
        np.ndarray,
        list,
        tuple[
            float,
            float,
            float,
        ],
    ],
    chi2: Union[
        float,
        np.ndarray,
        list,
        tuple[
            float,
            float,
            float,
        ],
    ],
    approximant: SupportedApproximants,
) -> tuple[float, float, float, float, float, float]:
    """
    Perform some checks on the dimensionless spin components of the black holes.

    Args:
        chi1 (Union[float,np.ndarray]): Dimensionless spin of primary.
            If ``float``, interpreted as z component
        chi2 (Union[float,np.ndarray]): Dimensionless spin of secondary.
            If ``float``, interpreted as z component

    Returns:
        Tuple: Values of the dimensionless spin components of the primary and
            secondary black holes
    """

    chi1_ = np.asarray(chi1, dtype=np.float64)
    chi2_ = np.asarray(chi2, dtype=np.float64)

    if len(chi1_.shape) == 0:
        if isinstance(chi1, bool) or not isinstance(
            chi1, numbers.Real
        ):  # chi1 is correct here
            raise ValueError("Boolean spin values unsupported")
        chi1_ = np.array([0, 0, chi1_], dtype=np.float64)
    elif chi1_.shape[0] > 3:
        raise ValueError("Incorrect number of spin elements.")
    elif any(
        isinstance(_, bool) or not isinstance(_, numbers.Real) for _ in chi1
    ):  # chi1 is correct here
        raise ValueError("Boolean spin values unsupported")

    if len(chi2_.shape) == 0:
        if isinstance(chi2, bool) or not isinstance(
            chi2, numbers.Real
        ):  # chi2 is correct here
            raise ValueError("Boolean spin values unsupported")
        chi2_ = np.array([0, 0, chi2_], dtype=np.float64)
    elif chi2_.shape[0] > 3:
        raise ValueError("Incorrect number of spin elements.")
    elif any(
        isinstance(_, bool) or not isinstance(_, numbers.Real) for _ in chi2
    ):  # chi2 is correct here
        raise ValueError("Boolean spin values unsupported")

    if approximant in ["SEOBNRv5HM", "SEOBNRv5EHM"] and (
        (np.abs(chi1_[:2]).max() > 1e-10) or (np.abs(chi2_[:2]).max() > 1e-10)
    ):
        raise ValueError(
            "In-plane spin components must be zero for calling the "
            "non-precessing approximant."
        )

    if np.linalg.norm(chi1_) > 1 or np.linalg.norm(chi2_) > 1:
        raise ValueError("Dimensionless spin magnitudes cannot be greater than 1!")

    return cast(tuple[float, float, float], tuple(float(_) for _ in chi1_)) + cast(
        tuple[float, float, float], tuple(float(_) for _ in chi2_)
    )


def generate_modes_opt(
    q: float,
    chi1: Union[float, np.ndarray],
    chi2: Union[float, np.ndarray],
    omega_start: float,
    omega_ref: float | None = None,
    eccentricity: float | None = None,
    rel_anomaly: float | None = None,
    approximant: SupportedApproximants = "SEOBNRv5HM",
    settings: Dict[Any, Any] = None,
    debug: bool = False,
) -> tuple[np.array, np.array] | tuple[np.array, np.array, Model]:
    """
    Compute the GW waveform modes for the given configuration and approximant.

    Args:
        q (float): Mass ratio >=1
        chi1 (Union[float,np.ndarray]): Dimensionless spin of primary.
            If ``float``, interpreted as z component
        chi2 (Union[float,np.ndarray]): Dimensionless spin of secondary.
            If ``float``, interpreted as z component
        omega_start (float): Starting orbital frequency, in geometric units
        omega_ref (float, optional): Reference orbital frequency, in
            geometric units. Defaults to None, in which case equals omega_start.
        eccentricity (float, optional): Eccentricity of the orbit,
            corresponding to the Keplerian (relativistic) parametrization.
            Defaults to 0 (no eccentricity)
        rel_anomaly (float, optional): Relativistic anomaly. Radial phase
            which parametrizes the orbit within the Keplerian (relativistic)
            parametrization. Defaults to 0 (periastron)
        approximant (SupportedApproximants, optional): The approximant to use.
            Defaults to "SEOBNRv5HM"
        settings (Dict[Any,Any], optional): Additional settings to pass to
            model. Defaults to ``None``
        debug (bool, optional): Run in debug mode. Defaults to False

    Raises:
        ValueError: If input parameters are not physical
        NotImplementedError: If the approximant requested does not exist

    Returns:
        Tuple: Either (time, dictionary of modes) or (time, dictionary of modes, model)
        if ``debug=True``
    """

    if q < 1.0:
        raise ValueError("mass-ratio has to be positive and with convention q>=1")
    if omega_start < 0:
        raise ValueError("omega_start has to be positive")

    if approximant != "SEOBNRv5EHM" and (
        (eccentricity is not None and eccentricity != 0)
        or (rel_anomaly is not None and rel_anomaly != 0)
    ):
        raise ValueError(
            "The selected approximant does not support input values for the "
            "eccentricity and relativistic anomaly."
        )

    chi1_x, chi1_y, chi1_z, chi2_x, chi2_y, chi2_z = _check_spins(
        chi1, chi2, approximant=approximant
    )

    # checks the anti-symmetric parameters
    if settings and settings.get("enable_antisymmetric_modes", False):
        if approximant != "SEOBNRv5PHM":
            raise ValueError(
                f"Antisymmetric modes not available for approximant {approximant}."
            )

    if approximant == "SEOBNRv5HM":
        RR_f = SEOBNRv5RRForce()
        model = SEOBNRv5HM.SEOBNRv5HM_opt(
            q,
            chi_1=chi1_z,
            chi_2=chi2_z,
            omega0=omega_start,
            H=Ham_aligned_opt,
            RR=RR_f,
            settings=settings,
        )
        model()

    elif approximant == "SEOBNRv5PHM":
        RR_f = SEOBNRv5RRForce()
        model = SEOBNRv5HM.SEOBNRv5PHM_opt(
            q,
            chi1_x=chi1_x,
            chi1_y=chi1_y,
            chi1_z=chi1_z,
            chi2_x=chi2_x,
            chi2_y=chi2_y,
            chi2_z=chi2_z,
            omega_start=omega_start,
            H=Ham_prec_pa_cy,
            RR=RR_f,
            omega_ref=omega_ref,
            settings=settings,
        )
        model()

    elif approximant == "SEOBNRv5EHM":
        if eccentricity is None:
            eccentricity = 0.0
        if rel_anomaly is None:
            rel_anomaly = 0.0
        if not 0.0 <= eccentricity < 1.0:
            raise ValueError(
                "The value of eccentricity must be inside the interval [0, 1)."
            )
        RR_force = SEOBNRv5RRForceEcc("Ecc")
        model = SEOBNRv5EHM.SEOBNRv5EHM_opt(
            q=q,
            chi_1=chi1_z,
            chi_2=chi2_z,
            omega_start=omega_start,
            eccentricity=eccentricity,
            rel_anomaly=rel_anomaly,
            H=Ham_aligned_opt,
            RR=RR_force,
            settings=settings,
        )
        model()

    else:
        raise NotImplementedError(f"Approximant '{approximant}' is not available")

    if debug:
        return model.t, model.waveform_modes, model
    else:
        return model.t, model.waveform_modes


class GenerateWaveform:
    """
    Class for generating modes, time-domain polarizations and frequency-domain
    polarizations following LAL conventions.
    """

    def __init__(self, parameters):
        """
        Initialize the class with the given ``parameters``.

        ``parameters`` is a dictionary whose keys are defined as follows:

        Parameters
        ----------
        float mass1:
            Mass of companion 1, in solar masses - Required
        float mass2:
            Mass of companion 2, in solar masses - Required
        float spin1x:
            x-component of dimensionless spin of companion 1 - Default: 0
        float spin1y:
            y-component of dimensionless spin of companion 1 - Default: 0
        float spin1z:
            z-component of dimensionless spin of companion 1 - Default: 0
        float spin2x:
            x-component of dimensionless spin of companion 2 - Default: 0
        float spin2y:
            y-component of dimensionless spin of companion 2 - Default: 0
        float spin2z:
            z-component of dimensionless spin of companion 2 - Default: 0
        float eccentricity:
            Keplerian eccentricity of the orbit - Default: 0
        float rel_anomaly:
            Relativistic anomaly - Default: 0
        float distance:
            Distance to the source, in Mpc - Default: 100 Mpc
        float inclination:
            Inclination of the source, in radians - Default: 0
        float phi_ref:
            Orbital phase at the reference frequency, in radians - Default: 0
        float f22_start:
            Starting waveform generation frequency, in Hz - Default: 20 Hz
        float f_ref:
            The reference frequency, in Hz - Default: ``f22_start``
        float deltaT:
            Time spacing, in seconds - Default: 1/2048 seconds
        float f_max:
            Maximum frequency, in Hz - Default: 1024 Hz
        float deltaF:
            Frequency spacing, in Hz - Default: 0.125
        list ModeArray:
        list mode_array:
            Mode content (only positive must be specified, e.g ``[(2,2),(2,1)]``).
            Defaults to ``None`` (all modes, see notes below).

        dict domega_dict:
            The non-GR fractional deviation to the frequencies for each mode.
        dict dA_dict:
            The non-GR fractional deviation to the merger amplitudes for each mode.
        dict dw_dict:
            The non-GR fractional deviation to the merger frequencies.
        dict dtau_dict:
            The non-GR fractional deviation to the damping times for each mode.
            Values should be :math:`> -1`.
        float dTpeak:
            The non-GR additive deviation to the amplitude's peak time.
        float da6:
            The non-GR additive deviation to the ``a6`` calibration parameter.
        float ddSO:
            The non-GR additive deviation to the dSO calibration parameter.
        bool deltaT_sampling:
            If set to ``True``, throws an error if the the attachment time
            induced by negative values of the deviation ``dTpeak`` is beyond the last
            point calculated from the dynamics. In those cases, the attachment time is
            set to be the last point of the dynamics.
            Setting the parameter to ``True`` prevents having incorrect posteriors when
            sampling over ``dTpeak``, as those parameters would be rejected.
            Needs to be set to ``True`` if ``dTpeak`` is non-zero. Defaults to ``False``.
        bool omega_prec_deviation:
            If ``True`` (default), the fractional deviations to the J-frame QNM frequencies
            are included into the precession rate computation (Eq. 13 in
            `arXiv:2301.06558 <https://arxiv.org/abs/2301.06558>`_ ).

        str initial_conditions:
            Possible values are ``adiabatic`` (default) and ``postadiabatic``. Used
            only for ``approximant="SEOBNRv5PHM"`` and ``postadiabatic`` is ``False``.
        str initial_conditions_postadiabatic_type:
            Possible values are ``analytic`` (default) and ``numeric``. Used together
            with ``initial_conditions``.
        bool postadiabatic:
            Defaults to ``True``.
        str postadiabatic_type:
            Either ``analytic`` (default) or ``numeric``. Used only when ``postadiabatic``
            is ``True``.
        float tol_PA:
            Tolerance for the root finding routine in case ``postadiabatic_type="numeric"``
            is used. Defaults to 1e-11.

        int lmax:
            If specified, will add all the modes compatible with the approximant up
            to :math:`\\ell =` ``lmax``. If provided together with ``mode_array``, the
            two parameters should be set consistently. See the notes for details.
        bool enable_antisymmetric_modes:
            when ``True`` will enable the handling of the anti-symmetric modes. Available
            only for the precessing approximant ``SEOBNRv5PHM``.
        list[tuple[int, int]] antisymmetric_modes:
            Select the modes for which the anti-symmetric processing is active.
        bool antisymmetric_modes_hm:
            if set to ``True``, will enable the higher modes ``(3,3)`` and ``(4,4)`` for the
            anti-symmetric modes corrections. Overrides ``antisymmetric_modes``.

        str approximant:
            * ``SEOBNRv5HM`` (default)
            * ``SEOBNRv5PHM``
            * ``SEOBNRv5EHM``

        float rtol_ode:
            Relative tolerance of the ODE integrators. Defaults to 1e-11.
        float atol_ode:
            Absolute tolerance of the ODE integrators. Defaults to 1e-12.

        bool gwsignal_environment:
            if set and evaluates to ``True``, indicates that the code is being run
            through ``GWSignal`` or environment requiring compliance with
            LVK reviews. Defaults to ``False``.

        Note
        ----
        The default modes are ``(2, 2)``, ``(2, 1)``, ``(3, 3)``, ``(3, 2)``,
        ``(4, 4)`` and ``(4, 3)``. In particular ``(5, 5)`` is not included
        and should be explicitly set through ``ModeArray``.

        All GR deviations default to 0.
        For the dictionaries ``domega_dict``, ``dA_dict``, ``dw_dict``, ``dtau_dict``,
        keys are the modes ``l,m`` as a string, for :math:`\\ell > 0`.
        For a description of the non-GR deviation parameters, see Section II of [Pompili2025]_.

        If ``lmax`` is set, it will activate all the supported modes for the selected
        approximant up until :math:`\\ell =` ``lmax``. If specified together with ``mode_array``,
        then the selected modes from ``lmax`` should be the same as the ones provides in
        ``mode_array``, otherwise a ``ValueError`` exception is raised. For instance,
        ``lmax=3`` and ``mode_array=[(2, 2), (3, 2)]`` will create inconsistencies for the
        selection of the modes ``(2, 1)`` and ``(3, 3)``.
        """

        self.swap_masses: bool = False
        parameters = self._validate_parameters(parameters)
        parameters = self._validate_eccentricity_parameters(parameters=parameters)
        self.parameters: dict[str, Any] | None = parameters

    @property
    def model(self):
        if not hasattr(self, "_model"):
            raise ValueError("A model object has not been created!")
        return self._model

    def _validate_parameters(self, parameters):
        if "mass1" not in parameters:
            raise ValueError("mass1 has to be specified!")

        if "mass2" not in parameters:
            raise ValueError("mass2 has to be specified!")

        if (
            not isinstance(parameters["mass1"], numbers.Real)
            or isinstance(parameters["mass1"], bool)
            or not isinstance(parameters["mass2"], numbers.Real)
            or isinstance(parameters["mass2"], bool)
        ):
            raise ValueError("Only 'float' and 'int' values of masses are supported.")

        if parameters["mass1"] > 0 and parameters["mass2"] > 0:
            mass1 = parameters["mass1"]
            mass2 = parameters["mass2"]
        else:
            raise ValueError("Masses have to be positive!")

        Mtot = mass1 + mass2
        if Mtot < 0.001 or Mtot > 1e12:
            raise ValueError("Unreasonable value for total mass, aborting.")

        if mass1 * mass2 / Mtot**2 < 100.0 / (1 + 100) ** 2:
            raise ValueError(
                "Internal function call failed: Input domain error. Model is "
                "only valid for systems with mass-ratio up to 100."
            )

        default_params: Final = {
            "spin1x": 0.0,
            "spin1y": 0.0,
            "spin1z": 0.0,
            "spin2x": 0.0,
            "spin2y": 0.0,
            "spin2z": 0.0,
            "eccentricity": 0.0,
            "rel_anomaly": 0.0,
            "distance": 100.0,
            "inclination": 0.0,
            "phi_ref": 0.0,
            "f22_start": 20.0,
            "deltaT": 1.0 / 2048.0,
            "deltaF": 0.0,
            "ModeArray": None,
            "mode_array": None,
            "approximant": "SEOBNRv5HM",
            "conditioning": 2,
            "polarizations_from_coprec": True,
            "initial_conditions": "adiabatic",
            "initial_conditions_postadiabatic_type": "analytic",
            "postadiabatic": True,
            "postadiabatic_type": "analytic",
            "r_size_input": 12,
            "dA_dict": {
                "2,2": 0.0,
                "2,1": 0.0,
                "3,3": 0.0,
                "3,2": 0.0,
                "4,4": 0.0,
                "4,3": 0.0,
                "5,5": 0.0,
            },
            "dw_dict": {
                "2,2": 0.0,
                "2,1": 0.0,
                "3,3": 0.0,
                "3,2": 0.0,
                "4,4": 0.0,
                "4,3": 0.0,
                "5,5": 0.0,
            },
            "dTpeak": 0.0,
            "da6": 0.0,
            "ddSO": 0.0,
            "domega_dict": {
                "2,2": 0.0,
                "2,1": 0.0,
                "3,3": 0.0,
                "3,2": 0.0,
                "4,4": 0.0,
                "4,3": 0.0,
                "5,5": 0.0,
            },
            "dtau_dict": {
                "2,2": 0.0,
                "2,1": 0.0,
                "3,3": 0.0,
                "3,2": 0.0,
                "4,4": 0.0,
                "4,3": 0.0,
                "5,5": 0.0,
            },
            "tol_PA": 1e-11,
            "rtol_ode": 1e-11,
            "atol_ode": 1e-12,
            "deltaT_sampling": False,
            "omega_prec_deviation": True,
            "enable_antisymmetric_modes": False,
            "antisymmetric_modes": [(2, 2)],
            "antisymmetric_modes_hm": False,
            "ivs_mrd": None,
            "antisymmetric_fits_version": 240417,
            "gwsignal_environment": False,
        }

        if "approximant" in parameters and parameters["approximant"] == "SEOBNRv5EHM":
            if "postadiabatic" not in parameters:
                parameters["postadiabatic"] = False
            if "conditioning" not in parameters:
                parameters["conditioning"] = 1

        # Fills the provided parameters over the default ones
        parameters = default_params | parameters

        if "f_ref" not in parameters.keys():
            parameters["f_ref"] = parameters["f22_start"]

        if parameters["approximant"] not in get_args(SupportedApproximants):
            raise ValueError("Approximant not implemented!")

        # Disable direct polarizations for aligned-spin model
        if parameters["approximant"] in ["SEOBNRv5HM", "SEOBNRv5EHM"]:
            parameters["polarizations_from_coprec"] = False

        (
            parameters["spin1x"],
            parameters["spin1y"],
            parameters["spin1z"],
            parameters["spin2x"],
            parameters["spin2y"],
            parameters["spin2z"],
        ) = _check_spins(
            chi1=(
                parameters["spin1x"],
                parameters["spin1y"],
                parameters["spin1z"],
            ),
            chi2=(
                parameters["spin2x"],
                parameters["spin2y"],
                parameters["spin2z"],
            ),
            approximant=parameters["approximant"],
        )

        for param in [
            "distance",
            "inclination",
            "phi_ref",
            "f22_start",
            "f_ref",
            "deltaT",
            "f_max",
            "f_min",
            "deltaF",
            "rel_anomaly",
            "eccentricity",
        ]:
            if param in parameters and (
                not isinstance(parameters[param], numbers.Real)
                or isinstance(parameters[param], bool)
            ):
                raise ValueError(f"{param} has to be a real number!")

        # Only integer parameters (excluding booleans)
        for param in ["conditioning"]:
            if param in parameters and (
                not isinstance(parameters[param], numbers.Integral)
                or isinstance(parameters[param], bool)
            ):
                raise ValueError(f"{param} has to be an integer number!")

        # Only boolean parameters
        for param in ["postadiabatic"]:
            if param in parameters and (not isinstance(parameters[param], bool)):
                raise ValueError(f"{param} has to be a boolean!")

        if parameters["initial_conditions"] not in ["adiabatic", "postadiabatic"]:
            raise ValueError("Unrecognised setting for initial conditions.")

        if parameters["initial_conditions_postadiabatic_type"] not in [
            "numeric",
            "analytic",
        ]:
            raise ValueError(
                "Unrecognised setting for initial conditions postadiabatic type."
            )

        if parameters["postadiabatic_type"] not in ["numeric", "analytic"]:
            raise ValueError("Unrecognised setting for dynamics postadiabatic type.")

        if parameters["ModeArray"] is not None and parameters["mode_array"] is not None:
            raise ValueError(
                "Only one setting can be specified between ModeArray and mode_array."
            )

        # Additional default value can be specified after the check
        # (would divide by 0 if deltaT is "False")
        if "f_max" not in parameters.keys():
            parameters["f_max"] = 0.5 / parameters["deltaT"]

        # Check frequencies
        for param in ["f22_start", "f_ref", "f_max", "deltaT", "distance"]:
            if parameters[param] <= 0:
                raise ValueError(f"{param} has to be positive!")

        # deltaF = 0 indicates that it has to be estimated from the
        # intrinsic parameters following lalsimulation
        if parameters["deltaF"] < 0:
            raise ValueError("deltaF has to be positive!")

        if (
            parameters["f_max"] < parameters["f22_start"]
            or parameters["f_max"] < parameters["f_ref"]
        ):
            raise ValueError("'f_max' cannot be smaller than 'f22_start' or 'f_ref'!")

        if mass2 > mass1:
            self.swap_masses = True
            aux_spin1 = [
                parameters["spin1x"],
                parameters["spin1y"],
                parameters["spin1z"],
            ]
            parameters["spin1x"], parameters["spin1y"], parameters["spin1z"] = [
                -parameters["spin2x"],
                -parameters["spin2y"],
                parameters["spin2z"],
            ]
            parameters["spin2x"], parameters["spin2y"], parameters["spin2z"] = [
                -aux_spin1[0],
                -aux_spin1[1],
                aux_spin1[2],
            ]
            aux_mass1 = parameters["mass1"]
            parameters["mass1"] = parameters["mass2"]
            parameters["mass2"] = aux_mass1

        if parameters["initial_conditions"] not in ["adiabatic", "postadiabatic"]:
            raise ValueError("Unrecognised setting for initial conditions.")

        if parameters["initial_conditions_postadiabatic_type"] not in [
            "numeric",
            "analytic",
        ]:
            raise ValueError(
                "Unrecognised setting for initial conditions postadiabatic type."
            )

        if parameters["postadiabatic_type"] not in ["numeric", "analytic"]:
            raise ValueError("Unrecognised setting for dynamics postadiabatic type.")

        if parameters["ModeArray"] is not None and parameters["mode_array"] is not None:
            raise ValueError(
                "Only one setting can be specified between ModeArray and mode_array."
            )

        if parameters.get("enable_antisymmetric_modes", False):
            if parameters["approximant"] != "SEOBNRv5PHM":
                raise ValueError(
                    "Antisymmetric modes not available for approximant "
                    f"{parameters['approximant']}."
                )

            # left for posterity
            # if parameters.get("gwsignal_environment", False):
            #     warnings.warn(
            #         "This code is currently UNREVIEWED, use with caution!!", UserWarning
            #     )

        return parameters

    @staticmethod
    def _validate_eccentricity_parameters(parameters):
        if parameters["approximant"] == "SEOBNRv5EHM":
            # - EccIC: Parameter determining the type of initial frequency.
            # EccIC = 0 for instantaneous initial orbital frequency,
            # and EccIC = 1 for orbit-averaged initial orbital frequency.
            # - t_backwards: The amount of time (in geometric units) to
            # step back in time a given binary system
            # - secular_bwd_int: If True, performs a backwards evolution
            # of a set of secular evolution equations in case the starting
            # separation is smaller than r = 10 M. If False, then raise
            # an error if the starting separation is smaller than r = 10 M.
            # - warning_bwd_int: If True, sends a warning to the user if
            # the backwards integration of the full EOMs is activated
            # (this happens if t_backwards != 0)
            # - warning_secular_bwd_int: If True, sends a warning to the
            # user if the secular backwards integration was activated
            default_eccentricity_parameters = {
                "EccIC": 1,
                "t_backwards": 0.0,
                "secular_bwd_int": True,
                "warning_bwd_int": True,
                "warning_secular_bwd_int": True,
            }
            parameters = default_eccentricity_parameters | parameters

            # Checks already done in the _validate_parameters. We keep
            # here only the consistency checks on the values to make
            # the error messages coherent
            assert not (
                isinstance(parameters["rel_anomaly"], bool)
                or isinstance(parameters["eccentricity"], bool)
            )

            if not 0.0 <= parameters["eccentricity"] < 1.0:
                raise ValueError("Eccentricity must be within the interval [0, 1).")

            if parameters["postadiabatic"]:
                raise ValueError(
                    "The approximant 'SEOBNRv5EHM' does not support a postadiabatic "
                    "evolution. Please, set 'postadiabatic':False in the input "
                    "parameters."
                )

            if parameters["conditioning"] != 1:
                raise ValueError(
                    "The approximant 'SEOBNRv5EHM' cannot be used with the "
                    f"setting 'conditioning':{parameters['conditioning']}. "
                    "Please, set 'conditioning':1 in the input parameters."
                )

            if parameters["f_ref"] != parameters["f22_start"]:
                raise ValueError(
                    "The approximant 'SEOBNRv5EHM' does not support the choice "
                    "for a reference frequency. Please, set 'f_ref' = 'f22_start'."
                )

        else:
            if parameters["eccentricity"] != 0.0 or parameters["rel_anomaly"] != 0.0:
                raise ValueError(
                    "Use the approximant 'SEOBNRv5EHM' for a system with "
                    "non-zero eccentricity or relativistic anomaly."
                )

        return parameters

    def generate_td_modes(self):
        """
        Generate dictionary of positive and negative m modes in physical units.
        """

        fmin, dt = self.parameters["f22_start"], self.parameters["deltaT"]
        f_ref = self.parameters.get("f_ref")
        m1, m2 = self.parameters["mass1"], self.parameters["mass2"]
        Mtot = m1 + m2
        eccentricity = self.parameters["eccentricity"]
        rel_anomaly = self.parameters["rel_anomaly"]
        dist = self.parameters["distance"]
        approx: SupportedApproximants = cast(
            SupportedApproximants, self.parameters["approximant"]
        )

        if approx in ["SEOBNRv5HM", "SEOBNRv5EHM"]:
            chi1 = self.parameters["spin1z"]
            chi2 = self.parameters["spin2z"]

        elif approx == "SEOBNRv5PHM":
            chi1 = np.array(
                [
                    self.parameters["spin1x"],
                    self.parameters["spin1y"],
                    self.parameters["spin1z"],
                ]
            )
            chi2 = np.array(
                [
                    self.parameters["spin2x"],
                    self.parameters["spin2y"],
                    self.parameters["spin2z"],
                ]
            )
        else:
            assert False, "Incorrect approximant"

        omega_start = np.pi * fmin * Mtot * lal.MTSUN_SI
        omega_ref = np.pi * f_ref * Mtot * lal.MTSUN_SI

        q = m1 / m2  # Model convention q=m1/m2>=1
        if q < 1.0:
            q = 1 / q

        # Generate internal models, in geometrized units
        settings = {
            "dt": dt,
            "M": Mtot,
            "beta_approx": None,
            "polarizations_from_coprec": False,
        }
        if approx == "SEOBNRv5EHM":
            settings = settings | {
                "EccIC": self.parameters["EccIC"],
                "t_backwards": self.parameters["t_backwards"],
                "secular_bwd_int": self.parameters["secular_bwd_int"],
                "warning_bwd_int": self.parameters["warning_bwd_int"],
                "warning_secular_bwd_int": self.parameters["warning_secular_bwd_int"],
            }

        if "postadiabatic" in self.parameters:
            settings.update(postadiabatic=self.parameters["postadiabatic"])

            if "postadiabatic_type" in self.parameters:
                settings.update(
                    postadiabatic_type=self.parameters["postadiabatic_type"]
                )

        if "r_size_input" in self.parameters:
            settings.update(r_size_input=self.parameters["r_size_input"])

        if "initial_conditions" in self.parameters:
            settings.update(initial_conditions=self.parameters["initial_conditions"])

            if "initial_conditions_postadiabatic_type" in self.parameters:
                settings.update(
                    initial_conditions_postadiabatic_type=self.parameters[
                        "initial_conditions_postadiabatic_type"
                    ]
                )

        # Select mode array
        if self.parameters["mode_array"] is not None:
            settings["return_modes"] = self.parameters["mode_array"]
        if self.parameters["ModeArray"] is not None:
            settings["return_modes"] = self.parameters["ModeArray"]

        if "lmax_nyquist" in self.parameters:
            settings.update(lmax_nyquist=self.parameters["lmax_nyquist"])

        if "dA_dict" in self.parameters:
            settings.update(dA_dict=self.parameters["dA_dict"])
        if "dw_dict" in self.parameters:
            settings.update(dw_dict=self.parameters["dw_dict"])
        if "dTpeak" in self.parameters:
            settings.update(dTpeak=self.parameters["dTpeak"])
        if "da6" in self.parameters:
            settings.update(da6=self.parameters["da6"])
        if "ddSO" in self.parameters:
            settings.update(ddSO=self.parameters["ddSO"])
        if "domega_dict" in self.parameters:
            settings.update(domega_dict=self.parameters["domega_dict"])
        if "dtau_dict" in self.parameters:
            settings.update(dtau_dict=self.parameters["dtau_dict"])
        if "tol_PA" in self.parameters:
            settings.update(tol_PA=self.parameters["tol_PA"])
        if "rtol_ode" in self.parameters:
            settings.update(rtol_ode=self.parameters["rtol_ode"])
        if "atol_ode" in self.parameters:
            settings.update(atol_ode=self.parameters["atol_ode"])
        if "deltaT_sampling" in self.parameters:
            settings.update(deltaT_sampling=self.parameters["deltaT_sampling"])
        if "omega_prec_deviation" in self.parameters:
            settings.update(
                omega_prec_deviation=self.parameters["omega_prec_deviation"]
            )

        if "lmax" in self.parameters:
            settings.update(lmax=self.parameters["lmax"])

        enable_antisymmetric_modes = False
        if "enable_antisymmetric_modes" in self.parameters:
            settings.update(
                enable_antisymmetric_modes=self.parameters["enable_antisymmetric_modes"]
            )
            enable_antisymmetric_modes = self.parameters["enable_antisymmetric_modes"]

        if enable_antisymmetric_modes:
            if "antisymmetric_modes" in self.parameters:
                settings.update(
                    antisymmetric_modes=self.parameters["antisymmetric_modes"]
                )

            if "ivs_mrd" in self.parameters:
                settings.update(ivs_mrd=self.parameters["ivs_mrd"])

            if self.parameters.get("antisymmetric_modes_hm", False):
                settings.update(antisymmetric_modes=[(2, 2), (3, 3), (4, 4)])

            if "antisymmetric_fits_version" in self.parameters:
                settings.update(
                    antisymmetric_fits_version=self.parameters[
                        "antisymmetric_fits_version"
                    ]
                )

        settings.update(f_ref=self.parameters["f_ref"])
        times, h, self._model = generate_modes_opt(
            q,
            chi1,
            chi2,
            omega_start,
            approximant=approx,
            omega_ref=omega_ref,
            eccentricity=eccentricity,
            rel_anomaly=rel_anomaly,
            settings=settings,
            debug=True,
        )

        # Convert to physical units and LAL convention
        Mpc_to_meters = lal.PC_SI * 1e6
        times *= Mtot * lal.MTSUN_SI  # Physical times
        fac = (
            -1 * Mtot * lal.MRSUN_SI / (dist * Mpc_to_meters)
        )  # Minus sign to satisfy LAL convention

        hlm_dict = {}
        for ellm, mode in h.items():
            ell = int(ellm[0])
            if ellm[2] == "-":
                emm = -int(ellm[3])
            else:
                emm = int(ellm[2])
            hlm_dict[(ell, emm)] = fac * mode

        # If aligned-spin model, compute negative modes using equatorial symmetry
        if approx in ["SEOBNRv5HM", "SEOBNRv5EHM"]:
            for ellm, mode in h.items():
                ell = int(ellm[0])
                emm = int(ellm[2])
                hlm_dict[(ell, -emm)] = pow(-1, ell) * fac * np.conj(mode)

        # If masses are swapped to satisfy the m1/m2>=1 convention, this implies a
        # pi rotation on the orbital plane, which translates into a minus sign for the odd modes.
        if self.swap_masses is True:
            for ell, emm in hlm_dict.keys():
                if np.abs(emm) % 2 != 0:
                    hlm_dict[(ell, emm)] *= -1.0

        return times, hlm_dict

    def generate_td_polarizations(self):
        """
        Generate time-domain polarizations, returned as LAL REAL8TimeSeries.
        """

        incl = self.parameters["inclination"]
        phi = self.parameters["phi_ref"]

        if self.parameters["polarizations_from_coprec"] is False:
            hpc = 0.0
            times, hlm_dict = self.generate_td_modes()
            for ell, emm in hlm_dict:
                hlm = hlm_dict[(ell, emm)]
                ylm = lal.SpinWeightedSphericalHarmonic(
                    incl, np.pi / 2 - phi, -2, ell, emm
                )
                hpc += ylm * hlm

        else:
            fmin, dt = self.parameters["f22_start"], self.parameters["deltaT"]
            f_ref = self.parameters.get("f_ref")
            m1, m2 = self.parameters["mass1"], self.parameters["mass2"]
            Mtot = m1 + m2
            chi1 = np.array(
                [
                    self.parameters["spin1x"],
                    self.parameters["spin1y"],
                    self.parameters["spin1z"],
                ]
            )
            chi2 = np.array(
                [
                    self.parameters["spin2x"],
                    self.parameters["spin2y"],
                    self.parameters["spin2z"],
                ]
            )
            dist = self.parameters["distance"]
            omega_start = np.pi * fmin * Mtot * lal.MTSUN_SI
            omega_ref = np.pi * f_ref * Mtot * lal.MTSUN_SI
            q = m1 / m2
            # Generate internal polarizations, in geometrized units
            settings = {"dt": dt, "M": Mtot, "beta_approx": None}
            if "postadiabatic" in self.parameters:
                settings.update(postadiabatic=self.parameters["postadiabatic"])

                if "postadiabatic_type" in self.parameters:
                    settings.update(
                        postadiabatic_type=self.parameters["postadiabatic_type"]
                    )
            if "initial_conditions" in self.parameters:
                settings.update(
                    initial_conditions=self.parameters["initial_conditions"]
                )

                if "initial_conditions_postadiabatic_type" in self.parameters:
                    settings.update(
                        initial_conditions_postadiabatic_type=self.parameters[
                            "initial_conditions_postadiabatic_type"
                        ]
                    )

            if "r_size_input" in self.parameters:
                settings.update(r_size_input=self.parameters["r_size_input"])

            # Select mode array
            if self.parameters["mode_array"] is not None:
                settings["return_modes"] = self.parameters["mode_array"]
            if self.parameters["ModeArray"] is not None:
                settings["return_modes"] = self.parameters["ModeArray"]

            if "lmax_nyquist" in self.parameters:
                settings.update(lmax_nyquist=self.parameters["lmax_nyquist"])

            if "dA_dict" in self.parameters:
                settings.update(dA_dict=self.parameters["dA_dict"])
            if "dw_dict" in self.parameters:
                settings.update(dw_dict=self.parameters["dw_dict"])
            if "dTpeak" in self.parameters:
                settings.update(dTpeak=self.parameters["dTpeak"])
            if "da6" in self.parameters:
                settings.update(da6=self.parameters["da6"])
            if "ddSO" in self.parameters:
                settings.update(ddSO=self.parameters["ddSO"])
            if "domega_dict" in self.parameters:
                settings.update(domega_dict=self.parameters["domega_dict"])
            if "dtau_dict" in self.parameters:
                settings.update(dtau_dict=self.parameters["dtau_dict"])
            if "tol_PA" in self.parameters:
                settings.update(tol_PA=self.parameters["tol_PA"])
            if "rtol_ode" in self.parameters:
                settings.update(rtol_ode=self.parameters["rtol_ode"])
            if "atol_ode" in self.parameters:
                settings.update(atol_ode=self.parameters["atol_ode"])
            if "deltaT_sampling" in self.parameters:
                settings.update(deltaT_sampling=self.parameters["deltaT_sampling"])
            if "omega_prec_deviation" in self.parameters:
                settings.update(
                    omega_prec_deviation=self.parameters["omega_prec_deviation"]
                )

            if "lmax" in self.parameters:
                settings.update(lmax=self.parameters["lmax"])

            enable_antisymmetric_modes = False
            if "enable_antisymmetric_modes" in self.parameters:
                settings.update(
                    enable_antisymmetric_modes=self.parameters[
                        "enable_antisymmetric_modes"
                    ]
                )
                enable_antisymmetric_modes = self.parameters[
                    "enable_antisymmetric_modes"
                ]

            if enable_antisymmetric_modes:
                if "antisymmetric_modes" in self.parameters:
                    settings.update(
                        antisymmetric_modes=self.parameters["antisymmetric_modes"]
                    )

                if "ivs_mrd" in self.parameters:
                    settings.update(ivs_mrd=self.parameters["ivs_mrd"])

                if self.parameters.get("antisymmetric_modes_hm", False):
                    settings.update(antisymmetric_modes=[(2, 2), (3, 3), (4, 4)])

                if "antisymmetric_fits_version" in self.parameters:
                    settings.update(
                        antisymmetric_fits_version=self.parameters[
                            "antisymmetric_fits_version"
                        ]
                    )

            settings.update(f_ref=self.parameters["f_ref"])
            Mpc_to_meters = lal.PC_SI * 1e6
            fac = (
                -1 * Mtot * lal.MRSUN_SI / (dist * Mpc_to_meters)
            )  # Minus sign to satisfy LAL convention

            # inclination and phiref for polarizations
            settings.update(inclination=incl)
            # If masses are swapped to satisfy the m1/m2>=1 convention,
            # this implies a pi rotation on the orbital plane.
            if self.swap_masses:
                phi += np.pi
            settings.update(phiref=np.pi / 2 - phi)
            times, hpc, self._model = generate_prec_hpc_opt(
                q,
                chi1,
                chi2,
                omega_start,
                omega_ref=omega_ref,
                settings=settings,
                debug=True,
            )
            hpc *= fac
            times *= Mtot * lal.MTSUN_SI

        hp = np.real(hpc)
        hc = -np.imag(hpc)
        epoch = lal.LIGOTimeGPS(times[0])

        hp_lal = lal.CreateREAL8TimeSeries(
            "hplus", epoch, 0, self.parameters["deltaT"], lal.DimensionlessUnit, len(hp)
        )
        hc_lal = lal.CreateREAL8TimeSeries(
            "hcross",
            epoch,
            0,
            self.parameters["deltaT"],
            lal.DimensionlessUnit,
            len(hp),
        )

        hp_lal.data.data = hp
        hc_lal.data.data = hc

        return hp_lal, hc_lal

    # Procedure as in v4PHM in SimInspiralFD
    def generate_td_polarizations_conditioned_1(self):
        """
        Generate time-domain polarizations, with tappering at the beginning
        of the waveform, returned as LAL REAL8TimeSeries.
        """

        hp_lal, hc_lal = self.generate_td_polarizations()
        lalsim.SimInspiralREAL8WaveTaper(hp_lal.data, 1)
        lalsim.SimInspiralREAL8WaveTaper(hc_lal.data, 1)

        return hp_lal, hc_lal

    # General SimInspiralFD procedure, with extra time at the beginning
    def generate_td_polarizations_conditioned_2(self):
        """
        Generate conditioned time-domain polarizations as in
        SimInspiralTDfromTD routine.
        """

        if self.parameters["approximant"] == "SEOBNRv5EHM":
            raise ValueError(
                "The approximant 'SEOBNRv5EHM' cannot be used with the "
                "conditioning routine "
                "'generate_td_polarizations_conditioned_2'. Please use "
                "instead the routine "
                "'generate_td_polarizations_conditioned_1'."
            )

        extra_time_fraction = (
            0.1  # fraction of waveform duration to add as extra time for tapering
        )
        extra_cycles = (
            3.0  # more extra time measured in cycles at the starting frequency
        )

        f_min = self.parameters["f22_start"]
        m1 = self.parameters["mass1"]
        m2 = self.parameters["mass2"]
        S1z = self.parameters["spin1z"]
        S2z = self.parameters["spin2z"]
        original_f_min = f_min

        fisco = 1.0 / (pow(9.0, 1.5) * np.pi * (m1 + m2) * lal.MTSUN_SI)
        if f_min > fisco:
            f_min = fisco

        # upper bound on the chirp time starting at f_min
        tchirp = lalsim.SimInspiralChirpTimeBound(
            f_min, m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, S1z, S2z
        )
        # upper bound on the final black hole spin */
        spinkerr = lalsim.SimInspiralFinalBlackHoleSpinBound(S1z, S2z)
        # upper bound on the final plunge, merger, and ringdown time */
        tmerge = lalsim.SimInspiralMergeTimeBound(
            m1 * lal.MSUN_SI, m2 * lal.MSUN_SI
        ) + lalsim.SimInspiralRingdownTimeBound((m1 + m2) * lal.MSUN_SI, spinkerr)

        # extra time to include for all waveforms to take care of situations where the
        # frequency is close to merger (and is sweeping rapidly): this is a few cycles
        # at the low frequency
        textra = extra_cycles / f_min
        # compute a new lower frequency
        fstart = lalsim.SimInspiralChirpStartFrequencyBound(
            (1.0 + extra_time_fraction) * tchirp + tmerge + textra,
            m1 * lal.MSUN_SI,
            m2 * lal.MSUN_SI,
        )

        self.parameters["f22_start"] = fstart

        hp_lal, hc_lal = self.generate_td_polarizations()
        self.parameters["f22_start"] = original_f_min

        # condition the time domain waveform by tapering in the extra time at the
        # beginning and high-pass filtering above original f_min
        lalsim.SimInspiralTDConditionStage1(
            hp_lal, hc_lal, extra_time_fraction * tchirp + textra, original_f_min
        )

        # final tapering at the beginning and at the end to remove filter transients
        # waveform should terminate at a frequency >= Schwarzschild ISCO
        # so taper one cycle at this frequency at the end; should not make
        # any difference to IMR waveforms */
        fisco = 1.0 / (pow(6.0, 1.5) * np.pi * (m1 + m2) * lal.MTSUN_SI)
        lalsim.SimInspiralTDConditionStage2(hp_lal, hc_lal, f_min, fisco)

        return hp_lal, hc_lal

    def generate_fd_polarizations(self):
        """
        Generate Fourier-domain polarizations, returned as LAL COMPLEX16FrequencySeries

        Routine similar to LAL ``SimInspiralFD``.

        This routine assumes that ``f_max`` is the Nyquist frequency of the time-domain waveform,
        so that ``deltaT = 0.5 / f_max``. If ``deltaF`` is set to ``0`` then this routine computes
        a ``deltaF`` that is small enough to represent the Fourier transform of the time-domain
        waveform while ensuring the length of the time-domain signal if a power of 2.

        If ``deltaF`` is specified but ``f_max / deltaF`` is not a power of 2, then ``f_max`` is
        increased so that ``f_max / deltaF`` is the next power of 2. (If the user wishes to discard
        the extra high frequency content, this must be done separately.)

        .. warning::

            The user should take care to ensure that ``deltaF`` is sufficiently small to contain the full
            signal (time series duration ``= 1 / deltaF``). If the provided ``deltaF`` is too large the
            signal will be abruptly truncated.

            Similarly, if the provided ``f_max`` is less than the ringdown frequency the underlying
            waveform generator may raise an error. If not, the frequency domain signal will be aliased
            in the frequency domain.

        """

        # Adjust deltaT depending on sampling rate
        fmax = self.parameters["f_max"]
        f_nyquist = fmax
        deltaF = 0
        if "deltaF" in self.parameters.keys():
            deltaF = self.parameters["deltaF"]

        if deltaF != 0:
            n = int(np.round(fmax / deltaF))
            if n & (n - 1):
                chirplen_exp = np.frexp(n)
                f_nyquist = np.ldexp(1, int(chirplen_exp[1])) * deltaF

        deltaT = 0.5 / f_nyquist
        self.parameters["deltaT"] = deltaT

        # Generate conditioned TD polarizations
        if (
            self.parameters["conditioning"] == 2
            and self.parameters["approximant"] != "SEOBNRv5EHM"
        ):
            hp, hc = self.generate_td_polarizations_conditioned_2()
        else:
            hp, hc = self.generate_td_polarizations_conditioned_1()

        # Adjust signal duration
        if deltaF == 0:
            chirplen = hp.data.length
            chirplen_exp = np.frexp(chirplen)
            chirplen = int(np.ldexp(1, chirplen_exp[1]))
            deltaF = 1.0 / (chirplen * deltaT)
            self.parameters["deltaF"] = deltaF

        else:
            chirplen = int(1.0 / (deltaF * deltaT))

        # resize waveforms to the required length
        lal.ResizeREAL8TimeSeries(hp, hp.data.length - chirplen, chirplen)
        lal.ResizeREAL8TimeSeries(hc, hc.data.length - chirplen, chirplen)

        # FFT - Using LAL routines
        hptilde = lal.CreateCOMPLEX16FrequencySeries(
            "FD H_PLUS",
            hp.epoch,
            0.0,
            deltaF,
            lal.DimensionlessUnit,
            int(chirplen / 2.0 + 1),
        )
        hctilde = lal.CreateCOMPLEX16FrequencySeries(
            "FD H_CROSS",
            hc.epoch,
            0.0,
            deltaF,
            lal.DimensionlessUnit,
            int(chirplen / 2.0 + 1),
        )

        plan = lal.CreateForwardREAL8FFTPlan(chirplen, 0)
        lal.REAL8TimeFreqFFT(hctilde, hc, plan)
        lal.REAL8TimeFreqFFT(hptilde, hp, plan)

        return hptilde, hctilde
