"""Near field to far field transformation plugin
"""
from __future__ import annotations
from typing import Dict, Tuple, Union, List
import numpy as np
import xarray as xr
import pydantic.v1 as pydantic

from rich.progress import track

from .data.data_array import FieldProjectionAngleDataArray, FieldProjectionCartesianDataArray
from .data.data_array import FieldProjectionKSpaceDataArray
from .data.monitor_data import FieldData
from .data.monitor_data import AbstractFieldProjectionData, FieldProjectionAngleData
from .data.monitor_data import FieldProjectionCartesianData, FieldProjectionKSpaceData
from .data.sim_data import SimulationData
from .monitor import FieldProjectionSurface
from .monitor import FieldMonitor, AbstractFieldProjectionMonitor, FieldProjectionAngleMonitor
from .monitor import FieldProjectionCartesianMonitor, FieldProjectionKSpaceMonitor
from .types import Direction, Coordinate, ArrayComplex4D
from .medium import MediumType
from .base import Tidy3dBaseModel, cached_property, skip_if_fields_missing
from ..exceptions import SetupError
from ..constants import C_0, MICROMETER, ETA_0, EPSILON_0, MU_0
from ..log import get_logging_console

# Default number of points per wavelength in the background medium to use for resampling fields.
PTS_PER_WVL = 10

# Numpy float array and related array types

ArrayLikeN2F = Union[float, Tuple[float, ...], ArrayComplex4D]



[docs]
class FieldProjector(Tidy3dBaseModel):
    """
    Projection of near fields to points on a given observation grid.

    .. TODO make images to illustrate this

    See Also
    --------

    :class:`FieldProjectionAngleMonitor
        :class:`Monitor` that samples electromagnetic near fields in the frequency domain
        and projects them at given observation angles.`

    **Notebooks**:
        * `Performing near field to far field projections <../../notebooks/FieldProjections.html>`_
    """

    sim_data: SimulationData = pydantic.Field(
        ...,
        title="Simulation data",
        description="Container for simulation data containing the near field monitors.",
    )

    surfaces: Tuple[FieldProjectionSurface, ...] = pydantic.Field(
        ...,
        title="Surface monitor with direction",
        description="Tuple of each :class:`.FieldProjectionSurface` to use as source of "
        "near field.",
    )

    pts_per_wavelength: Union[int, type(None)] = pydantic.Field(
        PTS_PER_WVL,
        title="Points per wavelength",
        description="Number of points per wavelength in the background medium with which "
        "to discretize the surface monitors for the projection. If ``None``, fields will "
        "will not resampled, but will still be colocated.",
    )

    origin: Coordinate = pydantic.Field(
        None,
        title="Local origin",
        description="Local origin used for defining observation points. If ``None``, uses the "
        "average of the centers of all surface monitors.",
        units=MICROMETER,
    )

    currents: Dict[str, xr.Dataset] = pydantic.Field(
        None,
        title="Surface current densities",
        description="Dictionary mapping monitor name to an ``xarray.Dataset`` storing the "
        "surface current densities.",
    )


[docs]
    @pydantic.validator("origin", always=True)
    @skip_if_fields_missing(["surfaces"])
    def set_origin(cls, val, values):
        """Sets .origin as the average of centers of all surface monitors if not provided."""
        if val is None:
            surfaces = values.get("surfaces")
            val = np.array([surface.monitor.center for surface in surfaces])
            return tuple(np.mean(val, axis=0))
        return val


    @cached_property
    def medium(self) -> MediumType:
        """Medium into which fields are to be projected."""
        sim = self.sim_data.simulation
        monitor = self.surfaces[0].monitor
        return sim.monitor_medium(monitor)

    @cached_property
    def frequencies(self) -> List[float]:
        """Return the list of frequencies associated with the field monitors."""
        return self.surfaces[0].monitor.freqs


[docs]
    @classmethod
    def from_near_field_monitors(
        cls,
        sim_data: SimulationData,
        near_monitors: List[FieldMonitor],
        normal_dirs: List[Direction],
        pts_per_wavelength: int = PTS_PER_WVL,
        origin: Coordinate = None,
    ):
        """Constructs :class:`FieldProjection` from a list of surface monitors and their directions.

        Parameters
        ----------
        sim_data : :class:`.SimulationData`
            Container for simulation data containing the near field monitors.
        near_monitors : List[:class:`.FieldMonitor`]
            Tuple of :class:`.FieldMonitor` objects on which near fields will be sampled.
        normal_dirs : List[:class:`.Direction`]
            Tuple containing the :class:`.Direction` of the normal to each surface monitor
            w.r.t. to the positive x, y or z unit vectors. Must have the same length as monitors.
        pts_per_wavelength : int = 10
            Number of points per wavelength with which to discretize the
            surface monitors for the projection. If ``None``, fields will not be resampled.
        origin : :class:`.Coordinate`
            Local origin used for defining observation points. If ``None``, uses the
            average of the centers of all surface monitors.
        """

        if len(near_monitors) != len(normal_dirs):
            raise SetupError(
                f"Number of monitors ({len(near_monitors)}) does not equal "
                f"the number of directions ({len(normal_dirs)})."
            )

        surfaces = [
            FieldProjectionSurface(monitor=monitor, normal_dir=normal_dir)
            for monitor, normal_dir in zip(near_monitors, normal_dirs)
        ]

        return cls(
            sim_data=sim_data,
            surfaces=surfaces,
            pts_per_wavelength=pts_per_wavelength,
            origin=origin,
        )


    @cached_property
    def currents(self):
        """Sets the surface currents."""
        sim_data = self.sim_data
        surfaces = self.surfaces
        pts_per_wavelength = self.pts_per_wavelength
        medium = self.medium

        surface_currents = {}
        for surface in surfaces:
            current_data = self.compute_surface_currents(
                sim_data, surface, medium, pts_per_wavelength
            )

            # shift source coordinates relative to the local origin
            for name, origin in zip(["x", "y", "z"], self.origin):
                current_data[name] = current_data[name] - origin

            surface_currents[surface.monitor.name] = current_data

        return surface_currents


[docs]
    @staticmethod
    def compute_surface_currents(
        sim_data: SimulationData,
        surface: FieldProjectionSurface,
        medium: MediumType,
        pts_per_wavelength: int = PTS_PER_WVL,
    ) -> xr.Dataset:
        """Returns resampled surface current densities associated with the surface monitor.

        Parameters
        ----------
        sim_data : :class:`.SimulationData`
            Container for simulation data containing the near field monitors.
        surface: :class:`.FieldProjectionSurface`
            :class:`.FieldProjectionSurface` to use as source of near field.
        medium : :class:`.MediumType`
            Background medium through which to project fields.
        pts_per_wavelength : int = 10
            Number of points per wavelength with which to discretize the
            surface monitors for the projection. If ``None``, fields will not be
            resampled, but will still be colocated.

        Returns
        -------
        xarray.Dataset
            Colocated surface current densities for the given surface.
        """

        monitor_name = surface.monitor.name
        if monitor_name not in sim_data.monitor_data.keys():
            raise SetupError(f"No data for monitor named '{monitor_name}' found in sim_data.")

        field_data = sim_data[monitor_name]

        currents = FieldProjector._fields_to_currents(field_data, surface)
        currents = FieldProjector._resample_surface_currents(
            currents, sim_data, surface, medium, pts_per_wavelength
        )

        return currents


    @staticmethod
    def _fields_to_currents(field_data: FieldData, surface: FieldProjectionSurface) -> FieldData:
        """Returns surface current densities associated with a given :class:`.FieldData` object.

        Parameters
        ----------
        field_data : :class:`.FieldData`
            Container for field data associated with the given near field surface.
        surface: :class:`.FieldProjectionSurface`
            :class:`.FieldProjectionSurface` to use as source of near field.

        Returns
        -------
        :class:`.FieldData`
            Surface current densities for the given surface.
        """

        # figure out which field components are tangential or normal to the monitor
        _, (cmp_1, cmp_2) = surface.monitor.pop_axis(("x", "y", "z"), axis=surface.axis)

        signs = np.array([-1, 1])
        if surface.axis % 2 != 0:
            signs *= -1
        if surface.normal_dir == "-":
            signs *= -1

        E1 = "E" + cmp_1
        E2 = "E" + cmp_2
        H1 = "H" + cmp_1
        H2 = "H" + cmp_2

        surface_currents = {}

        surface_currents[E2] = field_data.field_components[H1] * signs[1]
        surface_currents[E1] = field_data.field_components[H2] * signs[0]

        surface_currents[H2] = field_data.field_components[E1] * signs[0]
        surface_currents[H1] = field_data.field_components[E2] * signs[1]

        new_monitor = surface.monitor.copy(update=dict(fields=[E1, E2, H1, H2]))

        return FieldData(
            monitor=new_monitor,
            symmetry=field_data.symmetry,
            symmetry_center=field_data.symmetry_center,
            grid_expanded=field_data.grid_expanded,
            **surface_currents,
        )

    @staticmethod
    def _resample_surface_currents(
        currents: xr.Dataset,
        sim_data: SimulationData,
        surface: FieldProjectionSurface,
        medium: MediumType,
        pts_per_wavelength: int = PTS_PER_WVL,
    ) -> xr.Dataset:
        """Returns the surface current densities associated with the surface monitor.

        Parameters
        ----------
        currents : xarray.Dataset
            Surface currents defined on the original Yee grid.
        sim_data : :class:`.SimulationData`
            Container for simulation data containing the near field monitors.
        surface: :class:`.FieldProjectionSurface`
            :class:`.FieldProjectionSurface` to use as source of near field.
        medium : :class:`.MediumType`
            Background medium through which to project fields.
        pts_per_wavelength : int = 10
            Number of points per wavelength with which to discretize the
            surface monitors for the projection. If ``None``, fields will not be
            resampled, but will still be colocated.

        Returns
        -------
        xarray.Dataset
            Colocated surface current densities for the given surface.
        """

        # colocate surface currents on a regular grid of points on the monitor based on wavelength
        colocation_points = [None] * 3
        colocation_points[surface.axis] = surface.monitor.center[surface.axis]

        # use the highest frequency associated with the monitor to resample the surface currents
        frequency = max(surface.monitor.freqs)
        eps_complex = medium.eps_model(frequency)
        index_n, _ = medium.eps_complex_to_nk(eps_complex)
        wavelength = C_0 / frequency / index_n

        _, idx_uv = surface.monitor.pop_axis((0, 1, 2), axis=surface.axis)

        for idx in idx_uv:
            # pick sample points on the monitor and handle the possibility of an "infinite" monitor
            start = np.maximum(
                surface.monitor.center[idx] - surface.monitor.size[idx] / 2.0,
                sim_data.simulation.center[idx] - sim_data.simulation.size[idx] / 2.0,
            )
            stop = np.minimum(
                surface.monitor.center[idx] + surface.monitor.size[idx] / 2.0,
                sim_data.simulation.center[idx] + sim_data.simulation.size[idx] / 2.0,
            )

            if pts_per_wavelength is None:
                points = sim_data.simulation.grid.boundaries.to_list[idx]
                points[np.argwhere(points < start)] = start
                points[np.argwhere(points > stop)] = stop
                colocation_points[idx] = np.unique(points)
            else:
                size = stop - start
                num_pts = int(np.ceil(pts_per_wavelength * size / wavelength))
                points = np.linspace(start, stop, num_pts)
                colocation_points[idx] = points

        for idx, points in enumerate(colocation_points):
            if np.array(points).size == 1:
                colocation_points[idx] = None

        currents = currents.colocate(*colocation_points)
        return currents


[docs]
    def integrate_2d(
        self,
        function: np.ndarray,
        phase: np.ndarray,
        pts_u: np.ndarray,
        pts_v: np.ndarray,
    ):
        """Trapezoidal integration in two dimensions."""
        return np.trapz(np.trapz(np.squeeze(function) * phase, pts_u, axis=0), pts_v, axis=0)


    def _far_fields_for_surface(
        self,
        frequency: float,
        theta: ArrayLikeN2F,
        phi: ArrayLikeN2F,
        surface: FieldProjectionSurface,
        currents: xr.Dataset,
        medium: MediumType,
    ):
        """Compute far fields at an angle in spherical coordinates
        for a given set of surface currents and observation angles.

        Parameters
        ----------
        frequency : float
            Frequency to select from each :class:`.FieldMonitor` to use for projection.
            Must be a frequency stored in each :class:`FieldMonitor`.
        theta : Union[float, Tuple[float, ...], np.ndarray]
            Polar angles (rad) downward from x=y=0 line relative to the local origin.
        phi : Union[float, Tuple[float, ...], np.ndarray]
            Azimuthal (rad) angles from y=z=0 line relative to the local origin.
        surface: :class:`FieldProjectionSurface`
            :class:`FieldProjectionSurface` object to use as source of near field.
        currents : xarray.Dataset
            xarray Dataset containing surface currents associated with the surface monitor.
        medium : :class:`.MediumType`
            Background medium through which to project fields.

        Returns
        -------
        tuple(numpy.ndarray[float], ...)
            ``Er``, ``Etheta``, ``Ephi``, ``Hr``, ``Htheta``, ``Hphi`` for the given surface.
        """
        pts = [currents[name].values for name in ["x", "y", "z"]]

        try:
            currents_f = currents.sel(f=frequency)
        except Exception as e:
            raise SetupError(
                f"Frequency {frequency} not found in fields for monitor '{surface.monitor.name}'."
            ) from e

        idx_w, idx_uv = surface.monitor.pop_axis((0, 1, 2), axis=surface.axis)
        _, source_names = surface.monitor.pop_axis(("x", "y", "z"), axis=surface.axis)

        idx_u, idx_v = idx_uv
        cmp_1, cmp_2 = source_names

        theta = np.atleast_1d(theta)
        phi = np.atleast_1d(phi)

        sin_theta = np.sin(theta)
        cos_theta = np.cos(theta)
        sin_phi = np.sin(phi)
        cos_phi = np.cos(phi)

        J = np.zeros((3, len(theta), len(phi)), dtype=complex)
        M = np.zeros_like(J)

        phase = [None] * 3
        propagation_factor = -1j * AbstractFieldProjectionData.wavenumber(
            medium=medium, frequency=frequency
        )

        def integrate_for_one_theta(i_th: int):
            """Perform integration for a given theta angle index"""

            for j_ph in np.arange(len(phi)):
                phase[0] = np.exp(propagation_factor * pts[0] * sin_theta[i_th] * cos_phi[j_ph])
                phase[1] = np.exp(propagation_factor * pts[1] * sin_theta[i_th] * sin_phi[j_ph])
                phase[2] = np.exp(propagation_factor * pts[2] * cos_theta[i_th])

                phase_ij = phase[idx_u][:, None] * phase[idx_v][None, :] * phase[idx_w]

                J[idx_u, i_th, j_ph] = self.integrate_2d(
                    currents_f[f"E{cmp_1}"].values, phase_ij, pts[idx_u], pts[idx_v]
                )

                J[idx_v, i_th, j_ph] = self.integrate_2d(
                    currents_f[f"E{cmp_2}"].values, phase_ij, pts[idx_u], pts[idx_v]
                )

                M[idx_u, i_th, j_ph] = self.integrate_2d(
                    currents_f[f"H{cmp_1}"].values, phase_ij, pts[idx_u], pts[idx_v]
                )

                M[idx_v, i_th, j_ph] = self.integrate_2d(
                    currents_f[f"H{cmp_2}"].values, phase_ij, pts[idx_u], pts[idx_v]
                )

        if len(theta) < 2:
            integrate_for_one_theta(0)
        else:
            for i_th in track(
                np.arange(len(theta)),
                description=f"Processing surface monitor '{surface.monitor.name}'...",
                console=get_logging_console(),
            ):
                integrate_for_one_theta(i_th)

        cos_th_cos_phi = cos_theta[:, None] * cos_phi[None, :]
        cos_th_sin_phi = cos_theta[:, None] * sin_phi[None, :]

        # Ntheta (8.33a)
        Ntheta = J[0] * cos_th_cos_phi + J[1] * cos_th_sin_phi - J[2] * sin_theta[:, None]

        # Nphi (8.33b)
        Nphi = -J[0] * sin_phi[None, :] + J[1] * cos_phi[None, :]

        # Ltheta  (8.34a)
        Ltheta = M[0] * cos_th_cos_phi + M[1] * cos_th_sin_phi - M[2] * sin_theta[:, None]

        # Lphi  (8.34b)
        Lphi = -M[0] * sin_phi[None, :] + M[1] * cos_phi[None, :]

        eta = ETA_0 / np.sqrt(medium.eps_model(frequency))

        Etheta = -(Lphi + eta * Ntheta)
        Ephi = Ltheta - eta * Nphi
        Er = np.zeros_like(Ephi)
        Htheta = -Ephi / eta
        Hphi = Etheta / eta
        Hr = np.zeros_like(Hphi)

        return Er, Etheta, Ephi, Hr, Htheta, Hphi


[docs]
    @staticmethod
    def apply_window_to_currents(
        proj_monitor: AbstractFieldProjectionMonitor, currents: xr.Dataset
    ) -> xr.Dataset:
        """Apply windowing function to the surface currents."""
        if proj_monitor.size.count(0.0) == 0:
            return currents
        if proj_monitor.window_size == (0, 0):
            return currents

        pts = [currents[name].values for name in ["x", "y", "z"]]

        custom_bounds = [
            [pts[i][0] for i in range(3)],
            [pts[i][-1] for i in range(3)],
        ]

        window_size, window_minus, window_plus = proj_monitor.window_parameters(
            custom_bounds=custom_bounds
        )

        new_currents = currents.copy(deep=True)
        for dim, (dim_name, points) in enumerate(zip("xyz", pts)):
            window_fn = proj_monitor.window_function(
                points=points,
                window_size=window_size,
                window_minus=window_minus,
                window_plus=window_plus,
                dim=dim,
            )
            window_data = xr.DataArray(
                window_fn,
                dims=[dim_name],
                coords=[points],
            )
            new_currents *= window_data

        return new_currents



[docs]
    def project_fields(
        self, proj_monitor: AbstractFieldProjectionMonitor
    ) -> AbstractFieldProjectionData:
        """Compute projected fields.

        Parameters
        ----------
        proj_monitor : :class:`.AbstractFieldProjectionMonitor`
            Instance of :class:`.AbstractFieldProjectionMonitor` defining the projection
            observation grid.

        Returns
        -------
        :class:`.AbstractFieldProjectionData`
            Data structure with ``Er``, ``Etheta``, ``Ephi``, ``Hr``, ``Htheta``, ``Hphi``.
        """
        if isinstance(proj_monitor, FieldProjectionAngleMonitor):
            return self._project_fields_angular(proj_monitor)
        if isinstance(proj_monitor, FieldProjectionCartesianMonitor):
            return self._project_fields_cartesian(proj_monitor)
        return self._project_fields_kspace(proj_monitor)


    def _project_fields_angular(
        self, monitor: FieldProjectionAngleMonitor
    ) -> FieldProjectionAngleData:
        """Compute projected fields on an angle-based grid in spherical coordinates.

        Parameters
        ----------
        monitor : :class:`.FieldProjectionAngleMonitor`
            Instance of :class:`.FieldProjectionAngleMonitor` defining the projection
            observation grid.

        Returns
        -------
        :class:.`FieldProjectionAngleData`
            Data structure with ``Er``, ``Etheta``, ``Ephi``, ``Hr``, ``Htheta``, ``Hphi``.
        """
        freqs = np.atleast_1d(self.frequencies)
        theta = np.atleast_1d(monitor.theta)
        phi = np.atleast_1d(monitor.phi)

        # compute projected fields for the dataset associated with each monitor
        field_names = ("Er", "Etheta", "Ephi", "Hr", "Htheta", "Hphi")
        fields = [
            np.zeros((1, len(theta), len(phi), len(freqs)), dtype=complex) for _ in field_names
        ]

        medium = monitor.medium if monitor.medium else self.medium
        k = AbstractFieldProjectionData.wavenumber(medium=medium, frequency=freqs)
        phase = np.atleast_1d(
            AbstractFieldProjectionData.propagation_phase(dist=monitor.proj_distance, k=k)
        )

        for surface in self.surfaces:
            # apply windowing to currents
            currents = self.apply_window_to_currents(monitor, self.currents[surface.monitor.name])

            if monitor.far_field_approx:
                for idx_f, frequency in enumerate(freqs):
                    _fields = self._far_fields_for_surface(
                        frequency=frequency,
                        theta=theta,
                        phi=phi,
                        surface=surface,
                        currents=currents,
                        medium=medium,
                    )
                    for field, _field in zip(fields, _fields):
                        field[..., idx_f] += _field * phase[idx_f]
            else:
                iter_coords = [
                    ([_theta, _phi], [i, j])
                    for i, _theta in enumerate(theta)
                    for j, _phi in enumerate(phi)
                ]
                for (_theta, _phi), (i, j) in track(
                    iter_coords,
                    description=f"Processing surface monitor '{surface.monitor.name}'...",
                    console=get_logging_console(),
                ):
                    _x, _y, _z = monitor.sph_2_car(monitor.proj_distance, _theta, _phi)
                    _fields = self._fields_for_surface_exact(
                        x=_x, y=_y, z=_z, surface=surface, currents=currents, medium=medium
                    )
                    for field, _field in zip(fields, _fields):
                        field[0, i, j, :] += _field

        coords = {"r": np.atleast_1d(monitor.proj_distance), "theta": theta, "phi": phi, "f": freqs}
        fields = {
            name: FieldProjectionAngleDataArray(field, coords=coords)
            for name, field in zip(field_names, fields)
        }
        return FieldProjectionAngleData(
            monitor=monitor, projection_surfaces=self.surfaces, medium=medium, **fields
        )

    def _project_fields_cartesian(
        self, monitor: FieldProjectionCartesianMonitor
    ) -> FieldProjectionCartesianData:
        """Compute projected fields on a Cartesian grid in spherical coordinates.

        Parameters
        ----------
        monitor : :class:`.FieldProjectionCartesianMonitor`
            Instance of :class:`.FieldProjectionCartesianMonitor` defining the projection
            observation grid.

        Returns
        -------
        :class:.`FieldProjectionCartesianData`
            Data structure with ``Er``, ``Etheta``, ``Ephi``, ``Hr``, ``Htheta``, ``Hphi``.
        """
        freqs = np.atleast_1d(self.frequencies)
        x, y, z = monitor.unpop_axis(
            monitor.proj_distance, (monitor.x, monitor.y), axis=monitor.proj_axis
        )
        x, y, z = list(map(np.atleast_1d, [x, y, z]))

        # compute projected fields for the dataset associated with each monitor
        field_names = ("Er", "Etheta", "Ephi", "Hr", "Htheta", "Hphi")
        fields = [
            np.zeros((len(x), len(y), len(z), len(freqs)), dtype=complex) for _ in field_names
        ]

        medium = monitor.medium if monitor.medium else self.medium
        wavenumber = AbstractFieldProjectionData.wavenumber(medium=medium, frequency=freqs)

        # Zip together all combinations of observation points for better progress tracking
        iter_coords = [
            ([_x, _y, _z], [i, j, k])
            for i, _x in enumerate(x)
            for j, _y in enumerate(y)
            for k, _z in enumerate(z)
        ]

        for (_x, _y, _z), (i, j, k) in track(
            iter_coords, description="Computing projected fields", console=get_logging_console()
        ):
            r, theta, phi = monitor.car_2_sph(_x, _y, _z)
            phase = np.atleast_1d(
                AbstractFieldProjectionData.propagation_phase(dist=r, k=wavenumber)
            )

            for surface in self.surfaces:
                # apply windowing to currents
                currents = self.apply_window_to_currents(
                    monitor, self.currents[surface.monitor.name]
                )

                if monitor.far_field_approx:
                    for idx_f, frequency in enumerate(freqs):
                        _fields = self._far_fields_for_surface(
                            frequency=frequency,
                            theta=theta,
                            phi=phi,
                            surface=surface,
                            currents=currents,
                            medium=medium,
                        )
                        for field, _field in zip(fields, _fields):
                            field[i, j, k, idx_f] += _field * phase[idx_f]
                else:
                    _fields = self._fields_for_surface_exact(
                        x=_x, y=_y, z=_z, surface=surface, currents=currents, medium=medium
                    )
                    for field, _field in zip(fields, _fields):
                        field[i, j, k, :] += _field

        coords = {"x": x, "y": y, "z": z, "f": freqs}
        fields = {
            name: FieldProjectionCartesianDataArray(field, coords=coords)
            for name, field in zip(field_names, fields)
        }
        return FieldProjectionCartesianData(
            monitor=monitor, projection_surfaces=self.surfaces, medium=medium, **fields
        )

    def _project_fields_kspace(
        self, monitor: FieldProjectionKSpaceMonitor
    ) -> FieldProjectionKSpaceData:
        """Compute projected fields on a k-space grid in spherical coordinates.

        Parameters
        ----------
        monitor : :class:`.FieldProjectionKSpaceMonitor`
            Instance of :class:`.FieldProjectionKSpaceMonitor` defining the projection
            observation grid.

        Returns
        -------
        :class:.`FieldProjectionKSpaceData`
            Data structure with ``Er``, ``Etheta``, ``Ephi``, ``Hr``, ``Htheta``, ``Hphi``.
        """
        freqs = np.atleast_1d(self.frequencies)
        ux = np.atleast_1d(monitor.ux)
        uy = np.atleast_1d(monitor.uy)

        # compute projected fields for the dataset associated with each monitor
        field_names = ("Er", "Etheta", "Ephi", "Hr", "Htheta", "Hphi")
        fields = [np.zeros((len(ux), len(uy), 1, len(freqs)), dtype=complex) for _ in field_names]

        medium = monitor.medium if monitor.medium else self.medium
        k = AbstractFieldProjectionData.wavenumber(medium=medium, frequency=freqs)
        phase = np.atleast_1d(
            AbstractFieldProjectionData.propagation_phase(dist=monitor.proj_distance, k=k)
        )

        # Zip together all combinations of observation points for better progress tracking
        iter_coords = [([_ux, _uy], [i, j]) for i, _ux in enumerate(ux) for j, _uy in enumerate(uy)]

        for (_ux, _uy), (i, j) in track(
            iter_coords, description="Computing projected fields", console=get_logging_console()
        ):
            theta, phi = monitor.kspace_2_sph(_ux, _uy, monitor.proj_axis)

            for surface in self.surfaces:
                # apply windowing to currents
                currents = self.apply_window_to_currents(
                    monitor, self.currents[surface.monitor.name]
                )

                if monitor.far_field_approx:
                    for idx_f, frequency in enumerate(freqs):
                        _fields = self._far_fields_for_surface(
                            frequency=frequency,
                            theta=theta,
                            phi=phi,
                            surface=surface,
                            currents=currents,
                            medium=medium,
                        )
                        for field, _field in zip(fields, _fields):
                            field[i, j, 0, idx_f] += _field * phase[idx_f]

                else:
                    _x, _y, _z = monitor.sph_2_car(monitor.proj_distance, theta, phi)
                    _fields = self._fields_for_surface_exact(
                        x=_x, y=_y, z=_z, surface=surface, currents=currents, medium=medium
                    )
                    for field, _field in zip(fields, _fields):
                        field[i, j, 0, :] += _field

        coords = {
            "ux": np.array(monitor.ux),
            "uy": np.array(monitor.uy),
            "r": np.atleast_1d(monitor.proj_distance),
            "f": freqs,
        }
        fields = {
            name: FieldProjectionKSpaceDataArray(field, coords=coords)
            for name, field in zip(field_names, fields)
        }
        return FieldProjectionKSpaceData(
            monitor=monitor, projection_surfaces=self.surfaces, medium=medium, **fields
        )

    """Exact projections"""

    def _fields_for_surface_exact(
        self,
        x: float,
        y: float,
        z: float,
        surface: FieldProjectionSurface,
        currents: xr.Dataset,
        medium: MediumType,
    ):
        """Compute projected fields in spherical coordinates at a given projection point on a
        Cartesian grid for a given set of surface currents using the exact homogeneous medium
        Green's function without geometric approximations.

        Parameters
        ----------
        x : float
            Observation point x-coordinate (microns) relative to the local origin.
        y : float
            Observation point y-coordinate (microns) relative to the local origin.
        z : float
            Observation point z-coordinate (microns) relative to the local origin.
        surface: :class:`FieldProjectionSurface`
            :class:`FieldProjectionSurface` object to use as source of near field.
        currents : xarray.Dataset
            xarray Dataset containing surface currents associated with the surface monitor.
        medium : :class:`.MediumType`
            Background medium through which to project fields.

        Returns
        -------
        tuple(np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray)
            ``Er``, ``Etheta``, ``Ephi``, ``Hr``, ``Htheta``, ``Hphi`` projected fields for
            each frequency.
        """
        freqs = np.array(self.frequencies)
        i_omega = 1j * 2.0 * np.pi * freqs[None, None, None, :]
        wavenumber = AbstractFieldProjectionData.wavenumber(frequency=freqs, medium=medium)
        wavenumber = wavenumber[None, None, None, :]  # add space dimensions

        eps_complex = medium.eps_model(frequency=freqs)
        epsilon = EPSILON_0 * eps_complex[None, None, None, :]

        # source points
        pts = [currents[name].values for name in ["x", "y", "z"]]

        # transform the coordinate system so that the origin is at the source point
        # then the observation points in the new system are:
        x_new, y_new, z_new = (pt_obs - pt_src for pt_src, pt_obs in zip(pts, [x, y, z]))

        # tangential source components to use
        idx_w, idx_uv = surface.monitor.pop_axis((0, 1, 2), axis=surface.axis)
        _, source_names = surface.monitor.pop_axis(("x", "y", "z"), axis=surface.axis)

        idx_u, idx_v = idx_uv
        cmp_1, cmp_2 = source_names

        # set the surface current density Cartesian components
        J = [np.atleast_1d(0)] * 3
        M = [np.atleast_1d(0)] * 3

        J[idx_u] = currents[f"E{cmp_1}"].values
        J[idx_v] = currents[f"E{cmp_2}"].values
        J[idx_w] = np.zeros_like(J[idx_u])
        M[idx_u] = currents[f"H{cmp_1}"].values
        M[idx_v] = currents[f"H{cmp_2}"].values
        M[idx_w] = np.zeros_like(M[idx_u])

        # observation point in the new spherical system
        r, theta_obs, phi_obs = surface.monitor.car_2_sph(
            x_new[:, None, None, None], y_new[None, :, None, None], z_new[None, None, :, None]
        )

        # angle terms
        sin_theta = np.sin(theta_obs)
        cos_theta = np.cos(theta_obs)
        sin_phi = np.sin(phi_obs)
        cos_phi = np.cos(phi_obs)

        # Green's function and terms related to its derivatives
        ikr = 1j * wavenumber * r
        G = np.exp(ikr) / (4.0 * np.pi * r)
        dG_dr = G * (ikr - 1.0) / r
        d2G_dr2 = dG_dr * (ikr - 1.0) / r + G / (r**2)

        # operations between unit vectors and currents
        def r_x_current(current: Tuple[np.ndarray, ...]) -> Tuple[np.ndarray, ...]:
            """Cross product between the r unit vector and the current."""
            return [
                sin_theta * sin_phi * current[2] - cos_theta * current[1],
                cos_theta * current[0] - sin_theta * cos_phi * current[2],
                sin_theta * cos_phi * current[1] - sin_theta * sin_phi * current[0],
            ]

        def r_dot_current(current: Tuple[np.ndarray, ...]) -> np.ndarray:
            """Dot product between the r unit vector and the current."""
            return (
                sin_theta * cos_phi * current[0]
                + sin_theta * sin_phi * current[1]
                + cos_theta * current[2]
            )

        def r_dot_current_dtheta(current: Tuple[np.ndarray, ...]) -> np.ndarray:
            """Theta derivative of the dot product between the r unit vector and the current."""
            return (
                cos_theta * cos_phi * current[0]
                + cos_theta * sin_phi * current[1]
                - sin_theta * current[2]
            )

        def r_dot_current_dphi_div_sin_theta(current: Tuple[np.ndarray, ...]) -> np.ndarray:
            """Phi derivative of the dot product between the r unit vector and the current,
            analytically divided by sin theta."""
            return -sin_phi * current[0] + cos_phi * current[1]

        def grad_Gr_r_dot_current(current: Tuple[np.ndarray, ...]) -> Tuple[np.ndarray, ...]:
            """Gradient of the product of the gradient of the Green's function and the dot product
            between the r unit vector and the current."""
            temp = [
                d2G_dr2 * r_dot_current(current),
                dG_dr * r_dot_current_dtheta(current) / r,
                dG_dr * r_dot_current_dphi_div_sin_theta(current) / r,
            ]
            # convert to Cartesian coordinates
            return surface.monitor.sph_2_car_field(temp[0], temp[1], temp[2], theta_obs, phi_obs)

        def potential_terms(current: Tuple[np.ndarray, ...], const: complex):
            """Assemble vector potential and its derivatives."""
            r_x_c = r_x_current(current)
            pot = [const * item * G for item in current]
            curl_pot = [const * item * dG_dr for item in r_x_c]
            grad_div_pot = grad_Gr_r_dot_current(current)
            grad_div_pot = [const * item for item in grad_div_pot]
            return pot, curl_pot, grad_div_pot

        # magnetic vector potential terms
        A, curl_A, grad_div_A = potential_terms(J, MU_0)

        # electric vector potential terms
        F, curl_F, grad_div_F = potential_terms(M, epsilon)

        # assemble the electric field components (Taflove 8.24, 8.27)
        e_x_integrand, e_y_integrand, e_z_integrand = (
            i_omega * (a + grad_div_a / (wavenumber**2)) - curl_f / epsilon
            for a, grad_div_a, curl_f in zip(A, grad_div_A, curl_F)
        )

        # assemble the magnetic field components (Taflove 8.25, 8.28)
        h_x_integrand, h_y_integrand, h_z_integrand = (
            i_omega * (f + grad_div_f / (wavenumber**2)) + curl_a / MU_0
            for f, grad_div_f, curl_a in zip(F, grad_div_F, curl_A)
        )

        # integrate over the surface
        e_x = self.integrate_2d(e_x_integrand, 1.0, pts[idx_u], pts[idx_v])
        e_y = self.integrate_2d(e_y_integrand, 1.0, pts[idx_u], pts[idx_v])
        e_z = self.integrate_2d(e_z_integrand, 1.0, pts[idx_u], pts[idx_v])
        h_x = self.integrate_2d(h_x_integrand, 1.0, pts[idx_u], pts[idx_v])
        h_y = self.integrate_2d(h_y_integrand, 1.0, pts[idx_u], pts[idx_v])
        h_z = self.integrate_2d(h_z_integrand, 1.0, pts[idx_u], pts[idx_v])

        # observation point in the original spherical system
        _, theta_obs, phi_obs = surface.monitor.car_2_sph(x, y, z)

        # convert fields to the original spherical system
        e_r, e_theta, e_phi = surface.monitor.car_2_sph_field(e_x, e_y, e_z, theta_obs, phi_obs)
        h_r, h_theta, h_phi = surface.monitor.car_2_sph_field(h_x, h_y, h_z, theta_obs, phi_obs)

        return [e_r, e_theta, e_phi, h_r, h_theta, h_phi]