diff --git a/.gitignore b/.gitignore index 141e29e9e7..998fa9705a 100644 --- a/.gitignore +++ b/.gitignore @@ -42,3 +42,4 @@ qdrant_storage/ .roo/ tests/core/sensor/sensor_count.xml tests/core/xml/test.xml +tmp/ diff --git a/examples/6-satellite-sensors/amsu-a.py b/examples/6-satellite-sensors/amsu-a.py new file mode 100644 index 0000000000..2542684519 --- /dev/null +++ b/examples/6-satellite-sensors/amsu-a.py @@ -0,0 +1,177 @@ +""" +Build an AMSU-A-like microwave sensor from a simple channel sheet. + +This is an example of how to set up simulations for a custom sensor. +""" + +from dataclasses import dataclass +import os + +import matplotlib.pyplot as plt +import numpy as np +import pyarts3 as pa + + +@dataclass(frozen=True) +class ChannelSpec: + """ A simple data structure to hold the channel specifications. """ + number: int + spec_text: str + reference_hz: float + bandwidth_hz: float + polarization: str + nedt_k: float + lo_offsets_hz: tuple[float, ...] = () + + + +# fmt: off +CHANNELS = [ + ChannelSpec(1, "$23.8$", 23.8e9, 270e6, "Iv", 0.30), + ChannelSpec(2, "$31.4$", 31.4e9, 180e6, "Iv", 0.30), + ChannelSpec(3, "$50.3$", 50.3e9, 180e6, "Iv", 0.40), + ChannelSpec(4, "$52.8$", 52.8e9, 400e6, "Iv", 0.25), + ChannelSpec(5, "$53.596 \\pm 0.115$", 53.596e9, 170e6, "Ih", 0.25, (115e6,)), + ChannelSpec(6, "$54.4$", 54.4e9, 400e6, "Ih", 0.25), + ChannelSpec(7, "$54.94$", 54.94e9, 400e6, "Iv", 0.25), + ChannelSpec(8, "$55.5$", 55.5e9, 330e6, "Ih", 0.25), + ChannelSpec(9, "$f_0 = 57.290344$", 57.290344e9, 330e6, "Ih", 0.25), + ChannelSpec(10, "$f_0 \\pm 0.217$", 57.290344e9, 78e6, "Ih", 0.40, (217e6,)), + ChannelSpec(11, "$f_0 \\pm 0.3222 \\pm 0.048$", 57.290344e9, 6e6, "Ih", 0.40, (322.2e6, 48e6)), + ChannelSpec(12, "$f_0 \\pm 0.3222 \\pm 0.022$", 57.290344e9, 16e6, "Ih", 0.60, (322.2e6, 22e6)), + ChannelSpec(13, "$f_0 \\pm 0.3222 \\pm 0.01$", 57.290344e9, 8e6, "Ih", 0.80, (322.2e6, 10e6)), + ChannelSpec(14, "$f_0 \\pm 0.3222 \\pm 0.0045$", 57.290344e9, 3e6, "Ih", 1.20, (322.2e6, 4.5e6)), + ChannelSpec(15, "$89 \\pm 1$", 89e9, 1000e6, "Iv", 0.50, (1000e6,)), +] + +CHANNEL_COUNT = len(CHANNELS) +SAMPLES_PER_CHANNEL = 11 +POS = [817e3, 0.0, 0.0] +LOS = [130.0, 0.0] +# fmt: on + + +def sensor_channels(channels, n): + """ Helper to turn the table above into channel responses """ + out = [] + + for ch in channels: + # The range is [f_ref, 0], [f_ref, inf] + x = pa.arts.SensorHeterodyneFrequencyRange(ch.reference_hz, [0, np.inf]) + + # Zero or more mixes, which add more ranges follow. + for lo in ch.lo_offsets_hz: + x.mix(lo) + + # The channel response is a simple boxcar with the specified bandwidth. + m = pa.arts.SensorBoxChannel(0, ch.bandwidth_hz / 2, n) + v = x.channel_response(m) + + out.append(v) + + return out + + +def extract_channels(spectral_rad, measurement_sensor): + """ Helper to extract internal simulation grid for channels """ + out = [] + + for ch in measurement_sensor: + v = np.where(ch.weight_matrix.dense.flatten())[0] + out.append([np.array([ch.f_grid[i//4] for i in v]), spectral_rad.flatten()[v]]) + + return out + + +# %% Download data and set up workspace +ws = pa.workspace.Workspace() +pa.data.download() + +# %% Use built-in spectroscopic data. +ws.abs_speciesSet(species=["O2-PWR98", "H2O-PWR98"]) +ws.ReadCatalogData() +ws.spectral_propmat_agendaAuto() + +# %% Set up a simple atmosphere and surface. +ws.surf_fieldPlanet(option="Earth") +ws.surf_field[pa.arts.SurfaceKey("t")] = 300.0 +ws.atm_fieldRead(toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1) + +# %% Use a geometric path. +ws.ray_path_observer_agendaSetGeometric() + +# %% Simple sensor setup: position, line-of-sight, polarization, and channels. +sensor = pa.arts.SensorBuilder(channels=sensor_channels(CHANNELS, SAMPLES_PER_CHANNEL)) +ws.measurement_sensor, ws.measurement_sensor_meta = sensor(POS, LOS) + +for i in range(len(ws.measurement_sensor)): + ws.measurement_sensor[i].normalize(CHANNELS[i].polarization) + +# %% Collect on a single grid to help plotting and speed-up calculations - this sacrifices Jacobian flexibilities +ws.measurement_sensor.collect() + +# %% Compute the sensor response in Planck units +ws.spectral_rad_transform_operatorSet(option="Tb") +ws.measurement_vecFromSensor(kernel="High Performance") + +# %% Compute the internal frequency grid response for plotting +ws.freq_grid = ws.measurement_sensor[0].f_grid +ws.spectral_rad_observer_agendaExecute(obs_pos=POS, obs_los=LOS) +ws.spectral_radApplyUnitFromSpectralRadiance() + +# %% Plotting +ch = extract_channels(ws.spectral_rad, ws.measurement_sensor) +colors = plt.get_cmap('tab20')(range(CHANNEL_COUNT)) +fig_internal, ax_internal = plt.subplots(figsize=(6, 4)) +for i in range(CHANNEL_COUNT): + ax_internal.plot(ch[i][0]/1e9, + ch[i][1], + "o:", + markersize=2.5, + linewidth=1.1, + label=f"Ch {CHANNELS[i].number}", + color=colors[i], + ) + +ax_internal.set_xlabel("Frequency [GHz]") +ax_internal.set_ylabel("Brightness temperature (internal grid) [K]") +ax_internal.set_title("AMSU-A-like channels") +ax_internal.grid(True, alpha=0.3) +ax_internal.legend(ncol=3, fontsize=8) + +fig_channels, ax_channels = plt.subplots(figsize=(6, 4)) +ax_channels.errorbar( + [spec.number for spec in CHANNELS], + ws.measurement_vec, + yerr=[spec.nedt_k for spec in CHANNELS], + fmt="o", + capsize=3, + linewidth=1.0, +) +ax_channels.set_xticks([spec.number for spec in CHANNELS]) +ax_channels.set_xlabel("AMSU-A channel number") +ax_channels.set_ylabel("Channel brightness [K]") +ax_channels.set_title("AMSU-A-like channelized measurements") +ax_channels.grid(True, alpha=0.3) + +for x, spec, y in zip([spec.number for spec in CHANNELS], CHANNELS, ws.measurement_vec): + ax_channels.annotate( + spec.spec_text, + (x, y), + textcoords="offset points", + xytext=(4, -2.5), + ha="left", + fontsize=6, + ) + +fig_internal.tight_layout() +fig_channels.tight_layout() + +if "ARTS_HEADLESS" not in os.environ: + plt.show() + +ref = [294.813805535829, 297.18778420039365, 283.26228661223786, 262.0365042504266, 244.4836135697355, 226.2694109308727, 216.11284609241537, + 209.72414936956548, 210.3336073749979, 219.42825163097098, 230.06432799428234, 241.47556633001957, 252.6256887712745, 261.4587517221863, 292.20874260781585] + +assert np.allclose(ws.measurement_vec, ref), \ + f"Mismatch to reference simulations.\nReference: {ref}\nSimulations: {ws.measurement_vec:B,}" diff --git a/python/src/pyarts3/plots/SensorObsel.py b/python/src/pyarts3/plots/SensorObsel.py new file mode 100644 index 0000000000..6103517d99 --- /dev/null +++ b/python/src/pyarts3/plots/SensorObsel.py @@ -0,0 +1,438 @@ +""" Plotting routine for SensorObsel """ + +import numpy +import matplotlib +import matplotlib.tri as mtri +import numpy as np +import pyarts3 as pyarts + +from .common import default_fig_ax, select_flat_ax + +__all__ = [ + 'plot', +] + + +def _line_plot_axes(fig, ax, nkeys): + return default_fig_ax(fig, ax, 1, nkeys, fig_kwargs={'figsize': (10 * nkeys, 10)}) + + +def _geometry_plot_axes(fig, ax, polar): + return default_fig_ax(fig, + ax, + ax_kwargs={"subplot_kw": {'polar': polar}}, + fig_kwargs={'figsize': (10, 8) if polar else (10, 6)}) + + +def _azimuth_for_plot(azi: np.ndarray, wrap_azimuth: bool) -> np.ndarray: + if wrap_azimuth: + return ((azi + 180.0) % 360.0) - 180.0 + + return azi + + +def _los_to_enu(los: np.ndarray) -> np.ndarray: + zen = np.deg2rad(los[..., 0]) + azi = np.deg2rad(los[..., 1]) + sin_zen = np.sin(zen) + return np.stack((sin_zen * np.sin(azi), + sin_zen * np.cos(azi), + np.cos(zen)), + axis=-1) + + +def _antenna_basis(bore_los: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + zen = np.deg2rad(bore_los[0]) + azi = np.deg2rad(bore_los[1]) + + cza = np.cos(zen) + sza = np.sin(zen) + caa = np.cos(azi) + saa = np.sin(azi) + + vertical = np.array([-cza * saa, -cza * caa, sza], dtype=float) + horizontal = np.array([caa, -saa, 0.0], dtype=float) + bore = np.array([sza * saa, sza * caa, cza], dtype=float) + return vertical, horizontal, bore + + +def _local_los_from_enu(enu: np.ndarray, bore_los: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + vertical, horizontal, bore = _antenna_basis(bore_los) + + local_vertical = enu @ vertical + local_horizontal = enu @ horizontal + local_bore = np.clip(enu @ bore, -1.0, 1.0) + + zen = np.rad2deg(np.arccos(local_bore)) + azi = np.rad2deg(np.arctan2(local_horizontal, -local_vertical)) + azi = np.where(zen > 0.0, azi, 0.0) + return zen, azi, local_horizontal, -local_vertical + + +def _reference_los(data: pyarts.arts.SensorObsel, ifreq: int | None) -> np.ndarray: + weights = np.asarray(data.weight_matrix, dtype=float) + iweights = weights[:, :, 0] + + if ifreq is None: + ref_weights = np.abs(iweights).sum(axis=1) + else: + nfreq = iweights.shape[1] + if ifreq < 0: + ifreq += nfreq + ref_weights = np.abs(iweights[:, ifreq]) + + los = np.asarray(data.poslos[:, 3:5], dtype=float) + enu = _los_to_enu(los) + + if np.any(ref_weights > 0.0): + ref = ref_weights @ enu + else: + ref = np.sum(enu, axis=0) + + ref_norm = np.linalg.norm(ref) + if ref_norm == 0.0: + return los[0] + + ref /= ref_norm + ref_zen = np.rad2deg(np.arccos(np.clip(ref[2], -1.0, 1.0))) + ref_azi = np.rad2deg(np.arctan2(ref[0], ref[1])) + return np.array([ref_zen, ref_azi], dtype=float) + + +def _geometry_coordinates(data: pyarts.arts.SensorObsel, + frame: str, + ifreq: int | None, + wrap_azimuth: bool) -> tuple[np.ndarray, np.ndarray, str, str, np.ndarray | None]: + poslos = np.asarray(data.poslos, dtype=float) + zen = poslos[:, 3] + azi = poslos[:, 4] + + if frame == "global": + plot_azi = _azimuth_for_plot(azi, wrap_azimuth) + return zen, plot_azi, "Azimuth angle [deg]", "Zenith angle [deg]", None + + if frame != "local": + raise ValueError("frame must be 'local' or 'global'") + + bore_los = _reference_los(data, ifreq) + enu = _los_to_enu(poslos[:, 3:5]) + local_zen, local_azi, xoff, yoff = _local_los_from_enu(enu, bore_los) + + radial = local_zen + theta = _azimuth_for_plot(local_azi, wrap_azimuth) + _, _, bore = _antenna_basis(bore_los) + projected_bore = np.clip(enu @ bore, -1.0, 1.0) + xdeg = np.rad2deg(np.arctan2(xoff, projected_bore)) + ydeg = np.rad2deg(np.arctan2(yoff, projected_bore)) + + return radial, theta, "Cross-track offset [deg]", "Along-track offset [deg]", np.stack((xdeg, ydeg), axis=-1) + + +def _geometry_values(data: pyarts.arts.SensorObsel, + pol: pyarts.arts.Stokvec, + ifreq: int | None) -> tuple[np.ndarray, str]: + if ifreq is None: + values = np.asarray(data.weight_matrix.reduce(pol, + along_poslos=False, + along_freq=True), + dtype=float) + return values, "Integrated weight" + + weights = np.asarray(data.weight_matrix, dtype=float) + nfreq = weights.shape[1] + if ifreq < -nfreq or ifreq >= nfreq: + raise IndexError(f"ifreq {ifreq} out of range for {nfreq} frequencies") + + if ifreq < 0: + ifreq += nfreq + + pol_vec = np.asarray(pol, dtype=float) + values = weights[:, ifreq, :] @ pol_vec + return values, f"Weight at f[{ifreq}] = {float(data.f_grid[ifreq]):g}" + + +def _geometry_plot_type(point_spread: bool, type: str | None) -> str: + if type is None: + return "points" + + plot_type = str(type).lower() + if plot_type == "scatter": + plot_type = "points" + + valid_types = ("points", "cloud", "contour") + if plot_type not in valid_types: + choices = ", ".join(repr(choice) for choice in valid_types) + raise ValueError(f"type must be one of {choices}") + + if plot_type != "points" and not point_spread: + raise ValueError(f"type={type!r} requires point_spread=True") + + return plot_type + + +def _triangulated_point_spread(local_xy: np.ndarray, + values: np.ndarray) -> tuple[mtri.Triangulation, np.ndarray]: + x = np.asarray(local_xy[:, 0], dtype=float) + y = np.asarray(local_xy[:, 1], dtype=float) + values = np.asarray(values, dtype=float) + + finite = np.isfinite(x) & np.isfinite(y) & np.isfinite(values) + x = x[finite] + y = y[finite] + values = values[finite] + + if x.size < 3: + raise ValueError("Interpolated point-spread plots require at least 3 geometry samples") + + coords = np.stack((x, y), axis=-1) + unique_coords, inverse = np.unique(coords, axis=0, return_inverse=True) + if unique_coords.shape[0] != coords.shape[0]: + unique_values = np.zeros(unique_coords.shape[0], dtype=float) + counts = np.zeros(unique_coords.shape[0], dtype=int) + np.add.at(unique_values, inverse, values) + np.add.at(counts, inverse, 1) + x = unique_coords[:, 0] + y = unique_coords[:, 1] + values = unique_values / counts + + if x.size < 3: + raise ValueError("Interpolated point-spread plots require at least 3 unique geometry samples") + + triangulation = mtri.Triangulation(x, y) + if triangulation.triangles.size == 0: + raise ValueError("Interpolated point-spread plots require non-collinear geometry samples") + + return triangulation, values + + +def _contour_levels(values: np.ndarray, levels: int | np.ndarray) -> np.ndarray | int: + if np.ndim(levels) != 0: + return levels + + nlevels = max(int(levels), 2) + vmin = float(np.min(values)) + vmax = float(np.max(values)) + if np.isclose(vmin, vmax): + delta = 1e-6 if vmin == 0.0 else abs(vmin) * 1e-6 + return np.array((vmin - delta, vmax + delta), dtype=float) + + return np.linspace(vmin, vmax, nlevels) + + +def _plot_geometry(data: pyarts.arts.SensorObsel, + *, + fig: matplotlib.figure.Figure | None, + ax: matplotlib.axes.Axes | list[matplotlib.axes.Axes] | numpy.ndarray[matplotlib.axes.Axes] | None, + pol: pyarts.arts.Stokvec, + polar: bool, + point_spread: bool, + type: str | None, + ifreq: int | None, + colorbar: bool, + frame: str, + wrap_azimuth: bool, + **kwargs) -> tuple[matplotlib.figure.Figure, matplotlib.axes.Axes | list[matplotlib.axes.Axes] | numpy.ndarray[matplotlib.axes.Axes]]: + fig, ax = _geometry_plot_axes(fig, ax, polar) + axis = select_flat_ax(ax, 0) + + radial, theta, xlabel, ylabel, local_xy = _geometry_coordinates(data, + frame, + ifreq, + wrap_azimuth) + plot_type = _geometry_plot_type(point_spread, type) + + if plot_type == "points": + default_kwargs = {'marker': 'o', 'linestyle': 'None'} + for key, value in default_kwargs.items(): + kwargs.setdefault(key, value) + + artist = None + if point_spread: + values, label = _geometry_values(data, pol, ifreq) + kwargs.setdefault('cmap', 'viridis') + + if plot_type == "cloud" or plot_type == "contour": + if polar: + raise ValueError("Interpolated point-spread types require polar=False") + if frame != "local" or local_xy is None: + raise ValueError("Interpolated point-spread types require frame='local'") + + triangulation, values = _triangulated_point_spread(local_xy, values) + if plot_type == "cloud": + kwargs.setdefault('shading', 'gouraud') + artist = axis.tripcolor(triangulation, values, **kwargs) + else: + levels = _contour_levels(values, kwargs.pop('levels', 32)) + artist = axis.tricontourf(triangulation, values, levels=levels, **kwargs) + elif polar: + artist = axis.scatter(np.deg2rad(theta), radial, c=values, **kwargs) + elif frame == "local": + artist = axis.scatter(local_xy[:, 0], local_xy[:, 1], c=values, **kwargs) + else: + artist = axis.scatter(theta, radial, c=values, **kwargs) + if colorbar: + fig.colorbar(artist, ax=axis, label=label) + else: + if polar: + artist = axis.plot(np.deg2rad(theta), radial, **kwargs) + elif frame == "local": + artist = axis.plot(local_xy[:, 0], local_xy[:, 1], **kwargs) + else: + artist = axis.plot(theta, radial, **kwargs) + + if polar: + axis.set_theta_zero_location("N") + axis.set_theta_direction(-1) + axis.set_rlabel_position(225) + axis.set_ylim(0.0, np.max(radial) if radial.size else 1.0) + else: + axis.set_xlabel(xlabel) + axis.set_ylabel(ylabel) + if frame == "local": + axis.set_aspect('equal', adjustable='box') + + return fig, ax + + +def plot(data: pyarts.arts.SensorObsel, + *, + fig: matplotlib.figure.Figure | None = None, + ax: matplotlib.axes.Axes | list[matplotlib.axes.Axes] | numpy.ndarray[matplotlib.axes.Axes] | None = None, + keys: str | list = "f", + pol: str | pyarts.arts.Stokvec = "I", + polar: bool = False, + point_spread: bool = False, + type: str | None = None, + ifreq: int | None = None, + colorbar: bool = True, + frame: str = "global", + wrap_azimuth: bool = False, + **kwargs) -> tuple[matplotlib.figure.Figure, matplotlib.axes.Axes | list[matplotlib.axes.Axes] | numpy.ndarray[matplotlib.axes.Axes]]: + """Plot a sensor observational element. + + By default, this follows the same line-plot convention as + :mod:`pyarts3.plots.ArrayOfSensorObsel`: the selected quantity is reduced + over the non-plotted axis and drawn against either frequency or one sensor + geometry component. + + Setting ``polar=True`` switches to a geometry view. By default, geometry + plots use the stored ``SensorObsel.poslos`` zenith and azimuth values so + the plotted coordinates match the observation element directly. + Setting ``point_spread=True`` colors the geometry samples by their reduced + weights, integrated over all frequencies unless ``ifreq`` selects one + frequency column. The geometry ``type`` selects whether those weighted + samples appear as discrete points, a smooth triangulated cloud, or filled + triangulated contours. + + .. rubric:: Example + + .. plot:: + :include-source: + + import pyarts3 as pyarts + import numpy as np + + ant = pyarts.arts.SensorGaussianAiryAntenna(np.linspace(0, 1, 11), + np.linspace(0, 330, 12), + 0.3) + ch = pyarts.arts.SensorDiracChannel(100e9) + obsel = ant(ch, [1, 0, 0], [90, 0]) + + pyarts.plots.SensorObsel.plot(obsel, polar=True, point_spread=True) + pyarts.plots.SensorObsel.plot(obsel, point_spread=True, wrap_azimuth=True) + pyarts.plots.SensorObsel.plot(obsel, point_spread=True, frame="local") + pyarts.plots.SensorObsel.plot(obsel, point_spread=True, frame="local", type="cloud") + pyarts.plots.SensorObsel.plot(obsel, point_spread=True, frame="local", type="contour") + + Parameters + ---------- + data : ~pyarts3.arts.SensorObsel + A single sensor observation element. + fig : ~matplotlib.figure.Figure, optional + The matplotlib figure to draw on. Defaults to None for new figure. + ax : ~matplotlib.axes.Axes | list[~matplotlib.axes.Axes] | ~numpy.ndarray[~matplotlib.axes.Axes] | None, optional + The matplotlib axes to draw on. Defaults to None for new axes. + keys : str | list + The keys to use for line plotting. Options are in :class:`~pyarts3.arts.SensorKeyType`. + Ignored when ``polar`` or ``point_spread`` is enabled. + pol : str | ~pyarts3.arts.Stokvec + The polarization to use for plotting. Defaults to ``"I"``. + polar : bool, optional + If True, plot the observation geometry in polar coordinates. Defaults to False. + point_spread : bool, optional + If True, color the observation geometry by weight values. Defaults to False. + type : str | None, optional + Geometry rendering type. ``"points"`` shows discrete samples, + ``"cloud"`` draws a smooth triangulated color cloud, and + ``"contour"`` draws filled triangulated contours. Interpolated types + require ``point_spread=True``, ``frame="local"``, and ``polar=False``. + Defaults to ``None``, which behaves like ``"points"``. + ifreq : int | None, optional + Frequency index to use for point-spread values. If None, integrates over frequency. + colorbar : bool, optional + If True, add a colorbar to point-spread plots. Defaults to True. + frame : str, optional + Geometry frame for ``polar`` and ``point_spread`` plots. ``"global"`` + uses the stored zenith/azimuth coordinates from ``SensorObsel.poslos``, + while ``"local"`` uses a bore-centered frame inferred from the weighted + LOS samples. Defaults to ``"global"``. + wrap_azimuth : bool, optional + If True, wrap azimuths to ``[-180, 180)`` before geometry plotting. + Defaults to False. + **kwargs : keyword arguments + Additional keyword arguments to pass to the plotting functions. + + Returns + ------- + fig : + As input if input. Otherwise the created Figure. + ax : + As input if input. Otherwise the created Axes. + """ + + pol = pyarts.arts.Stokvec(pol) + + if type is not None and not (polar or point_spread): + _geometry_plot_type(point_spread, type) + + if polar or point_spread: + return _plot_geometry(data, + fig=fig, + ax=ax, + pol=pol, + polar=polar, + point_spread=point_spread, + type=type, + ifreq=ifreq, + colorbar=colorbar, + frame=frame, + wrap_azimuth=wrap_azimuth, + **kwargs) + + keys = [keys] if isinstance(keys, str) else keys + nkeys = len(keys) + if nkeys == 0: + return fig, ax + + fig, ax = _line_plot_axes(fig, ax, nkeys) + + key_map = { + pyarts.arts.SensorKeyType.f: None, + pyarts.arts.SensorKeyType.alt: 0, + pyarts.arts.SensorKeyType.lat: 1, + pyarts.arts.SensorKeyType.lon: 2, + pyarts.arts.SensorKeyType.zen: 3, + pyarts.arts.SensorKeyType.azi: 4, + } + + for isub, raw_key in enumerate(keys): + key = pyarts.arts.SensorKeyType(raw_key) + idx = key_map[key] + + values = data.weight_matrix.reduce(pol, + along_poslos=idx is None, + along_freq=idx is not None) + x = data.f_grid if idx is None else data.poslos[:, idx] + select_flat_ax(ax, isub).plot(x, values, **kwargs) + + return fig, ax \ No newline at end of file diff --git a/python/src/pyarts3/plots/__init__.py b/python/src/pyarts3/plots/__init__.py index 6b1e5b4287..73c7648e98 100644 --- a/python/src/pyarts3/plots/__init__.py +++ b/python/src/pyarts3/plots/__init__.py @@ -54,6 +54,7 @@ from . import SubsurfaceField from . import Sun from . import SurfaceField +from . import SensorObsel from . import Vector from . import ZenGrid diff --git a/src/core/sensor/CMakeLists.txt b/src/core/sensor/CMakeLists.txt index e5f48ad45e..7e0b2f8af3 100644 --- a/src/core/sensor/CMakeLists.txt +++ b/src/core/sensor/CMakeLists.txt @@ -1,3 +1,13 @@ -add_library(sensor STATIC obsel.cpp sensor_meta_info.cpp) +add_library(sensor STATIC + antenna_pattern.cpp + obsel.cpp + sensor_meta_info.cpp + sensor_builder.cpp + frequency_bandpass_filters.cpp + frequency_channel_selection.cpp + frequency_range_selection.cpp +) + target_link_libraries(sensor PUBLIC matpack rtepack arts_enum_options) + target_include_directories(sensor PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/core/sensor/antenna_pattern.cpp b/src/core/sensor/antenna_pattern.cpp new file mode 100644 index 0000000000..d25ac164cd --- /dev/null +++ b/src/core/sensor/antenna_pattern.cpp @@ -0,0 +1,476 @@ +#include "antenna_pattern.h" + +#include +#include +#include + +#include +#include +#include +#include + +namespace sensor { + +namespace { +constexpr Numeric gaussian_airy_hwhm_factor = + 3.8317059702075123156 / Constant::pi; + +struct AntennaBasis { + Vector3 v; + Vector3 h; + Vector3 k; +}; + +struct AntennaGeometrySample { + Size sample_index{}; + Numeric ant_zen{}; + Numeric ant_azi{}; + bool has_response{}; + bool is_representative{}; +}; + +struct AntennaGeometryLayout { + std::shared_ptr poslos_grid; + std::vector samples; +}; + +constexpr Numeric azimuth_degenerate_tol = 1e-12; + +[[nodiscard]] AntennaBasis antenna_basis(Vector2 bore_los) { + using Conversion::cosd, Conversion::sind; + + const Numeric cza = cosd(bore_los[0]); + const Numeric sza = sind(bore_los[0]); + const Numeric caa = cosd(bore_los[1]); + const Numeric saa = sind(bore_los[1]); + + return { + .v = {-cza * saa, -cza * caa, sza}, + .h = {caa, -saa, 0.0}, + .k = {sza * saa, sza * caa, cza}, + }; +} + +[[nodiscard]] Vector3 antenna_frame_los(Vector2 local_los) { + using Conversion::cosd, Conversion::sind; + + const Numeric cza = cosd(local_los[0]); + const Numeric sza = sind(local_los[0]); + const Numeric caa = cosd(local_los[1]); + const Numeric saa = sind(local_los[1]); + + return {-sza * caa, sza * saa, cza}; +} + +[[nodiscard]] bool antenna_azimuth_is_degenerate(const Vector3& local) { + return std::hypot(local[0], local[1]) <= azimuth_degenerate_tol; +} + +template +[[nodiscard]] AntennaGeometryLayout make_antenna_geometry_layout( + const ZenGridT& zen_grid, + const AziGridT& azi_grid, + const Vector3& pos, + const Vector2& bore_los) { + const auto basis = antenna_basis(bore_los); + + std::vector poslos; + poslos.reserve(zen_grid.size() * azi_grid.size()); + + std::vector samples; + samples.reserve(zen_grid.size() * azi_grid.size()); + + using Conversion::atan2d; + + for (Numeric izen : zen_grid) { + Size degenerate_sample_index = std::numeric_limits::max(); + + for (Numeric iazi : azi_grid) { + const Vector3 local = antenna_frame_los({izen, iazi}); + const bool degenerate = antenna_azimuth_is_degenerate(local); + + Size sample_index = degenerate_sample_index; + bool is_representative = false; + if (sample_index == std::numeric_limits::max()) { + const Vector3 enu = normalized(local[0] * basis.v + local[1] * basis.h + + local[2] * basis.k); + sample_index = poslos.size(); + poslos.push_back({.pos = pos, .los = enu2los(enu)}); + is_representative = true; + + if (degenerate) degenerate_sample_index = sample_index; + } + + auto& sample = samples.emplace_back(AntennaGeometrySample{ + .sample_index = sample_index, + .has_response = local[2] > 0.0, + .is_representative = is_representative}); + + if (sample.has_response) { + sample.ant_zen = atan2d(local[0], local[2]); + sample.ant_azi = atan2d(local[1], local[2]); + } + } + } + + auto poslos_grid = std::make_shared(poslos.size()); + for (Size isample = 0; isample < poslos.size(); ++isample) { + (*poslos_grid)[isample] = poslos[isample]; + } + + return {std::move(poslos_grid), std::move(samples)}; +} + +[[nodiscard]] AntennaPatternField make_gaussian_field(ZenGrid zen_grid, + AziGrid azi_grid, + Numeric zenith_std, + Numeric azimuth_std, + Stokvec weight) { + ARTS_USER_ERROR_IF(zenith_std <= 0.0, + "Gaussian antenna zenith_std must be positive") + ARTS_USER_ERROR_IF(azimuth_std <= 0.0, + "Gaussian antenna azimuth_std must be positive") + + AntennaPatternField out{ + .data_name = "gaussian"s, + .data = StokvecMatrix(zen_grid.size(), azi_grid.size()), + .grid_names = {"zenith"s, "azimuth"s}, + .grids = {std::move(zen_grid), std::move(azi_grid)}, + }; + + using Conversion::atan2d; + + for (Size izen = 0; izen < out.grid<0>().size(); ++izen) { + for (Size iazi = 0; iazi < out.grid<1>().size(); ++iazi) { + const Vector3 local = + antenna_frame_los({out.grid<0>()[izen], out.grid<1>()[iazi]}); + + if (local[2] <= 0.0) { + out[izen, iazi] = {0.0, 0.0, 0.0, 0.0}; + continue; + } + + const Numeric ant_zen = atan2d(local[0], local[2]); + const Numeric ant_azi = atan2d(local[1], local[2]); + const Numeric exponent = + -0.5 * ((ant_zen / zenith_std) * (ant_zen / zenith_std) + + (ant_azi / azimuth_std) * (ant_azi / azimuth_std)); + + out[izen, iazi] = std::exp(exponent) * weight; + } + } + + return out; +} + +[[nodiscard]] Numeric gaussian_airy_std(Numeric frequency, + Numeric aperture_diameter) { + const Numeric wavelength = Constant::speed_of_light / frequency; + const Numeric hwhm_deg = Conversion::rad2deg(gaussian_airy_hwhm_factor * + wavelength / aperture_diameter); + + return Conversion::hwhm2std(hwhm_deg); +} + +[[nodiscard]] Numeric gaussian_airy_response(Numeric ant_zen, + Numeric ant_azi, + Numeric inv_std_sq) { + return std::exp(-0.5 * (ant_zen * ant_zen + ant_azi * ant_azi) * inv_std_sq); +} + +std::shared_ptr make_single_poslos_grid( + const Vector3& pos, const Vector2& los) { + return std::make_shared(1, PosLos{.pos = pos, .los = los}); +} +} // namespace + +PencilBeamAntenna::PencilBeamAntenna(Stokvec weight) : weight(weight) {} + +Obsel PencilBeamAntenna::operator()(const Channel& channel, + const Vector3& pos, + const Vector2& bore_los) const { + const auto& channel_weights = channel.weights(); + const auto freq_grid = + std::make_shared(channel.freq_grid()); + SparseStokvecMatrix weight_matrix(1, channel_weights.size()); + + if (not weight.is_zero()) { + for (Size ifreq = 0; ifreq < channel_weights.size(); ++ifreq) { + if (channel_weights[ifreq] == 0.0) continue; + weight_matrix[0, ifreq] = channel_weights[ifreq] * weight; + } + } + + return {freq_grid, + make_single_poslos_grid(pos, bore_los), + std::move(weight_matrix)}; +} + +std::shared_ptr PencilBeamAntenna::clone() const { + return std::make_shared(*this); +} + +Obsel GriddedAntennaPattern::operator()(const Channel& channel, + const Vector3& pos, + const Vector2& bore_los) const { + ARTS_USER_ERROR_IF( + not data.ok(), + "SensorGriddedAntennaPattern data shape does not match its grids") + + const auto& zen_grid = data.grid<0>(); + const auto& azi_grid = data.grid<1>(); + const auto& channel_weights = channel.weights(); + const auto freq_grid = + std::make_shared(channel.freq_grid()); + + auto geometry = + make_antenna_geometry_layout(zen_grid, azi_grid, pos, bore_los); + SparseStokvecMatrix weight_matrix(geometry.poslos_grid->size(), + channel_weights.size()); + + Size igrid = 0; + + for (Size izen = 0; izen < zen_grid.size(); ++izen) { + for (Size iazi = 0; iazi < azi_grid.size(); ++iazi) { + const auto& antenna_weight = data[izen, iazi]; + if (not antenna_weight.is_zero()) { + const Size sample_index = geometry.samples[igrid].sample_index; + for (Size ifreq = 0; ifreq < channel_weights.size(); ++ifreq) { + if (channel_weights[ifreq] == 0.0) continue; + weight_matrix[sample_index, ifreq] += + channel_weights[ifreq] * antenna_weight; + } + } + + ++igrid; + } + } + + return {freq_grid, std::move(geometry.poslos_grid), std::move(weight_matrix)}; +} + +std::shared_ptr GriddedAntennaPattern::clone() const { + return std::make_shared(*this); +} + +GaussianAntenna::GaussianAntenna(ZenGrid zen_grid, + AziGrid azi_grid, + Numeric zenith_std, + Numeric azimuth_std, + Stokvec weight) { + data = make_gaussian_field(std::move(zen_grid), + std::move(azi_grid), + zenith_std, + azimuth_std, + weight); +} + +GaussianAntenna::GaussianAntenna(ZenGrid zen_grid, + AziGrid azi_grid, + Numeric std, + Stokvec weight) + : GaussianAntenna( + std::move(zen_grid), std::move(azi_grid), std, std, weight) {} + +std::shared_ptr GaussianAntenna::clone() const { + return std::make_shared(*this); +} + +GaussianAiryAntenna::GaussianAiryAntenna(ZenGrid zen_grid, + AziGrid azi_grid, + Numeric aperture_diameter, + Stokvec weight) + : zen_grid(std::move(zen_grid)), + azi_grid(std::move(azi_grid)), + aperture_diameter(aperture_diameter), + weight(weight) {} + +Obsel GaussianAiryAntenna::operator()(const Channel& channel, + const Vector3& pos, + const Vector2& bore_los) const { + ARTS_USER_ERROR_IF(aperture_diameter <= 0.0, + "Gaussian Airy antenna aperture_diameter must be positive") + + const auto& channel_weights = channel.weights(); + const auto freq_grid = + std::make_shared(channel.freq_grid()); + + ARTS_USER_ERROR_IF( + freq_grid->front() <= 0.0, + "Gaussian Airy antenna requires strictly positive channel frequencies because SensorBuilder only provides the channel frequency grid") + + auto geometry = + make_antenna_geometry_layout(zen_grid, azi_grid, pos, bore_los); + SparseStokvecMatrix weight_matrix(geometry.poslos_grid->size(), + channel_weights.size()); + + if (not weight.is_zero()) { + Vector inv_std_sq(channel_weights.size(), 0.0); + Vector normalization(channel_weights.size(), 0.0); + Vector retained_scale(channel_weights.size(), 0.0); + + for (Size ifreq = 0; ifreq < channel_weights.size(); ++ifreq) { + if (channel_weights[ifreq] == 0.0) continue; + + const Numeric airy_std = + gaussian_airy_std((*freq_grid)[ifreq], aperture_diameter); + inv_std_sq[ifreq] = 1.0 / (airy_std * airy_std); + + Numeric retained_mass = 0.0; + for (const auto& sample : geometry.samples) { + if (not sample.has_response) continue; + normalization[ifreq] += gaussian_airy_response( + sample.ant_zen, sample.ant_azi, inv_std_sq[ifreq]); + } + + if (normalization[ifreq] == 0.0) continue; + + for (const auto& sample : geometry.samples) { + if (not sample.has_response) continue; + if (not sample.is_representative) continue; + + const Numeric single_weight = + channel_weights[ifreq] * gaussian_airy_response( + sample.ant_zen, + sample.ant_azi, + inv_std_sq[ifreq]) / + normalization[ifreq]; + + retained_mass += single_weight; + } + + if (retained_mass > 0.0) { + retained_scale[ifreq] = channel_weights[ifreq] / retained_mass; + } + } + + // Keep sparse insertions row-major so they append instead of shifting. + for (const auto& sample : geometry.samples) { + if (not sample.has_response) continue; + if (not sample.is_representative) continue; + + for (Size ifreq = 0; ifreq < channel_weights.size(); ++ifreq) { + if (normalization[ifreq] == 0.0) continue; + if (retained_scale[ifreq] == 0.0) continue; + + const Numeric single_weight = + channel_weights[ifreq] * gaussian_airy_response( + sample.ant_zen, + sample.ant_azi, + inv_std_sq[ifreq]) / + normalization[ifreq]; + if (single_weight == 0.0) continue; + + weight_matrix[sample.sample_index, ifreq] += + (retained_scale[ifreq] * single_weight) * weight; + } + } + } + + return {freq_grid, std::move(geometry.poslos_grid), std::move(weight_matrix)}; +} + +std::shared_ptr GaussianAiryAntenna::clone() const { + return std::make_shared(*this); +} + +static_assert(AntennaPatternSelection); +static_assert(AntennaPatternSelection); +static_assert(AntennaPatternSelection); +static_assert(AntennaPatternSelection); +static_assert(AntennaPatternSelection); +} // namespace sensor + +void xml_io_stream::write( + std::ostream& os, + const sensor::GriddedAntennaPattern& n, + bofstream* pbofs, + std::string_view name) { + XMLTag tag(xml_io_stream_name_v, "name", name); + tag.write_to_stream(os); + + xml_write_to_stream(os, n.data, pbofs); + + tag.write_to_end_stream(os); +} + +void xml_io_stream::read( + std::istream& is, sensor::GriddedAntennaPattern& n, bifstream* pbifs) { + XMLTag tag{}; + tag.read_from_stream(is); + tag.check_name(xml_io_stream_name_v); + + xml_read_from_stream(is, n.data, pbifs); + + tag.read_from_stream(is); + tag.check_end_name(xml_io_stream_name_v); +} + +void xml_io_stream::write( + std::ostream& os, + const sensor::GaussianAiryAntenna& n, + bofstream* pbofs, + std::string_view name) { + XMLTag tag(xml_io_stream_name_v, "name", name); + tag.write_to_stream(os); + xml_write_to_stream(os, n.aperture_diameter, pbofs); + xml_write_to_stream(os, n.zen_grid, pbofs); + xml_write_to_stream(os, n.azi_grid, pbofs); + xml_write_to_stream(os, n.weight, pbofs); + tag.write_to_end_stream(os); +} + +void xml_io_stream::read( + std::istream& is, sensor::GaussianAiryAntenna& n, bifstream* pbifs) { + XMLTag tag{}; + tag.read_from_stream(is); + tag.check_name(xml_io_stream_name_v); + + xml_read_from_stream(is, n.aperture_diameter, pbifs); + xml_read_from_stream(is, n.zen_grid, pbifs); + xml_read_from_stream(is, n.azi_grid, pbifs); + xml_read_from_stream(is, n.weight, pbifs); + + tag.read_from_stream(is); + tag.check_end_name(xml_io_stream_name_v); +} + +void xml_io_stream::write( + std::ostream& os, + const sensor::PencilBeamAntenna& n, + bofstream* pbofs, + std::string_view name) { + XMLTag tag(xml_io_stream_name_v, "name", name); + tag.write_to_stream(os); + xml_write_to_stream(os, n.weight, pbofs); + tag.write_to_end_stream(os); +} + +void xml_io_stream::read( + std::istream& is, sensor::PencilBeamAntenna& n, bifstream* pbifs) { + XMLTag tag{}; + tag.read_from_stream(is); + tag.check_name(xml_io_stream_name_v); + xml_read_from_stream(is, n.weight, pbifs); + tag.read_from_stream(is); + tag.check_end_name(xml_io_stream_name_v); +} + +void xml_io_stream::write(std::ostream& os, + const sensor::AntennaPattern&, + bofstream*, + std::string_view name) { + XMLTag tag(xml_io_stream_name_v, "name", name); + tag.write_to_stream(os); + tag.write_to_end_stream(os); +} + +void xml_io_stream::read(std::istream& is, + sensor::AntennaPattern&, + bifstream*) { + XMLTag tag{}; + tag.read_from_stream(is); + tag.check_name(xml_io_stream_name_v); + tag.read_from_stream(is); + tag.check_end_name(xml_io_stream_name_v); +} diff --git a/src/core/sensor/antenna_pattern.h b/src/core/sensor/antenna_pattern.h new file mode 100644 index 0000000000..3edf62cf3e --- /dev/null +++ b/src/core/sensor/antenna_pattern.h @@ -0,0 +1,314 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include + +#include "frequency_channel_selection.h" +#include "obsel.h" + +namespace sensor { +//! A 2D angular antenna pattern on local zenith and azimuth offsets. +struct AntennaPattern; + +//! A 2D gridded antenna pattern on local zenith and azimuth offsets. +struct GriddedAntennaPattern; + +//! A 1x1 antenna pattern that samples only the bore line of sight. +struct PencilBeamAntenna; + +//! A 2D Gaussian antenna pattern on local zenith and azimuth offsets. +struct GaussianAntenna; + +//! A Gaussianized Airy antenna pattern with frequency-dependent width. +struct GaussianAiryAntenna; + +//! 2D gridded field of antenna weights on local zenith and azimuth offsets. +using AntennaPatternField = matpack::gridded_data_t; + +struct AntennaPattern { + virtual ~AntennaPattern() = default; + AntennaPattern() = default; + AntennaPattern(const AntennaPattern&) = default; + AntennaPattern& operator=(const AntennaPattern&) = default; + AntennaPattern(AntennaPattern&&) = default; + AntennaPattern& operator=(AntennaPattern&&) = default; + + //! Builds one sensor obsel for a channel at a sensor position and bore LOS. + [[nodiscard]] virtual Obsel operator()(const Channel& channel, + const Vector3& pos, + const Vector2& bore_los) const = 0; + + //! Creates an owning copy preserving the dynamic antenna type. + [[nodiscard]] virtual std::shared_ptr clone() const = 0; +}; + +struct GriddedAntennaPattern : AntennaPattern { + AntennaPatternField data; // center at [0, 0] + + [[nodiscard]] Obsel operator()(const Channel& channel, + const Vector3& pos, + const Vector2& bore_los) const override; + + [[nodiscard]] std::shared_ptr clone() const override; +}; + +struct PencilBeamAntenna final : AntennaPattern { + Stokvec weight{1.0, 0.0, 0.0, 0.0}; + + PencilBeamAntenna(Stokvec weight = {1.0, 0.0, 0.0, 0.0}); + + [[nodiscard]] Obsel operator()(const Channel& channel, + const Vector3& pos, + const Vector2& bore_los) const override; + + [[nodiscard]] std::shared_ptr clone() const override; +}; + +struct GaussianAntenna final : GriddedAntennaPattern { + GaussianAntenna() = default; + GaussianAntenna(const GaussianAntenna&) = default; + GaussianAntenna& operator=(const GaussianAntenna&) = default; + GaussianAntenna(GaussianAntenna&&) = default; + GaussianAntenna& operator=(GaussianAntenna&&) = default; + + GaussianAntenna(ZenGrid zen_grid, + AziGrid azi_grid, + Numeric zenith_std, + Numeric azimuth_std, + Stokvec weight = {1.0, 0.0, 0.0, 0.0}); + + GaussianAntenna(ZenGrid zen_grid, + AziGrid azi_grid, + Numeric std, + Stokvec weight = {1.0, 0.0, 0.0, 0.0}); + + [[nodiscard]] std::shared_ptr clone() const override; +}; + +struct GaussianAiryAntenna final : AntennaPattern { + ZenGrid zen_grid; + AziGrid azi_grid; + Numeric aperture_diameter; + Stokvec weight; + + GaussianAiryAntenna() = default; + GaussianAiryAntenna(const GaussianAiryAntenna&) = default; + GaussianAiryAntenna& operator=(const GaussianAiryAntenna&) = default; + GaussianAiryAntenna(GaussianAiryAntenna&&) = default; + GaussianAiryAntenna& operator=(GaussianAiryAntenna&&) = default; + + GaussianAiryAntenna(ZenGrid zen_grid, + AziGrid azi_grid, + Numeric aperture_diameter, + Stokvec weight = {1.0, 0.0, 0.0, 0.0}); + + [[nodiscard]] Obsel operator()(const Channel& channel, + const Vector3& pos, + const Vector2& bore_los) const override; + + [[nodiscard]] std::shared_ptr clone() const override; +}; + +//! Concept for types that can build a sensor obsel from channel geometry. +template +concept AntennaPatternSelection = requires(const T& antenna, + const Channel& channel, + const Vector3& pos, + const Vector2& bore_los) { + { antenna(channel, pos, bore_los) } -> std::same_as; +}; +} // namespace sensor + +// AntennaPattern format tags and XML I/O + +template <> +struct std::formatter { + format_tags tags; + + [[nodiscard]] constexpr auto& inner_fmt() { return *this; } + [[nodiscard]] constexpr auto& inner_fmt() const { return *this; } + + constexpr std::format_parse_context::iterator parse( + std::format_parse_context& ctx) { + return parse_format_tags(tags, ctx); + } + + template + FmtContext::iterator format(const sensor::AntennaPattern&, + FmtContext& ctx) const { + return tags.format(ctx, "SensorAntennaPattern"sv); + } +}; + +template <> +struct xml_io_stream_name { + static constexpr std::string_view name = "SensorAntennaPattern"; +}; + +template <> +struct xml_io_stream { + static constexpr std::string_view type_name = + xml_io_stream_name_v; + + static void write(std::ostream& os, + const sensor::AntennaPattern& n, + bofstream* pbofs = nullptr, + std::string_view name = ""sv); + + static void read(std::istream& is, + sensor::AntennaPattern& n, + bifstream* pbifs = nullptr); +}; + +// GriddedAntennaPattern format tags and XML I/O + +template <> +struct std::formatter { + format_tags tags; + + [[nodiscard]] constexpr auto& inner_fmt() { return *this; } + [[nodiscard]] constexpr auto& inner_fmt() const { return *this; } + + constexpr std::format_parse_context::iterator parse( + std::format_parse_context& ctx) { + return parse_format_tags(tags, ctx); + } + + template + FmtContext::iterator format(const sensor::GriddedAntennaPattern& v, + FmtContext& ctx) const { + return tags.format(ctx, v.data); + } +}; + +template <> +struct xml_io_stream_name { + static constexpr std::string_view name = "SensorGriddedAntennaPattern"; +}; + +template <> +struct xml_io_stream { + static constexpr std::string_view type_name = + xml_io_stream_name_v; + + static void write(std::ostream& os, + const sensor::GriddedAntennaPattern& n, + bofstream* pbofs = nullptr, + std::string_view name = ""sv); + + static void read(std::istream& is, + sensor::GriddedAntennaPattern& n, + bifstream* pbifs = nullptr); +}; + +// PencilBeamAntenna format tags and XML I/O + +template <> +struct std::formatter { + format_tags tags; + + [[nodiscard]] constexpr auto& inner_fmt() { return *this; } + [[nodiscard]] constexpr auto& inner_fmt() const { return *this; } + + constexpr std::format_parse_context::iterator parse( + std::format_parse_context& ctx) { + return parse_format_tags(tags, ctx); + } + + template + FmtContext::iterator format(const sensor::PencilBeamAntenna& v, + FmtContext& ctx) const { + return tags.format(ctx, v.weight); + } +}; + +template <> +struct xml_io_stream_name { + static constexpr std::string_view name = "SensorPencilBeamAntenna"; +}; + +template <> +struct xml_io_stream { + static constexpr std::string_view type_name = + xml_io_stream_name_v; + + static void write(std::ostream& os, + const sensor::PencilBeamAntenna& n, + bofstream* pbofs = nullptr, + std::string_view name = ""sv); + + static void read(std::istream& is, + sensor::PencilBeamAntenna& n, + bifstream* pbifs = nullptr); +}; + +// GaussianAntenna format tags and XML I/O + +template <> +struct std::formatter + : format_tag_inherit {}; + +template <> +struct xml_io_stream_name { + static constexpr std::string_view name = "SensorGaussianAntenna"; +}; + +template <> +struct xml_io_stream + : xml_io_stream_inherit {}; + +// GaussianAiryAntenna format tags and XML I/O + +template <> +struct std::formatter { + format_tags tags; + + [[nodiscard]] constexpr auto& inner_fmt() { return *this; } + [[nodiscard]] constexpr auto& inner_fmt() const { return *this; } + + constexpr std::format_parse_context::iterator parse( + std::format_parse_context& ctx) { + return parse_format_tags(tags, ctx); + } + + template + FmtContext::iterator format(const sensor::GaussianAiryAntenna& v, + FmtContext& ctx) const { + auto sep = tags.sep(); + return tags.format(ctx, + v.zen_grid, + sep, + v.azi_grid, + sep, + v.aperture_diameter, + sep, + v.weight); + } +}; + +template <> +struct xml_io_stream_name { + static constexpr std::string_view name = "SensorGaussianAiryAntenna"; +}; + +template <> +struct xml_io_stream { + static constexpr std::string_view type_name = + xml_io_stream_name_v; + + static void write(std::ostream& os, + const sensor::GaussianAiryAntenna& n, + bofstream* pbofs = nullptr, + std::string_view name = ""sv); + + static void read(std::istream& is, + sensor::GaussianAiryAntenna& n, + bifstream* pbifs = nullptr); +}; diff --git a/src/core/sensor/frequency_bandpass_filters.cpp b/src/core/sensor/frequency_bandpass_filters.cpp new file mode 100644 index 0000000000..74cb591578 --- /dev/null +++ b/src/core/sensor/frequency_bandpass_filters.cpp @@ -0,0 +1,307 @@ +#include "frequency_bandpass_filters.h" + +#include + +#include +#include +#include +#include +#include +#include + +namespace sensor { +namespace { +constexpr Numeric response_eps = 64 * std::numeric_limits::epsilon(); + +Numeric scale(Numeric x) { return std::max(1.0, std::abs(x)); } + +bool is_close(Numeric a, Numeric b) { + return std::abs(a - b) <= response_eps * std::max(scale(a), scale(b)); +} + +Numeric sample_filter(const SortedGriddedField1& filter, Numeric f) { + const auto& grid = filter.grid<0>(); + + if (grid.empty() or f < grid.front() or f > grid.back()) return 0.0; + + auto it = stdr::lower_bound(grid, f); + + if (it == grid.begin()) return filter.data.front(); + if (it == grid.end()) return filter.data.back(); + + const Size upper = static_cast(std::distance(grid.begin(), it)); + if (is_close(*it, f)) return filter.data[upper]; + + const Size lower = upper - 1; + const Numeric x0 = grid[lower]; + const Numeric x1 = grid[upper]; + const Numeric y0 = filter.data[lower]; + const Numeric y1 = filter.data[upper]; + const Numeric t = (f - x0) / (x1 - x0); + + return y0 + t * (y1 - y0); +} + +void add_support_points(std::vector& points, + const AscendingGrid& grid, + Numeric low, + Numeric high) { + if (grid.empty() or low > high) return; + + points.push_back(low); + if (not is_close(low, high)) points.push_back(high); + + auto lower = stdr::lower_bound(grid, low); + auto upper = stdr::upper_bound(grid, high); + for (auto it = lower; it != upper; ++it) points.push_back(*it); +} + +void sort_unique(std::vector& points) { + stdr::sort(points); + points.erase(std::unique(points.begin(), points.end(), &is_close), + points.end()); +} + +SortedGriddedField1 make_filter( + const std::vector>& samples, + std::string_view name) { + std::vector grid(samples.size()); + Vector data(samples.size()); + + for (Size i = 0; i < samples.size(); i++) { + grid[i] = samples[i].first; + data[i] = samples[i].second; + } + + return {.data_name = String{name}, + .data = std::move(data), + .grid_names = std::array{"frequency"s}, + .grids = std::array{AscendingGrid{grid}}}; +} + +void deduplicate_zero_frequency_images( + std::vector>& samples) { + std::vector> unique_samples; + unique_samples.reserve(samples.size()); + + for (const auto& sample : samples) { + const bool duplicate = stdr::any_of(unique_samples, [&](const auto& kept) { + return is_close(kept.first, sample.first); + }); + + if (not duplicate) unique_samples.push_back(sample); + } + + samples = std::move(unique_samples); +} + +using InterpolationPoint = + std::variant, + lagrange_interp::lag_t<1, lagrange_interp::identity>>; + +struct FoldedLocalPoint { + InterpolationPoint interpolation; + std::vector> folded_samples; +}; + +InterpolationPoint make_interpolation_point(const AscendingGrid& grid, + Numeric f) { + return lagrange_interp::variant_lag(grid, f); +} + +Numeric interpolate_channel_weight(const Vector& weights, + const InterpolationPoint& point) { + return std::visit( + [&weights](const auto& lag) { + return lagrange_interp::interp(weights, lag); + }, + point); +} + +std::vector collect_local_points(const FrequencyRange& range, + const AscendingGrid& channel_grid) { + std::vector local_points; + + for (const auto& path : range.paths()) { + if (channel_grid.empty()) continue; + + const Numeric local_low = + std::max(path.local_range[0], channel_grid.front()); + const Numeric local_high = + std::min(path.local_range[1], channel_grid.back()); + + if (local_low > local_high and not is_close(local_low, local_high)) { + continue; + } + + add_support_points(local_points, channel_grid, local_low, local_high); + for (const auto& filter : path.filters) { + add_support_points(local_points, filter.grid<0>(), local_low, local_high); + } + } + + sort_unique(local_points); + return local_points; +} + +std::vector> fold_unit_samples( + const FrequencyRange& range, Numeric local_frequency) { + std::vector> folded_samples; + folded_samples.reserve(range.size()); + + for (const auto& path : range.paths()) { + const Numeric path_weight = path.local_weight(local_frequency); + if (path_weight == 0.0) continue; + + folded_samples.emplace_back(path.map_to_global(local_frequency), + path_weight); + } + + if (folded_samples.empty()) return folded_samples; + + const Numeric fold_count = static_cast(folded_samples.size()); + for (auto& sample : folded_samples) { + sample.second /= fold_count; + } + + if (is_close(local_frequency, 0.0)) { + deduplicate_zero_frequency_images(folded_samples); + } + + return folded_samples; +} + +void combine_samples(std::vector>& samples) { + std::sort(samples.begin(), samples.end(), [](const auto& a, const auto& b) { + return a.first < b.first; + }); + + std::vector> combined; + combined.reserve(samples.size()); + for (const auto& sample : samples) { + if (combined.empty() or not is_close(combined.back().first, sample.first)) { + combined.push_back(sample); + } else { + combined.back().second += sample.second; + } + } + + samples = std::move(combined); +} + +std::vector build_channels( + const FrequencyRange& range, const std::span& channels) { + std::vector out; + out.reserve(channels.size()); + + for (Size ichan = 0; ichan < channels.size(); ichan++) { + const auto& channel = channels[ichan]; + const auto& channel_grid = channel.freq_grid(); + std::vector> samples; + std::vector local_points; + + for (const auto& path : range.paths()) { + if (channel_grid.empty()) continue; + + const Numeric local_low = + std::max(path.local_range[0], channel_grid.front()); + const Numeric local_high = + std::min(path.local_range[1], channel_grid.back()); + + if (local_low > local_high and not is_close(local_low, local_high)) + continue; + + add_support_points(local_points, channel_grid, local_low, local_high); + for (const auto& filter : path.filters) { + add_support_points( + local_points, filter.grid<0>(), local_low, local_high); + } + } + + sort_unique(local_points); + + for (Numeric local_frequency : local_points) { + const Numeric channel_weight = + sample_filter(channel.channel, local_frequency); + if (channel_weight == 0.0) continue; + + auto folded_samples = fold_unit_samples(range, local_frequency); + if (folded_samples.empty()) continue; + + for (const auto& sample : folded_samples) { + samples.emplace_back(sample.first, sample.second * channel_weight); + } + } + + combine_samples(samples); + out.push_back( + Channel{.channel = make_filter(samples, + std::format("channel-response-{}", + ichan))}); + } + + return out; +} + +std::vector build_synced_channels( + const FrequencyRange& range, const Spectrometer& spectrometer) { + const auto& channels = spectrometer.channels; + if (channels.empty()) return {}; + + if (not spectrometer.is_synced()) { + return build_channels( + range, std::span{channels.data(), channels.size()}); + } + + const auto& channel_grid = channels.front().freq_grid(); + const auto local_points = collect_local_points(range, channel_grid); + + std::vector folded_points; + folded_points.reserve(local_points.size()); + for (Numeric local_frequency : local_points) { + auto folded_samples = fold_unit_samples(range, local_frequency); + if (folded_samples.empty()) continue; + + folded_points.push_back( + {make_interpolation_point(channel_grid, local_frequency), + std::move(folded_samples)}); + } + + std::vector out; + out.reserve(channels.size()); + + for (Size ichan = 0; ichan < channels.size(); ++ichan) { + const auto& weights = channels[ichan].weights(); + std::vector> samples; + + for (const auto& point : folded_points) { + const Numeric channel_weight = + interpolate_channel_weight(weights, point.interpolation); + if (channel_weight == 0.0) continue; + + for (const auto& sample : point.folded_samples) { + samples.emplace_back(sample.first, sample.second * channel_weight); + } + } + + combine_samples(samples); + out.push_back( + Channel{.channel = make_filter(samples, + std::format("channel-response-{}", + ichan))}); + } + + return out; +} +} // namespace + +std::vector make_bandpass_channels( + const FrequencyRange& range, const std::span& channels) { + return build_channels(range, channels); +} + +std::vector make_bandpass_channels( + const FrequencyRange& range, const Spectrometer& spectrometer) { + return build_synced_channels(range, spectrometer); +} +} // namespace sensor diff --git a/src/core/sensor/frequency_bandpass_filters.h b/src/core/sensor/frequency_bandpass_filters.h new file mode 100644 index 0000000000..1035aa60ae --- /dev/null +++ b/src/core/sensor/frequency_bandpass_filters.h @@ -0,0 +1,17 @@ +#pragma once + +#include + +#include + +#include "frequency_channel_selection.h" +#include "frequency_range_selection.h" + +namespace sensor { +[[nodiscard]] std::vector make_bandpass_channels( + const FrequencyRange& range, const std::span& channels); + +[[nodiscard]] std::vector make_bandpass_channels( + const FrequencyRange& range, const Spectrometer& spectrometer); +} // namespace sensor + diff --git a/src/core/sensor/frequency_channel_selection.cpp b/src/core/sensor/frequency_channel_selection.cpp new file mode 100644 index 0000000000..5e28582a2f --- /dev/null +++ b/src/core/sensor/frequency_channel_selection.cpp @@ -0,0 +1,165 @@ +#include "frequency_channel_selection.h" + +#include + +#include +#include +#include +#include +#include + +#include "compare.h" + +namespace sensor { +const AscendingGrid& Channel::freq_grid() const { return channel.grid<0>(); } + +const Vector& Channel::weights() const { return channel.data; } + +DiracChannel::DiracChannel(Numeric f) + : Channel{.channel = { + .data_name = "dirac"s, + .data = Vector(1, 1.0), + .grid_names = std::array{"frequency"s}, + .grids = std::array{AscendingGrid{f}}}} {} + +DiracChannel::DiracChannel() : DiracChannel(0) {} + +BoxChannel::BoxChannel(AscendingGrid f) + : Channel{ + .channel = { + .data_name = "box"s, + .data = Vector(f.size(), 1.0 / static_cast(f.size())), + .grid_names = std::array{"frequency"s}, + .grids = std::array{std::move(f)}}} {} + +BoxChannel::BoxChannel(Numeric lower, Numeric upper, Size N) + : BoxChannel{nlinspace(lower, upper, N)} {} + +BoxChannel::BoxChannel(Numeric hw, Size N) : BoxChannel(-hw, hw, N) {} + +namespace { +SortedGriddedField1 gauss(AscendingGrid&& f, Numeric f0, Numeric std) { + using gauss = boost::math::normal_distribution; + + Vector data{std::from_range, + f | stdv::transform([dist = gauss(f0, std)](Numeric fi) { + return pdf(dist, fi); + })}; + + return {.data_name = "gaussian"s, + .data = std::move(data), + .grid_names = std::array{"frequency"s}, + .grids = std::array{std::move(f)}}; +} +} // namespace + +GaussianChannel::GaussianChannel(AscendingGrid f, Numeric f0, Numeric std) + : Channel{.channel = gauss(std::move(f), f0, std)} {} + +GaussianChannel::GaussianChannel(Numeric f0, Numeric std, Size N, Size M) + : GaussianChannel( + nlinspace( + -static_cast(M) * std, static_cast(M) * std, N), + f0, + std) {} + +GaussianChannel::GaussianChannel(AscendingGrid f, Numeric std) + : GaussianChannel(std::move(f), 0, std) {} + +GaussianChannel::GaussianChannel(Numeric std, Size N, Size M) + : GaussianChannel( + nlinspace( + -static_cast(M) * std, static_cast(M) * std, N), + std) {} + +static_assert(FrequencyChannelSelection); +static_assert(FrequencyChannelSelection); +static_assert(FrequencyChannelSelection); +static_assert(FrequencyChannelSelection); + +void Spectrometer::sync_frequency_grids() { + std::vector f; + f.reserve(std::transform_reduce(channels.begin(), + channels.end(), + Size{}, + std::plus<>{}, + [](const Channel& c) { return c.size(); })); + + for (const auto& channel : channels) { + const auto& freqs = channel.freq_grid(); + f.insert(f.end(), freqs.begin(), freqs.end()); + } + + stdr::sort(f); + f.erase(std::unique(f.begin(), f.end()), f.end()); + + const auto fs = AscendingGrid(std::move(f)); + Vector ws(fs.size(), 0); + + for (Channel& channel : channels) { + const auto& ch_fs = channel.freq_grid(); + const auto& ch_ws = channel.weights(); + const Size n = ch_fs.size(); + + const auto fs_ptr_first = stdr::lower_bound(fs, ch_fs.front()); + const Size OFFSET = std::distance(fs.begin(), fs_ptr_first); + const auto ws_ptr_first = ws.begin() + OFFSET; + + auto fs_ptr = fs_ptr_first; + auto ws_ptr = ws_ptr_first; + for (Size i = 0; i < n; i++) { + while (*fs_ptr < ch_fs[i]) { + ++fs_ptr; + ++ws_ptr; + } + *ws_ptr = ch_ws[i]; + } + + channel.channel.grid<0>() = fs; + channel.channel.data = ws; + + fs_ptr = fs_ptr_first; + ws_ptr = ws_ptr_first; + for (Size i = 0; i < n; i++) { + while (*fs_ptr < ch_fs[i]) { + ++fs_ptr; + ++ws_ptr; + } + *ws_ptr = 0.0; + } + } +} + +bool Spectrometer::is_synced() const { + return channels.empty() or stdr::all_of(channels, + Cmp::eq(channels.front().freq_grid()), + &Channel::freq_grid); +} + +Spectrometer::Spectrometer(const Channel& base_channel, + const AscendingGrid& freq_offsets) + : channels(freq_offsets.size(), base_channel) { + for (Size i = 0; i < channels.size(); ++i) { + Vector f = channels[i].channel.grid<0>(); + f += freq_offsets[i]; + channels[i].channel.grid<0>() = f; + } + + sync_frequency_grids(); +} +} // namespace sensor + +void xml_io_stream::write( + std::ostream& os, + const sensor::Spectrometer& spec, + bofstream* pbofs, + std::string_view name) { + xml_io_stream>::write( + os, spec.channels, pbofs, name); +} + +void xml_io_stream::read(std::istream& is, + sensor::Spectrometer& spec, + bifstream* pbifs) { + xml_io_stream>::read(is, spec.channels, pbifs); +} diff --git a/src/core/sensor/frequency_channel_selection.h b/src/core/sensor/frequency_channel_selection.h new file mode 100644 index 0000000000..cc58dc716e --- /dev/null +++ b/src/core/sensor/frequency_channel_selection.h @@ -0,0 +1,146 @@ +#pragma once + +#include +#include +#include + +namespace sensor { +//! Free-form channel struct. Others inherit from this. +struct Channel; + +//! A channel that is even between a lower and upper frequency, with equal weights. +struct BoxChannel; + +//! A channel with a single frequency and all weight on that frequency. +struct DiracChannel; + +//! A channel with Gaussian weights centered at the center frequency with given standard deviation. +struct GaussianChannel; + +//! Concept for selecting frequencies for a channel +template +concept FrequencyChannelSelection = std::derived_from; + +struct Channel { + SortedGriddedField1 channel; + + [[nodiscard]] const AscendingGrid& freq_grid() const; + [[nodiscard]] const Vector& weights() const; + [[nodiscard]] Size size() const { return channel.data.size(); } +}; + +struct BoxChannel final : Channel { + BoxChannel(Numeric lower, Numeric upper, Size N); // [lower, upper] + BoxChannel(Numeric hw, Size N); // [-hw, hw] + BoxChannel(AscendingGrid f); // f + BoxChannel() = default; // [0, 0] +}; + +struct DiracChannel final : Channel { + DiracChannel(Numeric f); // f + DiracChannel(); // f = 0 +}; + +struct GaussianChannel final : Channel { + GaussianChannel(AscendingGrid f, Numeric f0, Numeric std); // std around f0 + GaussianChannel(Numeric f0, Numeric std, Size N, Size M); // f0 +- M*std + GaussianChannel(AscendingGrid f, Numeric std); // f, f0 = 0 + GaussianChannel(Numeric std, Size N, Size M); // +-M*std, f0 = 0 + GaussianChannel() = default; // f0 = 0, std = 0 + + [[nodiscard]] AscendingGrid center_grid() const; +}; + +struct Spectrometer { + std::vector channels; + + private: + void sync_frequency_grids(); + + public: + [[nodiscard]] bool is_synced() const; + + /** Make a simple spectrometer + + The spectrometer will have the same channel response for each position + but will be shifted by the frequencies in freq_offsets. + + So if base_channel is f: [-1, 1], w: [0.5, 0.5], and freq_offset is [2, 3, 4], + the resulting channels will be [[1, 2], [2, 3], [3, 4]] with the weights [0.5, 0.5] for each channel. + + But, since this is supposed to simulate a spectrometer, the true frequency grid of each channel + will be [1, 2, 3, 4] for all channels but the weights will be [0.5, 0.5, 0, 0], [0, 0.5, 0.5, 0], and + [0, 0, 0.5, 0.5] for each channel respectively. + + @param[in] base_channel The base channel + @param[in] freq_offsets The frequency offsets for each channel +*/ + Spectrometer(const Channel& base_channel, const AscendingGrid& freq_offsets); + + Spectrometer(std::vector channels) : channels(std::move(channels)) {} + + operator const std::vector&() const { return channels; } +}; +} // namespace sensor + +// Channel format tags and XML I/O + +template <> +struct format_tag_aggregate { + constexpr static bool value = true; +}; + +template <> +struct xml_io_stream_name { + static constexpr std::string_view name = "SensorChannel"; +}; + +template <> +struct xml_io_stream_aggregate { + static constexpr bool value = true; +}; + +// BoxChannel format tags and XML I/O + +template <> +struct std::formatter + : format_tag_inherit {}; + +template <> +struct xml_io_stream + : xml_io_stream_inherit {}; + +// DiracChannel format tags and XML I/O + +template <> +struct std::formatter + : format_tag_inherit {}; + +template <> +struct xml_io_stream + : xml_io_stream_inherit {}; + +// GaussianChannel format tags and XML I/O + +template <> +struct std::formatter + : format_tag_inherit {}; + +template <> +struct xml_io_stream + : xml_io_stream_inherit {}; + +// Spectrometer format tags and XML I/O + +template <> +struct std::formatter + : std::formatter> {}; + +template <> +struct xml_io_stream { + static void write(std::ostream&, + const sensor::Spectrometer&, + bofstream* = nullptr, + std::string_view = ""sv); + static void read(std::istream&, sensor::Spectrometer&, bifstream* = nullptr); +}; diff --git a/src/core/sensor/frequency_range_selection.cpp b/src/core/sensor/frequency_range_selection.cpp new file mode 100644 index 0000000000..bf3a575acc --- /dev/null +++ b/src/core/sensor/frequency_range_selection.cpp @@ -0,0 +1,343 @@ +#include "frequency_range_selection.h" + +#include +#include +#include +#include + +#include "matpack_mdspan_cdata_t.h" + +namespace sensor { +Size FrequencyRange::size() const { return response_paths.size(); } + +const FrequencyResponsePath& FrequencyRange::path(Size index) const { + return response_paths.at(index); +} + +const std::vector& FrequencyRange::paths() const { + return response_paths; +} + +void FrequencyRange::sync_ranges() { + global_ranges.clear(); + local_ranges.clear(); + + global_ranges.reserve(response_paths.size()); + local_ranges.reserve(response_paths.size()); + + for (const auto& path : response_paths) { + global_ranges.push_back(path.global_range()); + local_ranges.push_back(path.local_range); + } +} + +namespace { +constexpr Numeric response_eps = 64 * std::numeric_limits::epsilon(); + +Numeric scale(Numeric x) { return std::max(1.0, std::abs(x)); } + +bool is_close(Numeric a, Numeric b) { + if (std::isinf(a) or std::isinf(b)) return a == b; + return std::abs(a - b) <= response_eps * std::max(scale(a), scale(b)); +} + +bool is_in_closed_interval(Numeric x, Numeric low, Numeric high) { + return x >= low - response_eps * scale(low) and + x <= high + response_eps * scale(high); +} + +void assert_frange(const Vector2& sideband) { + if (sideband[0] > sideband[1] or is_close(sideband[0], sideband[1]) or + sideband[0] < 0) { + throw std::invalid_argument(std::format( + "Frequency range must be unique, sorted, and non-negative. Got: [{}, {}].", + sideband[0], + sideband[1])); + } +} + +void assert_bandpass_ranges(const std::span& bandpasses) { + if (bandpasses.empty()) { + throw std::invalid_argument("At least one sideband must be provided."); + } + + for (auto sideband : bandpasses) assert_frange(sideband); +} + +void assert_lo_nonnegative(const Numeric& LO) { + if (LO < 0) { + throw std::invalid_argument( + std::format("LO frequency must be non-negative. Got: {}.", LO)); + } +} + +void assert_clocks(const std::span& clocks) { + if (clocks.empty()) { + throw std::invalid_argument("At least one LO frequency must be provided."); + } + + for (auto clock : clocks) assert_lo_nonnegative(clock); +} + +void assert_weighted_bandpass(const SortedGriddedField1& bandpass_filter) { + const auto& grid = bandpass_filter.grid<0>(); + + if (grid.empty()) { + throw std::invalid_argument("Bandpass filter grid must not be empty."); + } + + if (grid.size() != bandpass_filter.data.size()) { + throw std::invalid_argument(std::format( + "Bandpass filter grid and data must have the same size. Got {} grid points and {} weights.", + grid.size(), + bandpass_filter.data.size())); + } + + for (Size i = 1; i < grid.size(); i++) { + if (grid[i - 1] > grid[i] or is_close(grid[i - 1], grid[i])) { + throw std::invalid_argument( + "Bandpass filter grid must be unique and sorted."); + } + } +} + +Numeric sample_filter(const SortedGriddedField1& filter, Numeric f) { + const auto& grid = filter.grid<0>(); + + if (grid.empty() or f < grid.front() or f > grid.back()) return 0.0; + + auto it = std::lower_bound(grid.begin(), grid.end(), f); + + if (it == grid.begin()) return filter.data.front(); + if (it == grid.end()) return filter.data.back(); + + const Size upper = static_cast(std::distance(grid.begin(), it)); + if (is_close(*it, f)) return filter.data[upper]; + + const Size lower = upper - 1; + const Numeric x0 = grid[lower]; + const Numeric x1 = grid[upper]; + const Numeric y0 = filter.data[lower]; + const Numeric y1 = filter.data[upper]; + + const Numeric t = (f - x0) / (x1 - x0); + return y0 + t * (y1 - y0); +} + +std::vector transformed_filters( + const std::vector& filters, Numeric LO, bool upper) { + std::vector out; + out.reserve(filters.size()); + + for (const auto& filter : filters) { + const auto& grid = filter.grid<0>(); + std::vector mapped(grid.size()); + Vector weights(filter.data.size()); + + if (upper) { + for (Size i = 0; i < grid.size(); i++) { + mapped[i] = grid[i] - LO; + weights[i] = filter.data[i]; + } + } else { + for (Size i = 0; i < grid.size(); i++) { + const Size j = grid.size() - 1 - i; + mapped[i] = LO - grid[j]; + weights[i] = filter.data[j]; + } + } + + out.push_back( + {.data_name = filter.data_name, + .data = std::move(weights), + .grid_names = filter.grid_names, + .grids = std::array{AscendingGrid{mapped}}}); + } + + return out; +} + +void apply_interval_clip(std::vector& paths, + Numeric low, + Numeric high) { + std::erase_if(paths, [low, high](FrequencyResponsePath& path) { + path.local_range[0] = std::max(path.local_range[0], low); + path.local_range[1] = std::min(path.local_range[1], high); + return path.local_range[0] >= path.local_range[1] or + is_close(path.local_range[0], path.local_range[1]); + }); +} +} // namespace + +Vector2 FrequencyResponsePath::global_range() const { + return {map_to_global(local_range[0]), map_to_global(local_range[1])}; +} + +Numeric FrequencyResponsePath::map_to_global(Numeric local_frequency) const { + return intercept + slope * local_frequency; +} + +std::optional FrequencyResponsePath::map_to_local( + Numeric global_frequency) const { + const Numeric local_frequency = (global_frequency - intercept) / slope; + + if (not is_in_closed_interval( + local_frequency, local_range[0], local_range[1])) { + return std::nullopt; + } + + return std::clamp(local_frequency, local_range[0], local_range[1]); +} + +Numeric FrequencyResponsePath::local_weight(Numeric local_frequency) const { + if (not is_in_closed_interval( + local_frequency, local_range[0], local_range[1])) { + return 0.0; + } + + Numeric weight = 1.0; + + for (const auto& filter : filters) { + weight *= sample_filter(filter, local_frequency); + if (weight == 0.0) break; + } + + return weight; +} + +Numeric FrequencyResponsePath::global_weight(Numeric global_frequency) const { + const auto local_frequency = map_to_local(global_frequency); + return local_frequency.has_value() ? local_weight(*local_frequency) : 0.0; +} + +HeterodyneFrequencyRange::HeterodyneFrequencyRange( + const std::span& clocks, + const std::span& bandpasses) + : FrequencyRange() { + const Size N = clocks.size(); + + if (N != bandpasses.size()) { + throw std::invalid_argument(std::format( + "Number of clock frequencies and bandpasses must match. Got: {} clock " + "frequencies and {} bandpasses.", + clocks.size(), + bandpasses.size())); + } + + assert_bandpass_ranges(bandpasses); + assert_clocks(clocks); + + for (Size i = 0; i < N; i++) { + apply_bandpass(bandpasses[i]); + apply_mixer(clocks[i]); + } +} + +HeterodyneFrequencyRange::HeterodyneFrequencyRange(Numeric clock_frequency, + const Vector2& sideband) + : HeterodyneFrequencyRange(std::span{&clock_frequency, 1}, + std::span{&sideband, 1}) {} + +void HeterodyneFrequencyRange::apply_lowpass(Numeric upper_frequency) { + if (upper_frequency < 0) { + throw std::invalid_argument(std::format( + "Lowpass cutoff must be non-negative. Got: {}.", upper_frequency)); + } + + apply_interval_clip(response_paths, 0.0, upper_frequency); + sync_ranges(); +} + +void HeterodyneFrequencyRange::apply_highpass(Numeric lower_frequency) { + if (lower_frequency < 0) { + throw std::invalid_argument(std::format( + "Highpass cutoff must be non-negative. Got: {}.", lower_frequency)); + } + + apply_interval_clip(response_paths, lower_frequency, inf); + sync_ranges(); +} + +void HeterodyneFrequencyRange::apply_bandpass(const Vector2& bandpass_range) { + assert_frange(bandpass_range); + apply_interval_clip(response_paths, bandpass_range[0], bandpass_range[1]); + sync_ranges(); +} + +void HeterodyneFrequencyRange::apply_bandpass( + const SortedGriddedField1& bandpass_filter) { + assert_weighted_bandpass(bandpass_filter); + + apply_interval_clip(response_paths, + bandpass_filter.grid<0>().front(), + bandpass_filter.grid<0>().back()); + + for (auto& path : response_paths) path.filters.push_back(bandpass_filter); + + sync_ranges(); +} + +void HeterodyneFrequencyRange::apply_mixer(Numeric clock_frequency) { + assert_lo_nonnegative(clock_frequency); + + std::vector mixed_paths; + mixed_paths.reserve(2 * response_paths.size()); + + for (const auto& path : response_paths) { + const Numeric upper_low = std::max(path.local_range[0], clock_frequency); + const Numeric upper_high = path.local_range[1]; + + if (upper_low < upper_high and not is_close(upper_low, upper_high)) { + auto upper = path; + upper.intercept += upper.slope * clock_frequency; + upper.local_range = {upper_low - clock_frequency, + upper_high - clock_frequency}; + upper.filters = transformed_filters(path.filters, clock_frequency, true); + mixed_paths.push_back(std::move(upper)); + } + + const Numeric lower_low = path.local_range[0]; + const Numeric lower_high = std::min(path.local_range[1], clock_frequency); + + if (lower_low < lower_high and not is_close(lower_low, lower_high)) { + auto lower = path; + lower.intercept += lower.slope * clock_frequency; + lower.slope *= -1; + lower.local_range = {clock_frequency - lower_high, + clock_frequency - lower_low}; + lower.filters = transformed_filters(path.filters, clock_frequency, false); + mixed_paths.push_back(std::move(lower)); + } + } + + response_paths = std::move(mixed_paths); + sync_ranges(); +} + +Vector HeterodyneFrequencyRange::local_response( + ConstVectorView local_frequency_grid, Size path_index) const { + const auto& response_path = path(path_index); + Vector response(local_frequency_grid.size()); + + for (Size i = 0; i < response.size(); i++) { + response[i] = response_path.local_weight(local_frequency_grid[i]); + } + + return response; +} + +Vector HeterodyneFrequencyRange::global_response( + ConstVectorView global_frequency_grid, Size path_index) const { + const auto& response_path = path(path_index); + Vector response(global_frequency_grid.size()); + + for (Size i = 0; i < response.size(); i++) { + response[i] = response_path.global_weight(global_frequency_grid[i]); + } + + return response; +} + +static_assert(FrequencyRangeSelection); +static_assert(FrequencyRangeSelection); +} // namespace sensor diff --git a/src/core/sensor/frequency_range_selection.h b/src/core/sensor/frequency_range_selection.h new file mode 100644 index 0000000000..7e284d0fce --- /dev/null +++ b/src/core/sensor/frequency_range_selection.h @@ -0,0 +1,117 @@ +#pragma once + +#include + +#include +#include + +namespace sensor { +//! Free-form frequency range struct. Others inherit from this. +struct FrequencyRange; + +//! A frequency range for heterodyne channels. True frequencies are computed. +struct HeterodyneFrequencyRange; + +//! A single affine path through a staged heterodyne chain. +struct FrequencyResponsePath; + +//! Concept for selecting frequency ranges for a set of channels +template +concept FrequencyRangeSelection = std::derived_from; + +struct FrequencyResponsePath { + static constexpr Numeric inf = std::numeric_limits::infinity(); + + Numeric intercept{0.0}; + Numeric slope{1.0}; + Vector2 local_range{0.0, inf}; + std::vector filters{}; + + [[nodiscard]] Vector2 global_range() const; + [[nodiscard]] Numeric map_to_global(Numeric local_frequency) const; + [[nodiscard]] std::optional map_to_local( + Numeric global_frequency) const; + [[nodiscard]] Numeric local_weight(Numeric local_frequency) const; + [[nodiscard]] Numeric global_weight(Numeric global_frequency) const; +}; + +struct FrequencyRange { + static constexpr Numeric inf = FrequencyResponsePath::inf; + + std::vector response_paths{{}}; + + std::vector global_ranges{{0, inf}}; + std::vector local_ranges{{0, inf}}; + + [[nodiscard]] Size size() const; + [[nodiscard]] const FrequencyResponsePath& path(Size index) const; + [[nodiscard]] const std::vector& paths() const; + + void sync_ranges(); +}; + +struct HeterodyneFrequencyRange final : FrequencyRange { + HeterodyneFrequencyRange() = default; + HeterodyneFrequencyRange(Numeric clock_frequency, + const Vector2& bandpass_range); + HeterodyneFrequencyRange(const std::span& clock_frequencies, + const std::span& bandpass_ranges); + + void apply_lowpass(Numeric upper_frequency); + void apply_highpass(Numeric lower_frequency); + void apply_bandpass(const Vector2& bandpass_range); + void apply_bandpass(const SortedGriddedField1& bandpass_filter); + void apply_mixer(Numeric clock_frequency); + + [[nodiscard]] Vector local_response(ConstVectorView local_frequency_grid, + Size path_index = 0) const; + [[nodiscard]] Vector global_response(ConstVectorView global_frequency_grid, + Size path_index = 0) const; +}; +} // namespace sensor + +// FrequencyResponsePath format tags and XML I/O + +template <> +struct format_tag_aggregate { + constexpr static bool value = true; +}; + +template <> +struct xml_io_stream_name { + static constexpr std::string_view name = "SensorFrequencyResponsePath"; +}; + +template <> +struct xml_io_stream_aggregate { + static constexpr bool value = true; +}; + +// FrequencyRange format tags and XML I/O + +template <> +struct format_tag_aggregate { + constexpr static bool value = true; +}; + +template <> +struct xml_io_stream_name { + static constexpr std::string_view name = "SensorFrequencyRange"; +}; + +template <> +struct xml_io_stream_aggregate { + static constexpr bool value = true; +}; + +// HeterodyneFrequencyRange format tags and XML I/O + +template <> +struct std::formatter + : format_tag_inherit {}; + +template <> +struct xml_io_stream + : xml_io_stream_inherit {}; diff --git a/src/core/sensor/obsel.cpp b/src/core/sensor/obsel.cpp index 5836656f4a..1c1cedc554 100644 --- a/src/core/sensor/obsel.cpp +++ b/src/core/sensor/obsel.cpp @@ -1,11 +1,13 @@ #include "obsel.h" +#include #include #include #include #include #include +#include #include bool SensorKey::operator==(const SensorKey& other) const { @@ -516,7 +518,7 @@ void make_exclusive(std::span obsels) { for (const auto& w : elem.weight_matrix()) { const Index iposlos = stdr::distance(stdr::begin(gs), stdr::find(gs, old_gs[w.irow])); - const Index ifreq = stdr::distance(stdr::begin(fs), + const Index ifreq = stdr::distance(stdr::begin(fs), stdr::lower_bound(fs, old_fs[w.icol])); ws[iposlos, ifreq] = w.data; } @@ -525,6 +527,41 @@ void make_exclusive(std::span obsels) { } } +void collect_frequency_grids(std::span obsels) { + if (obsels.empty()) return; + + std::set> freq_set; + for (const auto& obsel : obsels) freq_set.insert(obsel.f_grid_ptr()); + + if (freq_set.size() == 1) return; + + Size total_size = 0; + for (const auto& freqs : freq_set) total_size += freqs->size(); + + std::vector fs; + fs.reserve(total_size); + for (const auto& freqs : freq_set) fs.append_range(*freqs); + + stdr::sort(fs); + fs.erase(std::unique(fs.begin(), fs.end()), fs.end()); + + const auto freq_ptr = std::make_shared(std::move(fs)); + +#pragma omp parallel for if (not arts_omp_in_parallel()) + for (auto& obsel : obsels) { + sensor::SparseStokvecMatrix w(obsel.poslos_grid().size(), freq_ptr->size()); + + for (const auto& we : obsel.weight_matrix()) { + const Index ifreq = + stdr::distance(freq_ptr->begin(), + stdr::lower_bound(*freq_ptr, obsel.f_grid()[we.icol])); + w[we.irow, ifreq] = we.data; + } + + obsel = {freq_ptr, obsel.poslos_grid_ptr(), std::move(w)}; + } +} + namespace { void set_frq(const SensorObsel& v, ArrayOfSensorObsel& sensor, @@ -537,7 +574,7 @@ void set_frq(const SensorObsel& v, const auto xs = std::make_shared( x.begin(), x.end(), [](auto& x) { return x; }); - // Must copy, as we may change the shared_ptr later + // Copy the shared_ptr value before updating matching obsels. const auto fs = v.f_grid_ptr(); for (auto& elem : sensor) { @@ -573,12 +610,12 @@ void set_poslos(const SensorObsel& v, const auto xs = std::make_shared(std::move(xsv)); - // Must copy, as we may change the shared_ptr later + // Copy the shared_ptr value before updating matching obsels. const auto ps = v.poslos_grid_ptr(); for (auto& elem : sensor) { if (elem.poslos_grid_ptr() == ps) { - elem.set_poslos_grid_ptr(ps); // may change here + elem.set_poslos_grid_ptr(xs); // may change here } } } diff --git a/src/core/sensor/obsel.h b/src/core/sensor/obsel.h index f11a4c0e43..8674588ad9 100644 --- a/src/core/sensor/obsel.h +++ b/src/core/sensor/obsel.h @@ -7,8 +7,6 @@ #include #include -#include -#include #include "matpack_mdspan_helpers_grid_t.h" @@ -265,6 +263,8 @@ void make_exhaustive(std::span obsels); */ void make_exclusive(std::span obsels); +void collect_frequency_grids(std::span obsels); + std::vector unique_frequency_grids( const std::span& obsels); diff --git a/src/core/sensor/sensor_builder.cpp b/src/core/sensor/sensor_builder.cpp new file mode 100644 index 0000000000..f1404b1a3a --- /dev/null +++ b/src/core/sensor/sensor_builder.cpp @@ -0,0 +1,120 @@ +#include "sensor_builder.h" + +#include +#include + +#include +#include + +namespace sensor { +namespace { +std::shared_ptr clone_antenna( + const std::shared_ptr& antenna) { + return antenna ? antenna->clone() : nullptr; +} + +SensorMetaInfo make_meta_info(Size nchannels, Size geometry_index) { + SortedGriddedField1 gf; + gf.data_name = std::format("sensor-builder-{}", geometry_index); + gf.gridname<0>() = "channel"; + Vector channel_axis(nchannels); + for (Size i = 0; i < nchannels; ++i) { + channel_axis[i] = static_cast(i); + } + gf.grid<0>() = AscendingGrid{std::move(channel_axis)}; + gf.data.resize(nchannels); + gf.data = 0.0; + + return SensorMetaInfo{std::move(gf)}; +} +} // namespace + +SensorBuilder::SensorBuilder() : antenna(PencilBeamAntenna{}.clone()) {} + +SensorBuilder::SensorBuilder(std::vector channels, + std::shared_ptr antenna) + : channels(std::move(channels)), antenna(std::move(antenna)) {} + +SensorBuilder::SensorBuilder(const Spectrometer& spectrometer, + std::shared_ptr antenna) + : SensorBuilder(std::vector{spectrometer.channels}, + std::move(antenna)) { + preserve_common_frequency_grid = spectrometer.is_synced(); +} + +SensorBuilder::SensorBuilder(const Spectrometer& spectrometer, + const FrequencyRange& backend, + std::shared_ptr antenna) + : SensorBuilder(make_bandpass_channels(backend, spectrometer), + std::move(antenna)) { + preserve_common_frequency_grid = spectrometer.is_synced(); +} + +SensorBuilder& SensorBuilder::operator=(const SensorBuilder& other) { + if (this != &other) { + channels = other.channels; + antenna = clone_antenna(other.antenna); + preserve_common_frequency_grid = other.preserve_common_frequency_grid; + } + + return *this; +} + +std::pair SensorBuilder::operator()( + std::span pos, std::span los) const { + ARTS_USER_ERROR_IF(channels.empty(), + "SensorBuilder requires at least one channel") + ARTS_USER_ERROR_IF(not antenna, "SensorBuilder requires an antenna pattern") + ARTS_USER_ERROR_IF(pos.empty(), + "SensorBuilder requires at least one sensor position") + ARTS_USER_ERROR_IF(los.empty(), + "SensorBuilder requires at least one bore LOS") + ARTS_USER_ERROR_IF( + pos.size() != los.size(), + "SensorBuilder requires matching position and LOS counts. Got {} positions and {} LOS values.", + pos.size(), + los.size()) + + std::vector> freq_grids; + freq_grids.reserve(channels.size()); + for (const auto& channel : channels) { + freq_grids.push_back( + std::make_shared(channel.freq_grid())); + } + + ArrayOfSensorObsel out; + out.reserve(pos.size() * channels.size()); + + ArrayOfSensorMetaInfo meta; + meta.reserve(pos.size()); + + const auto append_geometry = [&](const Vector3& sensor_pos, + const Vector2& bore_los, + Size geometry_index) { + std::shared_ptr poslos_grid; + + for (Size ichan = 0; ichan < channels.size(); ++ichan) { + auto obsel = antenna->operator()(channels[ichan], sensor_pos, bore_los); + obsel.set_f_grid_ptr(freq_grids[ichan]); + + if (not poslos_grid) { + poslos_grid = obsel.poslos_grid_ptr(); + } else { + obsel.set_poslos_grid_ptr(poslos_grid); + } + + out.emplace_back(std::move(obsel)); + } + + meta.push_back(make_meta_info(channels.size(), geometry_index)); + }; + + for (Size i = 0; i < pos.size(); ++i) append_geometry(pos[i], los[i], i); + + if (preserve_common_frequency_grid) collect_frequency_grids(out); + + return {std::move(out), std::move(meta)}; +} + +static_assert(SensorBuilderSelection); +} // namespace sensor diff --git a/src/core/sensor/sensor_builder.h b/src/core/sensor/sensor_builder.h new file mode 100644 index 0000000000..b376934c3b --- /dev/null +++ b/src/core/sensor/sensor_builder.h @@ -0,0 +1,44 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "antenna_pattern.h" +#include "frequency_channel_selection.h" +#include "obsel.h" +#include "sensor_meta_info.h" + +namespace sensor { +//! Combines channels with an antenna pattern and bore geometries into obsels. +struct SensorBuilder; +struct FrequencyRange; + +//! Concept for selecting sensor builders. +template +concept SensorBuilderSelection = std::derived_from; + +struct SensorBuilder { + std::vector channels; + std::shared_ptr antenna; + bool preserve_common_frequency_grid{false}; + + SensorBuilder(); + SensorBuilder(std::vector channels, + std::shared_ptr antenna); + SensorBuilder(const Spectrometer& spectrometer, + std::shared_ptr antenna); + SensorBuilder(const Spectrometer& spectrometer, + const FrequencyRange& backend, + std::shared_ptr antenna); + SensorBuilder(const SensorBuilder& other) = default; + SensorBuilder(SensorBuilder&&) noexcept = default; + SensorBuilder& operator=(const SensorBuilder& other); + SensorBuilder& operator=(SensorBuilder&&) noexcept = default; + + [[nodiscard]] std::pair operator()( + std::span pos, std::span los) const; +}; +} // namespace sensor diff --git a/src/core/util/aggregate_helper.h b/src/core/util/aggregate_helper.h new file mode 100644 index 0000000000..801201a82e --- /dev/null +++ b/src/core/util/aggregate_helper.h @@ -0,0 +1,448 @@ +#pragma once + +#include +#include + +namespace { +struct aggregate_init { + template + constexpr operator T() const noexcept; +}; + +template +concept aggregate_0 = std::is_aggregate_v and requires { T{}; }; + +template +concept aggregate_1 = aggregate_0 and requires { T{aggregate_init{}}; }; + +template +concept aggregate_2 = + aggregate_1 and requires { T{aggregate_init{}, aggregate_init{}}; }; + +template +concept aggregate_3 = aggregate_2 and requires { + T{aggregate_init{}, aggregate_init{}, aggregate_init{}}; +}; + +template +concept aggregate_4 = aggregate_3 and requires { + T{aggregate_init{}, aggregate_init{}, aggregate_init{}, aggregate_init{}}; +}; + +template +concept aggregate_5 = aggregate_4 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_6 = aggregate_5 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_7 = aggregate_6 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_8 = aggregate_7 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_9 = aggregate_8 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_10 = aggregate_9 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_11 = aggregate_10 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_12 = aggregate_11 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_13 = aggregate_12 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_14 = aggregate_13 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_15 = aggregate_14 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_16 = aggregate_15 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_17 = aggregate_16 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_18 = aggregate_17 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_19 = aggregate_18 and requires { + T{aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}, + aggregate_init{}}; +}; + +template +concept aggregate_20 = aggregate_19 and requires { + T{aggregate_init{}, aggregate_init{}, aggregate_init{}, aggregate_init{}, + aggregate_init{}, aggregate_init{}, aggregate_init{}, aggregate_init{}, + aggregate_init{}, aggregate_init{}, aggregate_init{}, aggregate_init{}, + aggregate_init{}, aggregate_init{}, aggregate_init{}, aggregate_init{}, + aggregate_init{}, aggregate_init{}, aggregate_init{}, aggregate_init{}}; +}; + +template +concept aggregate_0_exact = aggregate_0 and not aggregate_1; + +template +concept aggregate_1_exact = aggregate_1 and not aggregate_2; + +template +concept aggregate_2_exact = aggregate_2 and not aggregate_3; + +template +concept aggregate_3_exact = aggregate_3 and not aggregate_4; + +template +concept aggregate_4_exact = aggregate_4 and not aggregate_5; + +template +concept aggregate_5_exact = aggregate_5 and not aggregate_6; + +template +concept aggregate_6_exact = aggregate_6 and not aggregate_7; + +template +concept aggregate_7_exact = aggregate_7 and not aggregate_8; + +template +concept aggregate_8_exact = aggregate_8 and not aggregate_9; + +template +concept aggregate_9_exact = aggregate_9 and not aggregate_10; + +template +concept aggregate_10_exact = aggregate_10 and not aggregate_11; + +template +concept aggregate_11_exact = aggregate_11 and not aggregate_12; + +template +concept aggregate_12_exact = aggregate_12 and not aggregate_13; + +template +concept aggregate_13_exact = aggregate_13 and not aggregate_14; + +template +concept aggregate_14_exact = aggregate_14 and not aggregate_15; + +template +concept aggregate_15_exact = aggregate_15 and not aggregate_16; + +template +concept aggregate_16_exact = aggregate_16 and not aggregate_17; + +template +concept aggregate_17_exact = aggregate_17 and not aggregate_18; + +template +concept aggregate_18_exact = aggregate_18 and not aggregate_19; + +template +concept aggregate_19_exact = aggregate_19 and not aggregate_20; + +template +concept aggregate_20_exact = aggregate_20; +} // namespace + +auto as_tuple(aggregate_0_exact auto&) { return std::tuple<>(); } + +auto as_tuple(aggregate_1_exact auto& x) { + auto&& [a] = x; + return std::tie(a); +} + +auto as_tuple(aggregate_2_exact auto& x) { + auto&& [a, b] = x; + return std::tie(a, b); +} + +auto as_tuple(aggregate_3_exact auto& x) { + auto&& [a, b, c] = x; + return std::tie(a, b, c); +} + +auto as_tuple(aggregate_4_exact auto& x) { + auto&& [a, b, c, d] = x; + return std::tie(a, b, c, d); +} + +auto as_tuple(aggregate_5_exact auto& x) { + auto&& [a, b, c, d, e] = x; + return std::tie(a, b, c, d, e); +} + +auto as_tuple(aggregate_6_exact auto& x) { + auto&& [a, b, c, d, e, f] = x; + return std::tie(a, b, c, d, e, f); +} + +auto as_tuple(aggregate_7_exact auto& x) { + auto&& [a, b, c, d, e, f, g] = x; + return std::tie(a, b, c, d, e, f, g); +} + +auto as_tuple(aggregate_8_exact auto& x) { + auto&& [a, b, c, d, e, f, g, h] = x; + return std::tie(a, b, c, d, e, f, g, h); +} + +auto as_tuple(aggregate_9_exact auto& x) { + auto&& [a, b, c, d, e, f, g, h, i] = x; + return std::tie(a, b, c, d, e, f, g, h, i); +} + +auto as_tuple(aggregate_10_exact auto& x) { + auto&& [a, b, c, d, e, f, g, h, i, j] = x; + return std::tie(a, b, c, d, e, f, g, h, i, j); +} + +auto as_tuple(aggregate_11_exact auto& x) { + auto&& [a, b, c, d, e, f, g, h, i, j, k] = x; + return std::tie(a, b, c, d, e, f, g, h, i, j, k); +} + +auto as_tuple(aggregate_12_exact auto& x) { + auto&& [a, b, c, d, e, f, g, h, i, j, k, l] = x; + return std::tie(a, b, c, d, e, f, g, h, i, j, k, l); +} + +auto as_tuple(aggregate_13_exact auto& x) { + auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m] = x; + return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m); +} + +auto as_tuple(aggregate_14_exact auto& x) { + auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m, n] = x; + return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m, n); +} + +auto as_tuple(aggregate_15_exact auto& x) { + auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o] = x; + return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o); +} + +auto as_tuple(aggregate_16_exact auto& x) { + auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p] = x; + return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p); +} + +auto as_tuple(aggregate_17_exact auto& x) { + auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q] = x; + return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q); +} + +auto as_tuple(aggregate_18_exact auto& x) { + auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r] = x; + return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r); +} + +auto as_tuple(aggregate_19_exact auto& x) { + auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s] = x; + return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s); +} + +auto as_tuple(aggregate_20_exact auto& x) { + auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t] = x; + return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t); +} + +template +concept arts_aggregate = requires(T a) { as_tuple(a); }; diff --git a/src/core/util/arts_conversions.h b/src/core/util/arts_conversions.h index d174173cbe..70cd8f783b 100644 --- a/src/core/util/arts_conversions.h +++ b/src/core/util/arts_conversions.h @@ -166,6 +166,12 @@ constexpr auto angstrom2meter(auto x) noexcept { return x * 1e-10; } /** Conversion from meter to Å **/ constexpr auto meter2angstrom(auto x) noexcept { return x * 1e10; } +/** Conversion from HWHM to STD */ +constexpr auto hwhm2std(auto x) noexcept { return x / (Constant::sqrt_ln_2 * Constant::sqrt_2); } + +/** Conversion from FWHM to STD */ +constexpr auto fwhm2std(auto x) noexcept { return 0.5 * hwhm2std(x); } + //! Converts the number to a metric prefix (kilo:=k, Mega:=M, nano:=n, etc) // And empty char (' ') is used when no conversion happens. The function returns a pair // containing the prefix character and the scaled value. diff --git a/src/core/util/format_tags.h b/src/core/util/format_tags.h index be49ba8143..1d92d39eb1 100644 --- a/src/core/util/format_tags.h +++ b/src/core/util/format_tags.h @@ -17,6 +17,8 @@ #include #include +#include "aggregate_helper.h" + using namespace std::literals; template @@ -604,3 +606,55 @@ struct std::formatter> { return tags.format(ctx, ""); } }; + +template +struct format_tag_aggregate { + constexpr static bool value = false; +}; + +template +constexpr bool format_tag_aggregate_v = format_tag_aggregate::value; + +template +concept format_tag_aggratable = format_tag_aggregate_v; + +template +struct std::formatter { + static_assert(arts_aggregate, + "format_tag_aggratable types must be arts_aggregate"); + + format_tags tags; + + [[nodiscard]] constexpr auto& inner_fmt() { return *this; } + + [[nodiscard]] constexpr const auto& inner_fmt() const { return *this; } + + constexpr std::format_parse_context::iterator parse( + std::format_parse_context& ctx) { + return parse_format_tags(tags, ctx); + } + + template + FmtContext::iterator format(const T& v, FmtContext& ctx) const { + return tags.format(ctx, as_tuple(v)); + } +}; + +template U> +struct format_tag_inherit { + format_tags tags; + + [[nodiscard]] constexpr auto& inner_fmt() { return *this; } + + [[nodiscard]] constexpr const auto& inner_fmt() const { return *this; } + + constexpr std::format_parse_context::iterator parse( + std::format_parse_context& ctx) { + return parse_format_tags(tags, ctx); + } + + template + FmtContext::iterator format(const U& v, FmtContext& ctx) const { + return tags.format(ctx, reinterpret_cast(v)); + } +}; diff --git a/src/core/xml/xml.h b/src/core/xml/xml.h index 9b73d00949..9e16178134 100644 --- a/src/core/xml/xml.h +++ b/src/core/xml/xml.h @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include diff --git a/src/core/xml/xml_io_stream.h b/src/core/xml/xml_io_stream.h index 52e172656a..ff3509c06e 100644 --- a/src/core/xml/xml_io_stream.h +++ b/src/core/xml/xml_io_stream.h @@ -103,28 +103,3 @@ template concept xml_io_parseable = requires(std::span b, std::istream& is) { xml_io_stream::parse(b, is); }; - -//! Explicitly inherit and cast between spans -template T> -struct xml_io_stream_inherit : xml_io_stream { - static void get(std::span b, bifstream* pbifs) - requires(xml_io_binary) - { - xml_io_stream::get(std::span{reinterpret_cast(b.data()), b.size()}, - pbifs); - } - - static void put(std::span b, bofstream* pbofs) - requires(xml_io_binary) - { - xml_io_stream::put( - std::span{reinterpret_cast(b.data()), b.size()}, pbofs); - } - - static void parse(std::span b, std::istream& is) - requires(xml_io_parseable) - { - xml_io_stream::parse(std::span{reinterpret_cast(b.data()), b.size()}, - is); - } -}; diff --git a/src/core/xml/xml_io_stream_aggregate.h b/src/core/xml/xml_io_stream_aggregate.h index 3b9a121096..ab27216fb7 100644 --- a/src/core/xml/xml_io_stream_aggregate.h +++ b/src/core/xml/xml_io_stream_aggregate.h @@ -1,218 +1,10 @@ #pragma once -#include +#include #include "xml_io_base.h" #include "xml_io_stream.h" -template -concept aggregate_0 = std::is_aggregate_v and requires { T{}; }; - -template -concept aggregate_1 = aggregate_0 and requires { T({}); }; - -template -concept aggregate_2 = aggregate_1 and requires { T({}, {}); }; - -template -concept aggregate_3 = aggregate_2 and requires { T({}, {}, {}); }; - -template -concept aggregate_4 = aggregate_3 and requires { T({}, {}, {}, {}); }; - -template -concept aggregate_5 = aggregate_4 and requires { T({}, {}, {}, {}, {}); }; - -template -concept aggregate_6 = - aggregate_5 and requires { T({}, {}, {}, {}, {}, {}); }; - -template -concept aggregate_7 = - aggregate_6 and requires { T({}, {}, {}, {}, {}, {}, {}); }; - -template -concept aggregate_8 = - aggregate_7 and requires { T({}, {}, {}, {}, {}, {}, {}, {}); }; - -template -concept aggregate_9 = - aggregate_8 and requires { T({}, {}, {}, {}, {}, {}, {}, {}, {}); }; - -template -concept aggregate_10 = - aggregate_9 and requires { T({}, {}, {}, {}, {}, {}, {}, {}, {}, {}); }; - -template -concept aggregate_11 = aggregate_10 and requires { - T({}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}); -}; - -template -concept aggregate_12 = aggregate_11 and requires { - T({}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}); -}; - -template -concept aggregate_13 = aggregate_12 and requires { - T({}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}); -}; - -template -concept aggregate_14 = aggregate_13 and requires { - T({}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}); -}; - -template -concept aggregate_15 = aggregate_14 and requires { - T({}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}); -}; - -template -concept aggregate_16 = aggregate_15 and requires { - T({}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}); -}; - -template -concept aggregate_17 = aggregate_16 and requires { - T({}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}); -}; - -template -concept aggregate_18 = aggregate_17 and requires { - T({}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}); -}; - -template -concept aggregate_19 = aggregate_18 and requires { - T({}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}); -}; - -template -concept aggregate_20 = aggregate_19 and requires { - T({}, - {}, - {}, - {}, - {}, - {}, - {}, - {}, - {}, - {}, - {}, - {}, - {}, - {}, - {}, - {}, - {}, - {}, - {}, - {}); -}; - -auto as_tuple(aggregate_0 auto&) { return std::tuple<>(); } -auto as_tuple(aggregate_1 auto& x) { - auto&& [a] = x; - return std::tie(a); -} - -auto as_tuple(aggregate_2 auto& x) { - auto&& [a, b] = x; - return std::tie(a, b); -} - -auto as_tuple(aggregate_3 auto& x) { - auto&& [a, b, c] = x; - return std::tie(a, b, c); -} - -auto as_tuple(aggregate_4 auto& x) { - auto&& [a, b, c, d] = x; - return std::tie(a, b, c, d); -} - -auto as_tuple(aggregate_5 auto& x) { - auto&& [a, b, c, d, e] = x; - return std::tie(a, b, c, d, e); -} - -auto as_tuple(aggregate_6 auto& x) { - auto&& [a, b, c, d, e, f] = x; - return std::tie(a, b, c, d, e, f); -} - -auto as_tuple(aggregate_7 auto& x) { - auto&& [a, b, c, d, e, f, g] = x; - return std::tie(a, b, c, d, e, f, g); -} - -auto as_tuple(aggregate_8 auto& x) { - auto&& [a, b, c, d, e, f, g, h] = x; - return std::tie(a, b, c, d, e, f, g, h); -} - -auto as_tuple(aggregate_9 auto& x) { - auto&& [a, b, c, d, e, f, g, h, i] = x; - return std::tie(a, b, c, d, e, f, g, h, i); -} - -auto as_tuple(aggregate_10 auto& x) { - auto&& [a, b, c, d, e, f, g, h, i, j] = x; - return std::tie(a, b, c, d, e, f, g, h, i, j); -} - -auto as_tuple(aggregate_11 auto& x) { - auto&& [a, b, c, d, e, f, g, h, i, j, k] = x; - return std::tie(a, b, c, d, e, f, g, h, i, j, k); -} - -auto as_tuple(aggregate_12 auto& x) { - auto&& [a, b, c, d, e, f, g, h, i, j, k, l] = x; - return std::tie(a, b, c, d, e, f, g, h, i, j, k, l); -} - -auto as_tuple(aggregate_13 auto& x) { - auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m] = x; - return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m); -} - -auto as_tuple(aggregate_14 auto& x) { - auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m, n] = x; - return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m, n); -} - -auto as_tuple(aggregate_15 auto& x) { - auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o] = x; - return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o); -} - -auto as_tuple(aggregate_16 auto& x) { - auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p] = x; - return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p); -} - -auto as_tuple(aggregate_17 auto& x) { - auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q] = x; - return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q); -} - -auto as_tuple(aggregate_18 auto& x) { - auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r] = x; - return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r); -} - -auto as_tuple(aggregate_19 auto& x) { - auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s] = x; - return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s); -} - -auto as_tuple(aggregate_20 auto& x) { - auto&& [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t] = x; - return std::tie(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t); -} - template struct xml_io_stream_aggregate { constexpr static bool value = false; @@ -222,13 +14,15 @@ template constexpr bool xml_io_stream_aggregate_v = xml_io_stream_aggregate::value; template -concept xml_io_aggregratable = - xml_io_stream_aggregate_v and requires(T a) { as_tuple(a); }; +concept xml_io_aggregratable = xml_io_stream_aggregate_v; template struct xml_io_stream { constexpr static std::string_view type_name = xml_io_stream_name_v; + static_assert(arts_aggregate, + "xml_io_aggregratable types must be arts_aggregate"); + using _lambda = decltype([](auto& v) { return as_tuple(v); }); using ctup_t = std::invoke_result_t<_lambda, const T&>; using mtup_t = std::invoke_result_t<_lambda, T&>; @@ -238,8 +32,8 @@ struct xml_io_stream { const T& t, bofstream* pbofs = nullptr, std::string_view name = ""sv) { - XMLTag tag(type_name, "name", name); - tag.write_to_stream(os); + XMLTag tag(type_name, "name", name); + tag.write_to_stream(os); std::apply( [&os, &pbofs](const Ts&... v) { @@ -247,7 +41,7 @@ struct xml_io_stream { }, as_tuple(t)); - tag.write_to_end_stream(os); + tag.write_to_end_stream(os); } static void read(std::istream& is, T& t, bifstream* pbifs = nullptr) try { diff --git a/src/core/xml/xml_io_stream_inherit.h b/src/core/xml/xml_io_stream_inherit.h new file mode 100644 index 0000000000..5574013732 --- /dev/null +++ b/src/core/xml/xml_io_stream_inherit.h @@ -0,0 +1,70 @@ +#pragma once + +#include +#include + +#include "xml_io_base.h" +#include "xml_io_stream.h" + +using namespace std::literals; + +template U> + requires(not std::same_as) +struct xml_io_stream_inherit { + static constexpr std::string_view type_name = xml_io_stream_name_v; + + struct _no_name_ {}; + static_assert(type_name != xml_io_stream_name_v<_no_name_>, + "xml_io_stream_inherit requires that the base type has a name"); + + static void write(std::ostream& os, + const U& n, + bofstream* pbofs = nullptr, + std::string_view name = ""sv) { + if constexpr (type_name != xml_io_stream_name_v) { + XMLTag tag(type_name, "name", name); + tag.write_to_stream(os); + xml_io_stream::write(os, reinterpret_cast(n), pbofs); + tag.write_to_end_stream(os); + } else { + xml_io_stream::write(os, reinterpret_cast(n), pbofs, name); + } + } + + static void read(std::istream& is, U& n, bifstream* pbifs = nullptr) { + if constexpr (type_name != xml_io_stream_name_v) { + XMLTag tag{}; + tag.read_from_stream(is); + tag.check_name(type_name); + } + + xml_io_stream::read(is, reinterpret_cast(n), pbifs); + + if constexpr (type_name != xml_io_stream_name_v) { + XMLTag tag{}; + tag.read_from_stream(is); + tag.check_end_name(type_name); + } + } + + static void get(std::span b, bifstream* pbifs) + requires(xml_io_binary) + { + xml_io_stream::get(std::span{reinterpret_cast(b.data()), b.size()}, + pbifs); + } + + static void put(std::span b, bofstream* pbofs) + requires(xml_io_binary) + { + xml_io_stream::put( + std::span{reinterpret_cast(b.data()), b.size()}, pbofs); + } + + static void parse(std::span b, std::istream& is) + requires(xml_io_parseable) + { + xml_io_stream::parse(std::span{reinterpret_cast(b.data()), b.size()}, + is); + } +}; diff --git a/src/python_interface/hpy_arts.h b/src/python_interface/hpy_arts.h index c307ef89d6..34e87aa68d 100644 --- a/src/python_interface/hpy_arts.h +++ b/src/python_interface/hpy_arts.h @@ -17,8 +17,8 @@ namespace Python { namespace py = nanobind; using namespace py::literals; -template -void xml_interface(py::class_& c) { +template +void xml_interface(py::class_& c) { c.def( "savexml", [](const T& x, @@ -260,8 +260,8 @@ void value_holder_interface(py::class_>& c) { } } -template -void str_interface(py::class_& c) { +template +void str_interface(py::class_& c) { c.def("__format__", [](const T& x, std::string fmt) { if constexpr (std::formattable) { fmt = std::format("{}{}{}", "{:"sv, fmt, "}"sv); diff --git a/src/python_interface/py_griddedfield.cpp b/src/python_interface/py_griddedfield.cpp index 0d9a1a7904..b1a29dbe02 100644 --- a/src/python_interface/py_griddedfield.cpp +++ b/src/python_interface/py_griddedfield.cpp @@ -197,6 +197,13 @@ void py_griddedfield(py::module_& m) try { m, "ArrayOfNamedGriddedField2"); generic_interface(d2); vector_interface(d2); + + auto vsgf1num = py::bind_vector, + py::rv_policy::reference_internal>( + m, "ArrayOfSortedGriddedField1"); + vsgf1num.doc() = "A list of :class:`~pyarts3.arts.SortedGriddedField1`"; + generic_interface(vsgf1num); + vector_interface(vsgf1num); } catch (std::exception& e) { throw std::runtime_error( std::format("DEV ERROR:\nCannot initialize gridded field\n{}", e.what())); diff --git a/src/python_interface/py_sensor.cpp b/src/python_interface/py_sensor.cpp index 84905031ce..605e3399b2 100644 --- a/src/python_interface/py_sensor.cpp +++ b/src/python_interface/py_sensor.cpp @@ -1,12 +1,18 @@ +#include #include +#include +#include +#include #include #include +#include #include #include #include #include #include #include +#include #include #include @@ -14,9 +20,10 @@ #include #include "hpy_arts.h" +#include "hpy_matpack.h" #include "hpy_numpy.h" #include "hpy_vector.h" -#include "rtepack.h" +#include "xml_io_stream.h" namespace Python { void py_sensor(py::module_& m) try { @@ -263,7 +270,11 @@ Numeric, Vector, or Matrix .def_prop_ro( "poslos", &SensorObsel::poslos_grid, - "Position and line of sight grid\n\n.. :class:`SensorPosLosVector`"); + "Position and line of sight grid\n\n.. :class:`SensorPosLosVector`") + .def("normalize", + &SensorObsel::normalize, + "pol"_a = Stokvec{1., 0., 0., 0.}, + R"(Normalize the weights so that they sum to pol.)"); auto a0 = py::bind_vector, py::rv_policy::reference_internal>( @@ -298,49 +309,7 @@ Numeric, Vector, or Matrix a1.def( "collect_freq_grids", - [](ArrayOfSensorObsel& x) { - if (x.empty()) return; - - const auto tmp = [&x]() { - std::set> freq_set; - for (const auto& obsel : x) freq_set.insert(obsel.f_grid_ptr()); - if (freq_set.size() == 1) - return std::make_pair(true, *freq_set.begin()); - - std::vector fs; - fs.reserve(std::accumulate( - freq_set.begin(), - freq_set.end(), - Size(0), - [](Size a, const auto& b) { return a + b->size(); })); - for (auto& freqs : freq_set) fs.append_range(*freqs); - - stdr::sort(fs); - fs.erase(std::unique(fs.begin(), fs.end()), fs.end()); - return std::make_pair( - false, std::make_shared(std::move(fs))); - }(); - - //! Clang bug workaround - const auto& early = tmp.first; - const auto& freq_ptr = tmp.second; - - if (early) return; - -#pragma omp parallel for if (not arts_omp_in_parallel()) - for (Size i = 0; i < x.size(); i++) { - sensor::SparseStokvecMatrix w(x[i].poslos_grid().size(), - freq_ptr->size()); - for (const auto& we : x[i].weight_matrix()) { - const Index ifreq = stdr::distance( - freq_ptr->begin(), - stdr::lower_bound(*freq_ptr, x[i].f_grid()[we.icol])); - w[we.irow, ifreq] = we.data; - } - - x[i] = {freq_ptr, x[i].poslos_grid_ptr(), std::move(w)}; - } - }, + [](ArrayOfSensorObsel& x) { collect_frequency_grids(x); }, R"(Collect all frequency grids in a single grid. .. tip:: @@ -410,6 +379,30 @@ Numeric, Vector, or Matrix of one channel will affect all channels. )"); + a1.def( + "normalize", + [](ArrayOfSensorObsel& x, Stokvec pol) { + for (auto& obsel : x) obsel.normalize(pol); + }, + R"(Normalize the weights of each obsel so that they sum to pol. + +See :meth:`SensorObsel.normalize` for details. + +.. note:: + The polarization of each channel will be the same. This is a + catch-all method that simply helps when you want something up-and-running. + Each channel must be normalized separately if they have different polarizations. +)", + "pol"_a = Stokvec{1., 0., 0., 0.}); + + a1.def( + "collect", + [](py::object& self) { + self.attr("collect_freq_grids")(); + self.attr("collect_poslos_grids")(); + }, + R"(Calls the collect_freq_grids and collect_poslos_grids methods in said order.)"); + py::class_ smi(m, "SensorMetaInfo"); generic_interface(smi); smi.def_rw( @@ -428,6 +421,456 @@ Numeric, Vector, or Matrix generic_interface(a2); vector_interface(a2); + py::class_ sapf(m, "SensorAntennaPatternField"); + sapf.doc() = + "A 2D gridded antenna response on local zenith and azimuth offsets."; + gridded_data_interface(sapf); + generic_interface(sapf); + + auto sap = py::class_(m, "SensorAntennaPattern"); + sap.doc() = + "Base class for angular antenna responses defined around a bore line of sight."; + sap.def("__call__", + &sensor::AntennaPattern::operator(), + "channel"_a, + "pos"_a, + "bore_los"_a, + R"(Map the antenna pattern onto one observation element.)"); + // generic_interface(sap); -- this should break things. The class is abstract + + auto sgp = py::class_( + m, "SensorGriddedAntennaPattern"); + sgp.doc() = + "A gridded angular antenna response on local zenith and azimuth offsets."; + sgp.def(py::init<>(), "Construct an empty gridded antenna pattern.") + .def( + "__init__", + [](sensor::GriddedAntennaPattern* self, + const sensor::AntennaPatternField& response) { + new (self) sensor::GriddedAntennaPattern{}; + self->data = response; + }, + "response"_a, + "Construct a gridded antenna pattern from a local zenith/azimuth response field.") + .def_rw( + "response", + &sensor::GriddedAntennaPattern::data, + py::rv_policy::reference_internal, + "Local antenna response field.\n\n.. :class:`~pyarts3.arts.SensorAntennaPatternField`"); + generic_interface(sgp); + + auto spp = py::class_( + m, "SensorPencilBeamAntenna"); + spp.doc() = "A 1x1 antenna response centered on the bore line of sight."; + spp.def( + py::init(), + "weight"_a = Stokvec{1.0, 0.0, 0.0, 0.0}, + "Construct a pencil-beam antenna with one Stokes weight at the bore LOS."); + generic_interface(spp); + + auto sga = py::class_( + m, "SensorGaussianAntenna"); + sga.doc() = + "A Gaussian antenna response defined on a local zenith/azimuth grid."; + sga.def(py::init(), + "zen_grid"_a, + "azi_grid"_a, + "zenith_std"_a, + "azimuth_std"_a, + "weight"_a = Stokvec{1.0, 0.0, 0.0, 0.0}, + "Construct a Gaussian antenna response on the supplied local grids."); + generic_interface(sga); + + auto sgairy = py::class_( + m, "SensorGaussianAiryAntenna"); + sgairy.doc() = + "A Gaussianized Airy antenna whose width scales with wavelength and aperture diameter."; + sgairy.def( + py::init(), + "dzen_grid"_a, + "azi_grid"_a, + "aperture_diameter"_a, + "weight"_a = Stokvec{1.0, 0.0, 0.0, 0.0}, + "Construct a Gaussian Airy antenna on the supplied local grids." + " The channel frequency grid must be strictly positive because the" + " current builder path does not carry a separate reference frequency."); + sgairy.def_rw( + "aperture_diameter", + &sensor::GaussianAiryAntenna::aperture_diameter, + "Aperture diameter in the same length units as the zenith grid.\n\n.. :class:`Numeric`"); + sgairy.def_rw( + "weight", + &sensor::GaussianAiryAntenna::weight, + "Stokes weights of the antenna response.\n\n.. :class:`Stokvec`"); + sgairy.def_rw( + "zen_grid", + &sensor::GaussianAiryAntenna::zen_grid, + "Local zenith grid of the antenna response.\n\n.. :class:`ZenGrid`"); + sgairy.def_rw( + "azi_grid", + &sensor::GaussianAiryAntenna::azi_grid, + "Local azimuth grid of the antenna response.\n\n.. :class:`AziGrid`"); + generic_interface(sgairy); + static_assert(arts_xml_ioable); + + auto sch = py::class_(m, "SensorChannel"); + sch.doc() = "Base class for relative spectrometer channel responses."; + generic_interface(sch); + sch.def(py::init_implicit()); + sch.def_prop_ro( + "response", + [](const sensor::Channel& self) -> const SortedGriddedField1& { + return self.channel; + }, + py::rv_policy::reference_internal, + "Relative channel response as a gridded field." + "\n\n.. :class:`~pyarts3.arts.SortedGriddedField1`") + .def_prop_ro("freq_grid", + &sensor::Channel::freq_grid, + py::rv_policy::reference_internal, + "Relative frequency grid.\n\n.. :class:`AscendingGrid`") + .def_prop_ro("weights", + &sensor::Channel::weights, + py::rv_policy::reference_internal, + "Channel weights on the relative frequency grid." + "\n\n.. :class:`Vector`"); + + auto sbox = + py::class_(m, "SensorBoxChannel"); + sbox.doc() = + "A channel with uniform weights across a finite relative-frequency interval."; + sbox.def(py::init(), + "lower"_a, + "upper"_a, + "n"_a, + "Construct a box channel on ``[lower, upper]`` with ``n`` points.") + .def( + py::init(), + "half_width"_a, + "n"_a, + "Construct a symmetric box channel on ``[-half_width, half_width]``.") + .def(py::init(), + "freq_grid"_a, + "Construct a box channel directly from a relative frequency grid."); + generic_interface(sbox); + + auto sdirac = py::class_( + m, "SensorDiracChannel"); + sdirac.doc() = "A single-frequency relative channel."; + sdirac + .def(py::init(), + "frequency"_a, + "Construct a Dirac channel at one relative frequency.") + .def(py::init<>(), "Construct a Dirac channel at 0 relative frequency."); + generic_interface(sdirac); + + auto sgauss = py::class_( + m, "SensorGaussianChannel"); + sgauss.doc() = "A Gaussian relative channel response."; + sgauss + .def(py::init(), + "freq_grid"_a, + "center"_a, + "std"_a, + "Construct a Gaussian channel on a custom relative frequency grid.") + .def( + py::init(), + "center"_a, + "std"_a, + "n"_a, + "m"_a, + "Construct a Gaussian channel on ``center +/- m * std`` with ``n`` points.") + .def(py::init(), + "freq_grid"_a, + "std"_a, + "Construct a zero-centered Gaussian channel on a custom grid.") + .def(py::init(), + "std"_a, + "n"_a, + "m"_a, + "Construct a zero-centered Gaussian channel on ``+/- m * std``."); + generic_interface(sgauss); + + auto sspec = py::class_(m, "SensorSpectrometer"); + sspec.doc() = + "A spectrometer made from channels arranged on one shared relative-frequency grid."; + sspec + .def( + py::init(), + "base_channel"_a, + "freq_offsets"_a, + "Construct a spectrometer by shifting one base channel across the supplied relative-frequency offsets.") + .def(py::init>(), + "channels"_a, + "Construct a spectrometer from explicit channels.") + .def_prop_ro( + "channels", + [](const sensor::Spectrometer& self) + -> const std::vector& { return self.channels; }, + py::rv_policy::reference_internal, + "Spectrometer channels.\n\n.. :class:`list[~pyarts3.arts.SensorChannel]`") + .def( + "is_synced", + &sensor::Spectrometer::is_synced, + "Return whether all channels share the same relative-frequency grid.") + .def( + "__len__", + [](const sensor::Spectrometer& self) { return self.channels.size(); }) + .def( + "__getitem__", + [](const sensor::Spectrometer& self, + Size i) -> const sensor::Channel& { return self.channels.at(i); }, + py::rv_policy::reference_internal) + .def( + "__iter__", + [](const sensor::Spectrometer& self) { + return py::make_iterator(py::type(), + "sensor-spectrometer-iterator", + self.channels.begin(), + self.channels.end()); + }, + py::rv_policy::reference_internal, + "Iterate over the spectrometer channels."); + generic_interface(sspec); + + auto ssb = py::class_(m, "SensorBuilder"); + ssb.doc() = + "Combines channels and an antenna pattern into sensor obsels for paired position/LOS samples."; + ssb.def(py::init<>(), "Construct an empty sensor builder.") + .def( + "__init__", + [](sensor::SensorBuilder* self, + const std::vector& channels, + const sensor::AntennaPattern& antenna) { + new (self) sensor::SensorBuilder{channels, antenna.clone()}; + }, + "channels"_a, + "antenna"_a = sensor::PencilBeamAntenna{}, + "Construct a sensor builder from channels and one antenna pattern (defaults to pencil-beam).") + .def( + "__init__", + [](sensor::SensorBuilder* self, + const sensor::Spectrometer& spectrometer, + const sensor::AntennaPattern& antenna) { + new (self) sensor::SensorBuilder{spectrometer, antenna.clone()}; + }, + "spectrometer"_a, + "antenna"_a = sensor::PencilBeamAntenna{}, + "Construct a sensor builder from a spectrometer and one antenna pattern (defaults to pencil-beam).") + .def( + "__init__", + [](sensor::SensorBuilder* self, + const sensor::AntennaPattern& antenna, + const sensor::Spectrometer& spectrometer, + const sensor::HeterodyneFrequencyRange& backend) { + new (self) + sensor::SensorBuilder{spectrometer, backend, antenna.clone()}; + }, + "antenna"_a, + "spectrometer"_a, + "backend"_a, + "Construct a sensor builder from an antenna, a spectrometer, and a backend frequency response.") + .def_rw( + "channels", + &sensor::SensorBuilder::channels, + "Spectrometer channels.\n\n.. :class:`list[~pyarts3.arts.SensorChannel]`") + .def_prop_rw( + "antenna", + [](const sensor::SensorBuilder& self) + -> std::shared_ptr { + if (not self.antenna) + throw std::runtime_error("SensorBuilder has no antenna pattern"); + return self.antenna->clone(); + }, + [](sensor::SensorBuilder& self, + const sensor::AntennaPattern& antenna) { + self.antenna = antenna.clone(); + }, + py::rv_policy::reference_internal, + "Angular antenna response.\n\n.. :class:`~pyarts3.arts.SensorAntennaPattern`") + .def( + "__call__", + [](const sensor::SensorBuilder& self, + const std::vector& pos, + const std::vector& los) { + if (pos.size() != los.size()) { + throw std::invalid_argument(std::format( + "SensorBuilder expects matching position and LOS counts. Got {} positions and {} LOS values.", + pos.size(), + los.size())); + } + + return self(pos, los); + }, + "pos"_a, + "los"_a, + R"(Build sensor obsels and metadata from paired position and bore-LOS sequences. + +Each ``pos[i]`` is combined with ``los[i]``. The returned obsels are ordered by +geometry first and channel second, and the returned value is +``(measurement_sensor, measurement_sensor_meta)``.)") + .def( + "__call__", + [](const sensor::SensorBuilder& self, + const Vector3 pos, + const Vector2 los) { + return self(std::span{&pos, 1}, std::span{&los, 1}); + }, + "pos"_a, + "los"_a, + R"(Build sensor obsels and metadata from paired position and bore-LOS.)") + .def( + "__call__", + [](const sensor::SensorBuilder& self, + const SensorPosLosVector& poslos) { + std::vector pos(poslos.size()); + std::vector los(poslos.size()); + + for (Size i = 0; i < poslos.size(); ++i) { + pos[i] = poslos[i].pos; + los[i] = poslos[i].los; + } + + return self(pos, los); + }, + "poslos"_a, + "Build sensor obsels and metadata from paired position/LOS samples."); + generic_interface(ssb); + + auto shdfr = py::class_( + m, "SensorHeterodyneFrequencyRange"); + shdfr.def(py::init<>(), + R"(Construct an empty staged heterodyne response. + + Stages can then be applied in sequence using :func:`lowpass`, :func:`highpass`, + :func:`bandpass`, :func:`filter`, and :func:`mix`.)"); + shdfr.def( + py::init(), + "lo"_a, + "bandpass"_a = Vector2{0.0, std::numeric_limits::infinity()}, + R"(Construct a heterodyne response from one ideal bandpass and one LO stage. + + This is shorthand for creating an empty object, applying :func:`bandpass`, and + then applying :func:`mix`.)"); + shdfr.def( + "__init__", + [](sensor::HeterodyneFrequencyRange* v, + const std::vector& clock_frequencies, + const std::vector& bandpasses) { + new (v) sensor::HeterodyneFrequencyRange(clock_frequencies, bandpasses); + }, + "lo"_a, + "bandpasses"_a, + R"(Construct a heterodyne response from a sequence of ideal bandpass and LO stages. + + The sequence is applied as ``bandpass[0] -> lo[0] -> bandpass[1] -> lo[1] -> ...``.)"); + shdfr + .def_ro("global_ranges", + &sensor::HeterodyneFrequencyRange::global_ranges, + "Global frequency range\n\n.. :class:`list[Vector2]`") + .def_ro("local_ranges", + &sensor::HeterodyneFrequencyRange::local_ranges, + "Local frequency range\n\n.. :class:`list[Vector2]`") + .def("lowpass", + &sensor::HeterodyneFrequencyRange::apply_lowpass, + "upper"_a, + "Apply an ideal lowpass filter on the current local frequency axis.") + .def( + "highpass", + &sensor::HeterodyneFrequencyRange::apply_highpass, + "lower"_a, + "Apply an ideal highpass filter on the current local frequency axis.") + .def( + "bandpass", + [](sensor::HeterodyneFrequencyRange& self, const Vector2& bandpass) { + self.apply_bandpass(bandpass); + }, + "bandpass"_a, + "Apply an ideal bandpass filter on the current local frequency axis.") + .def( + "filter", + [](sensor::HeterodyneFrequencyRange& self, + const SortedGriddedField1& bandpass_filter) { + self.apply_bandpass(bandpass_filter); + }, + "bandpass_filter"_a, + R"(Apply a weighted bandpass filter on the current local frequency axis. + + The filter weights are interpreted on the filter's relative frequency grid and + are zero outside that grid.)") + .def("mix", + &sensor::HeterodyneFrequencyRange::apply_mixer, + "lo"_a, + "Apply one heterodyne LO mixing stage.") + .def( + "local_response", + [](const sensor::HeterodyneFrequencyRange& self, + const Vector& f, + Size path_index) { return self.local_response(f, path_index); }, + "f"_a, + "path_index"_a = 0, + "Evaluate one path response on the current local frequency axis.") + .def( + "global_response", + [](const sensor::HeterodyneFrequencyRange& self, + const Vector& f, + Size path_index) { return self.global_response(f, path_index); }, + "f"_a, + "path_index"_a = 0, + "Evaluate one path response on the original real-frequency axis.") + .def( + "channel_response", + [](const sensor::HeterodyneFrequencyRange& self, + const sensor::Channel& channel) { + return sensor::make_bandpass_channels( + self, std::span{&channel, 1}) + .front() + .channel; + }, + "channel"_a, + R"(Compute the real-frequency response for one spectrometer channel. + + The returned gridded field is aggregated across all active mixer paths.)") + .def( + "channel_responses", + [](const sensor::HeterodyneFrequencyRange& self, + const std::vector& channels) { + auto response_channels = sensor::make_bandpass_channels( + self, + std::span{channels.data(), + channels.size()}); + std::vector out; + out.reserve(response_channels.size()); + for (auto& response : response_channels) { + out.push_back(std::move(response.channel)); + } + return out; + }, + "channels"_a, + R"(Compute the real-frequency response for multiple spectrometer channels. + + Each returned gridded field is aggregated across all active mixer paths for the + matching input channel.)") + .def( + "channel_responses", + [](const sensor::HeterodyneFrequencyRange& self, + const sensor::Spectrometer& spectrometer) { + auto response_channels = + sensor::make_bandpass_channels(self, spectrometer); + std::vector out; + out.reserve(response_channels.size()); + for (auto& response : response_channels) { + out.push_back(std::move(response.channel)); + } + return out; + }, + "spectrometer"_a, + R"(Compute the real-frequency response for all channels in a spectrometer. + + Each returned gridded field is aggregated across all active mixer paths for the + matching spectrometer channel.)"); + shdfr.doc() = "A staged heterodyne mixer and filter response builder."; + generic_interface(shdfr); } catch (std::exception& e) { throw std::runtime_error( std::format("DEV ERROR:\nCannot initialize sensors\n{}", e.what())); diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index 6e92ad0b6b..f7a73be6ac 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -171,6 +171,18 @@ add_executable(test_rng test_rng.cc) target_link_libraries(test_rng PUBLIC artscore) target_include_directories(test_rng PRIVATE ${ARTS_SOURCE_DIR}/src) +# #### +add_executable(test_antenna_pattern test_antenna_pattern.cc) +target_link_libraries(test_antenna_pattern PUBLIC artsworkspace) +add_test(NAME "cpp.fast.test_antenna_pattern" COMMAND test_antenna_pattern) +add_dependencies(check-deps test_antenna_pattern) + +# #### +add_executable(test_sensor_builder test_sensor_builder.cc) +target_link_libraries(test_sensor_builder PUBLIC artsworkspace) +add_test(NAME "cpp.fast.test_sensor_builder" COMMAND test_sensor_builder) +add_dependencies(check-deps test_sensor_builder) + # #### add_executable(test_rtepack test_rtepack.cc) target_link_libraries(test_rtepack PUBLIC artstime rtepack rng) diff --git a/src/tests/test_antenna_pattern.cc b/src/tests/test_antenna_pattern.cc new file mode 100644 index 0000000000..9ea5f46cae --- /dev/null +++ b/src/tests/test_antenna_pattern.cc @@ -0,0 +1,220 @@ +#include + +#include + +#include +#include +#include + +namespace { +void assert_stokvec(Stokvec actual, + Stokvec expected, + Numeric tol, + std::string_view name) { + for (Size i = 0; i < 4; ++i) { + ARTS_USER_ERROR_IF(std::abs(actual[i] - expected[i]) > tol, + "{} mismatch at {}: actual {} expected {}", + name, + i, + actual, + expected) + } +} + +Stokvec sum_column_weights(const sensor::Obsel& obsel, Size ifreq) { + Stokvec sum{0.0, 0.0, 0.0, 0.0}; + for (Size isample = 0; isample < obsel.poslos_grid().size(); ++isample) { + sum += obsel.weight_matrix()[isample, ifreq]; + } + return sum; +} + +void test_gaussian_initialization_uses_antenna_frame_offsets() { + const Stokvec peak_weight{2.0, 1.0, 0.0, 0.0}; + const sensor::GaussianAntenna ant{ + ZenGrid{{0.0, 1.0}}, AziGrid{{0.0, 180.0}}, 1.0, 2.0, peak_weight}; + + assert_stokvec(ant.data[0, 0], peak_weight, 1e-12, "gaussian peak weight"); + + const Numeric scale = std::exp(-0.5); + const Stokvec expected = scale * peak_weight; + assert_stokvec(ant.data[1, 0], expected, 1e-12, "gaussian local [1, 0]"); + assert_stokvec(ant.data[1, 1], expected, 1e-12, "gaussian local [1, 180]"); + + const sensor::GaussianAntenna def{ZenGrid{{0.0}}, AziGrid{{0.0}}, 1.0, 1.0}; + assert_stokvec(def.data[0, 0], + {1.0, 0.0, 0.0, 0.0}, + 1e-12, + "default gaussian peak weight"); +} + +Numeric gaussian_airy_expected_gain(Numeric zenith_deg, + Numeric frequency, + Numeric aperture_diameter) { + constexpr Numeric gaussian_airy_hwhm_factor = + 3.8317059702075123156 / Constant::pi; + const Numeric wavelength = Constant::speed_of_light / frequency; + const Numeric hwhm_deg = Conversion::rad2deg( + gaussian_airy_hwhm_factor * wavelength / aperture_diameter); + const Numeric ratio = zenith_deg / hwhm_deg; + + return std::exp(-std::log(2.0) * ratio * ratio); +} + +Numeric gaussian_airy_expected_std(Numeric frequency, + Numeric aperture_diameter) { + constexpr Numeric gaussian_airy_hwhm_factor = + 3.8317059702075123156 / Constant::pi; + const Numeric wavelength = Constant::speed_of_light / frequency; + const Numeric hwhm_deg = Conversion::rad2deg( + gaussian_airy_hwhm_factor * wavelength / aperture_diameter); + + return hwhm_deg / std::sqrt(2.0 * std::log(2.0)); +} + +void test_gaussian_airy_is_frequency_dependent() { + const Stokvec peak_weight{2.0, 0.0, 0.0, 0.0}; + const sensor::GaussianAiryAntenna ant{ + ZenGrid{{0.0, 0.2}}, AziGrid{{0.0}}, 1.0, peak_weight}; + const sensor::BoxChannel channel{AscendingGrid{100.0e9, 200.0e9}}; + + const auto obsel = ant(channel, {600e3, 10.0, 20.0}, {45.0, 30.0}); + + const Numeric low_gain = gaussian_airy_expected_gain(0.2, 100.0e9, 1.0); + const Numeric high_gain = gaussian_airy_expected_gain(0.2, 200.0e9, 1.0); + + const Numeric low_norm = 1.0 + low_gain; + const Numeric high_norm = 1.0 + high_gain; + + assert_stokvec(obsel.weight_matrix()[0, 0], + (0.5 / low_norm) * peak_weight, + 1e-12, + "gaussian airy bore low frequency"); + assert_stokvec(obsel.weight_matrix()[0, 1], + (0.5 / high_norm) * peak_weight, + 1e-12, + "gaussian airy bore high frequency"); + + assert_stokvec(obsel.weight_matrix()[1, 0], + (0.5 * low_gain / low_norm) * peak_weight, + 1e-12, + "gaussian airy off-axis low frequency"); + assert_stokvec(obsel.weight_matrix()[1, 1], + (0.5 * high_gain / high_norm) * peak_weight, + 1e-12, + "gaussian airy off-axis high frequency"); + + assert_stokvec( + sum_column_weights(obsel, 0), 0.5 * peak_weight, 1e-12, "gaussian airy normalized low frequency"); + assert_stokvec( + sum_column_weights(obsel, 1), 0.5 * peak_weight, 1e-12, "gaussian airy normalized high frequency"); + + const auto low_weight = obsel.weight_matrix()[1, 0]; + const auto high_weight = obsel.weight_matrix()[1, 1]; + ARTS_USER_ERROR_IF(low_weight[0] <= high_weight[0], + "Higher frequency Gaussian Airy sample must be narrower") +} + +void test_gaussian_airy_matches_frequency_specific_gaussian_pattern() { + const ZenGrid zen_grid{{0.1, 0.2}}; + const AziGrid azi_grid{{0.0, 0.2}}; + const Stokvec peak_weight{2.0, 0.0, 0.0, 0.0}; + const sensor::GaussianAiryAntenna ant{zen_grid, azi_grid, 1.0, peak_weight}; + const sensor::BoxChannel channel{AscendingGrid{100.0e9, 200.0e9}}; + + const auto obsel = ant(channel, {600e3, 10.0, 20.0}, {45.0, 30.0}); + + for (Size ifreq = 0; ifreq < channel.freq_grid().size(); ++ifreq) { + const Numeric airy_std = + gaussian_airy_expected_std(channel.freq_grid()[ifreq], 1.0); + const sensor::GaussianAntenna gaussian{ + zen_grid, azi_grid, airy_std, airy_std, peak_weight}; + + Numeric normalization = 0.0; + for (Size izen = 0; izen < zen_grid.size(); ++izen) { + for (Size iazi = 0; iazi < azi_grid.size(); ++iazi) { + normalization += gaussian.data[izen, iazi][0] / peak_weight[0]; + } + } + + Size isample = 0; + for (Size izen = 0; izen < zen_grid.size(); ++izen) { + for (Size iazi = 0; iazi < azi_grid.size(); ++iazi) { + assert_stokvec(obsel.weight_matrix()[isample, ifreq], + (channel.weights()[ifreq] / normalization) * + gaussian.data[izen, iazi], + 1e-12, + "gaussian airy per-frequency gaussian match"); + ++isample; + } + } + + assert_stokvec(sum_column_weights(obsel, ifreq), + channel.weights()[ifreq] * peak_weight, + 1e-12, + "gaussian airy per-frequency column normalization"); + } +} + + void test_degenerate_azimuth_ring_is_collapsed() { + const sensor::GaussianAntenna gaussian{ + ZenGrid{{0.0, 1.0}}, AziGrid{{0.0, 120.0, 240.0}}, 1.0, 1.0}; + const auto gaussian_obsel = + gaussian(sensor::DiracChannel{100.0e9}, {600e3, 10.0, 20.0}, {45.0, 30.0}); + + ARTS_USER_ERROR_IF(gaussian_obsel.poslos_grid().size() != 4, + "Expected zero-zenith azimuth collapse to leave 4 LOS samples, got {}", + gaussian_obsel.poslos_grid().size()) + assert_stokvec(gaussian_obsel.weight_matrix()[0, 0], + {3.0, 0.0, 0.0, 0.0}, + 1e-12, + "collapsed gaussian zero-zenith ring weight"); + + const sensor::GaussianAiryAntenna airy{ + ZenGrid{{0.0, 0.2}}, AziGrid{{0.0, 120.0, 240.0}}, 1.0}; + const auto airy_obsel = + airy(sensor::DiracChannel{100.0e9}, {600e3, 10.0, 20.0}, {45.0, 30.0}); + const Numeric gain = gaussian_airy_expected_gain(0.2, 100.0e9, 1.0); + const Numeric retained_scale = 3.0 * (1.0 + gain) / (1.0 + 3.0 * gain); + + ARTS_USER_ERROR_IF(airy_obsel.poslos_grid().size() != 4, + "Expected Gaussian Airy zero-zenith azimuth collapse to leave 4 LOS samples, got {}", + airy_obsel.poslos_grid().size()) + assert_stokvec(airy_obsel.weight_matrix()[0, 0], + {retained_scale / (3.0 * (1.0 + gain)), 0.0, 0.0, 0.0}, + 1e-6, + "collapsed gaussian airy center is renormalized with retained points"); + assert_stokvec(airy_obsel.weight_matrix()[1, 0], + {retained_scale * gain / (3.0 * (1.0 + gain)), 0.0, 0.0, 0.0}, + 1e-6, + "collapsed gaussian airy off-axis weight is renormalized uniformly"); + assert_stokvec(sum_column_weights(airy_obsel, 0), + {1.0, 0.0, 0.0, 0.0}, + 1e-6, + "collapsed gaussian airy column renormalizes to unity"); + } + +void test_gaussian_airy_rejects_nonpositive_frequencies() { + const sensor::GaussianAiryAntenna ant{ZenGrid{{0.0}}, AziGrid{{0.0}}, 1.0}; + + bool threw = false; + try { + static_cast( + ant(sensor::DiracChannel{}, {600e3, 10.0, 20.0}, {0.0, 0.0})); + } catch (const std::runtime_error&) { + threw = true; + } + + ARTS_USER_ERROR_IF(not threw, + "Gaussian Airy antenna must reject nonpositive channel frequencies") +} +} // namespace + +int main() { + test_gaussian_initialization_uses_antenna_frame_offsets(); + test_gaussian_airy_is_frequency_dependent(); + test_gaussian_airy_matches_frequency_specific_gaussian_pattern(); + test_degenerate_azimuth_ring_is_collapsed(); + test_gaussian_airy_rejects_nonpositive_frequencies(); + return 0; +} \ No newline at end of file diff --git a/src/tests/test_sensor_builder.cc b/src/tests/test_sensor_builder.cc new file mode 100644 index 0000000000..34e5965ba4 --- /dev/null +++ b/src/tests/test_sensor_builder.cc @@ -0,0 +1,200 @@ +#include +#include +#include + +#include +#include +#include +#include + +#include "arts_constants.h" + +namespace { +void assert_close(Numeric actual, + Numeric expected, + Numeric tol, + std::string_view name) { + ARTS_USER_ERROR_IF(std::abs(actual - expected) > tol, + "{} mismatch: actual {} expected {}", + name, + actual, + expected) +} + +void assert_stokvec(Stokvec actual, + Stokvec expected, + Numeric tol, + std::string_view name) { + for (Size i = 0; i < 4; ++i) { + ARTS_USER_ERROR_IF(std::abs(actual[i] - expected[i]) > tol, + "{} mismatch at {}: actual {} expected {}", + name, + i, + actual, + expected) + } +} + +Numeric gaussian_airy_expected_gain(Numeric zenith_deg, + Numeric frequency, + Numeric aperture_diameter) { + constexpr Numeric gaussian_airy_hwhm_factor = + 3.8317059702075123156 / Constant::pi; + const Numeric wavelength = Constant::speed_of_light / frequency; + const Numeric hwhm_deg = Conversion::rad2deg(gaussian_airy_hwhm_factor * + wavelength / aperture_diameter); + const Numeric ratio = zenith_deg / hwhm_deg; + + return std::exp(-Constant::ln_2 * ratio * ratio); +} + +void test_sensor_builder_returns_meta_per_geometry() { + sensor::SensorBuilder builder( + {sensor::BoxChannel{AscendingGrid{100.0, 101.0}}, + sensor::DiracChannel{200.0}}, + std::make_shared()); + + const std::array pos{{{600e3, 10.0, 20.0}, {601e3, 11.0, 21.0}}}; + const std::array los{{{20.0, 30.0}, {40.0, 50.0}}}; + + const auto [obsels, meta] = builder(pos, los); + + ARTS_USER_ERROR_IF( + obsels.size() != 4, "Expected 4 obsels, got {}", obsels.size()) + ARTS_USER_ERROR_IF( + meta.size() != 2, "Expected 2 meta entries, got {}", meta.size()) + ARTS_USER_ERROR_IF(meta[0].count() != 2 or meta[1].count() != 2, + "Each meta block must describe one 2-channel geometry") + + const auto& gf0 = std::get(meta[0].data); + const auto& gf1 = std::get(meta[1].data); + + ARTS_USER_ERROR_IF( + gf0.gridname<0>() != "channel" or gf1.gridname<0>() != "channel", + "SensorBuilder meta grid must be the channel axis") + assert_close(gf0.grid<0>()[0], 0.0, 0.0, "meta[0] channel index 0"); + assert_close(gf0.grid<0>()[1], 1.0, 0.0, "meta[0] channel index 1"); + + ARTS_USER_ERROR_IF(obsels[0].poslos_grid()[0].pos != pos[0], + "First geometry position mismatch") + ARTS_USER_ERROR_IF(obsels[2].poslos_grid()[0].pos != pos[1], + "Second geometry position mismatch") + ARTS_USER_ERROR_IF(not obsels[0].same_poslos(obsels[1]), + "Obsels from the same geometry must share poslos") + ARTS_USER_ERROR_IF(not obsels[2].same_poslos(obsels[3]), + "Obsels from the same geometry must share poslos") + ARTS_USER_ERROR_IF(not obsels[0].same_freqs(obsels[2]), + "Obsels for the same channel must share frequencies") + ARTS_USER_ERROR_IF(not obsels[1].same_freqs(obsels[3]), + "Obsels for the same channel must share frequencies") +} + +void test_sensor_builder_rejects_mismatched_geometry_counts() { + sensor::SensorBuilder builder({sensor::DiracChannel{}}, + std::make_shared()); + + const std::array pos{{{600e3, 10.0, 20.0}}}; + const std::array los{{{20.0, 30.0}, {40.0, 50.0}}}; + + bool threw = false; + try { + static_cast(builder(pos, los)); + } catch (const std::runtime_error&) { + threw = true; + } + + ARTS_USER_ERROR_IF( + not threw, + "SensorBuilder must reject mismatching position and LOS counts") +} + +void test_sensor_builder_uses_gaussian_airy_frequency_dependence() { + const Stokvec peak_weight{2.0, 0.0, 0.0, 0.0}; + sensor::SensorBuilder builder( + {sensor::BoxChannel{AscendingGrid{100.0e9, 200.0e9}}}, + std::make_shared( + ZenGrid{{0.0, 0.2}}, AziGrid{{0.0}}, 1.0, peak_weight)); + + const std::array pos{{{600e3, 10.0, 20.0}}}; + const std::array los{{{45.0, 30.0}}}; + + const auto [obsels, meta] = builder(pos, los); + + ARTS_USER_ERROR_IF( + obsels.size() != 1, "Expected 1 obsel, got {}", obsels.size()) + ARTS_USER_ERROR_IF( + meta.size() != 1, "Expected 1 meta entry, got {}", meta.size()) + + const Numeric low_gain = gaussian_airy_expected_gain(0.2, 100.0e9, 1.0); + const Numeric high_gain = gaussian_airy_expected_gain(0.2, 200.0e9, 1.0); + const Numeric low_norm = 1.0 + low_gain; + const Numeric high_norm = 1.0 + high_gain; + + assert_stokvec(obsels[0].weight_matrix()[1, 0], + 0.5 * low_gain * peak_weight / low_norm, + 1e-12, + "builder gaussian airy off-axis low frequency"); + assert_stokvec(obsels[0].weight_matrix()[1, 1], + 0.5 * high_gain * peak_weight / high_norm, + 1e-12, + "builder gaussian airy off-axis high frequency"); +} + +void test_unflatten_updates_shared_poslos_grids() { + auto freq_grid = std::make_shared(AscendingGrid{100.0}); + + SensorPosLosVector poslos_grid(1); + poslos_grid[0] = {.pos = {600e3, 10.0, 20.0}, .los = {30.0, 40.0}}; + auto shared_poslos = + std::make_shared(std::move(poslos_grid)); + + sensor::SparseStokvecMatrix weights(1, 1); + weights[0, 0] = {1.0, 0.0, 0.0, 0.0}; + + ArrayOfSensorObsel obsels(2); + obsels[0] = {freq_grid, shared_poslos, weights}; + obsels[1] = {freq_grid, shared_poslos, weights}; + + const auto original_poslos = obsels[0].poslos_grid_ptr(); + const Vector zeniths{15.0}; + + unflatten(obsels, zeniths, obsels[0], SensorKeyType::zen); + + ARTS_USER_ERROR_IF(obsels[0].poslos_grid_ptr() == original_poslos, + "unflatten must replace the shared poslos grid") + ARTS_USER_ERROR_IF(not obsels[0].same_poslos(obsels[1]), + "obsels sharing a geometry must still share it after unflatten") + assert_close(obsels[0].poslos_grid()[0].los[0], + 15.0, + 0.0, + "updated zenith for first obsel"); + assert_close(obsels[1].poslos_grid()[0].los[0], + 15.0, + 0.0, + "updated zenith for second obsel"); + assert_close(obsels[0].poslos_grid()[0].los[1], + 40.0, + 0.0, + "azimuth must remain unchanged"); + + const auto original_freq = obsels[0].f_grid_ptr(); + const Vector freqs{101.0}; + + unflatten(obsels, freqs, obsels[0], SensorKeyType::freq); + + ARTS_USER_ERROR_IF(obsels[0].f_grid_ptr() == original_freq, + "unflatten must replace the shared frequency grid") + ARTS_USER_ERROR_IF(not obsels[0].same_freqs(obsels[1]), + "obsels sharing frequencies must still share them after unflatten") + assert_close(obsels[0].f_grid()[0], 101.0, 0.0, "updated frequency for first obsel"); + assert_close(obsels[1].f_grid()[0], 101.0, 0.0, "updated frequency for second obsel"); +} +} // namespace + +int main() { + test_sensor_builder_returns_meta_per_geometry(); + test_sensor_builder_rejects_mismatched_geometry_counts(); + test_sensor_builder_uses_gaussian_airy_frequency_dependence(); + test_unflatten_updates_shared_poslos_grids(); + return 0; +} \ No newline at end of file diff --git a/tests/core/sensor/30cm-antenna.py b/tests/core/sensor/30cm-antenna.py new file mode 100644 index 0000000000..72d626664a --- /dev/null +++ b/tests/core/sensor/30cm-antenna.py @@ -0,0 +1,36 @@ +import pyarts3 as pa +import numpy as np + +NZA = 20 +NAZ = 30 +POS = [1, 0, 0] +LOS = [45, 30] + +zen_grid = np.geomspace(1, 3, NZA) - 1 +azi_grid = np.linspace(0, 360, NAZ + 1)[:-1] + +ant1 = pa.arts.SensorGaussianAiryAntenna(zen_grid, azi_grid, 30e-2, "Iv") +ant2 = pa.arts.SensorGaussianAiryAntenna(zen_grid + 1e-3, azi_grid, 30e-2, "I") +ch = pa.arts.SensorBoxChannel(1e9, 1e12, 1001) + +obsel1 = ant1(ch, POS, LOS) +obsel2 = ant2(ch, POS, LOS) + +poslos1 = np.asarray(obsel1.poslos) +poslos2 = np.asarray(obsel2.poslos) +weights1 = np.asarray(obsel1.weight_matrix.dense) +weights2 = np.asarray(obsel2.weight_matrix.dense) + +assert np.allclose(weights1.sum(), 2.0), \ + "Should be normalized to 2 (since polarized)" +assert np.allclose(weights2.sum(), 1.0), \ + "Should be normalized to 1 (since non-polarized)" +assert len(poslos1) == (NZA * NAZ - NAZ + 1), \ + "Should have one poslos per non-zero DZA, but only one for zero DZA - has zero DZA" +assert len(poslos2) == NZA * NAZ, \ + "Should have one poslos per non-zero DZA, but only one for zero DZA - has no zero DZA" + +ant = pa.arts.SensorGaussianAiryAntenna(zen_grid, azi_grid, 30e-2, "I") + +pa.plots.plot(obsel1, point_spread=True) +pa.plots.plot(obsel2, point_spread=True, frame='local', type='contour', levels=10) diff --git a/tests/core/sensor/600ghz-upper-sideband-sensor.py b/tests/core/sensor/600ghz-upper-sideband-sensor.py new file mode 100644 index 0000000000..2e3f6fd2c8 --- /dev/null +++ b/tests/core/sensor/600ghz-upper-sideband-sensor.py @@ -0,0 +1,56 @@ +import numpy as np +import pyarts3 as pyarts +import time + +start_time = time.time() + +NCHANNELS = 10_000 +IF_LOW = 5.5e9 +IF_HIGH = 6.5e9 +LO = 600e9 +APERTURE_DIAMETER = 30e-2 +POS = [600e3, 0.0, 0.0] +LOS = [45.0, 0.0] +ZEN = np.linspace(0.0, 0.3, 7) +AZI = np.linspace(0.0, 360.0, 9)[:-1] + + +if_offsets = pyarts.arts.AscendingGrid(np.linspace(IF_LOW, IF_HIGH, NCHANNELS)) +print(f"if_offsets: {time.time() - start_time:.2f} seconds") + +spectrometer = pyarts.arts.SensorSpectrometer( + pyarts.arts.SensorDiracChannel(), + if_offsets, +) +print(f"spectrometer: {time.time() - start_time:.2f} seconds") + +backend = pyarts.arts.SensorHeterodyneFrequencyRange( + LO, + np.array([LO + IF_LOW, LO + IF_HIGH]), +) +print(f"backend: {time.time() - start_time:.2f} seconds") + +antenna = pyarts.arts.SensorGaussianAiryAntenna( + ZEN, + AZI, + APERTURE_DIAMETER, + "I", +) +print(f"antenna: {time.time() - start_time:.2f} seconds") + +sensor = pyarts.arts.SensorBuilder(antenna, spectrometer, backend) +print(f"sensor: {time.time() - start_time:.2f} seconds") + +ws = pyarts.Workspace() + +ws.measurement_sensor, ws.measurement_sensor_meta = sensor(POS, LOS) + +assert len(ws.measurement_sensor.unique_freq_grids()) == 1 +assert len(ws.measurement_sensor[0].f_grid) == NCHANNELS + +np.testing.assert_allclose( + np.asarray(ws.measurement_sensor[0].f_grid)[[0, -1]], + np.array([LO + IF_LOW, LO + IF_HIGH]), + atol=0.0, + rtol=0.0, +) diff --git a/tests/core/sensor/heterodyne_frequency_response.py b/tests/core/sensor/heterodyne_frequency_response.py new file mode 100644 index 0000000000..a1c35b9e36 --- /dev/null +++ b/tests/core/sensor/heterodyne_frequency_response.py @@ -0,0 +1,330 @@ +import numpy as np +import pyarts3 as pyarts + + +def db_to_lin(db): + return 10.0 ** (db / 10.0) + + +def ranges_as_lists(ranges): + return [list(map(float, pair)) for pair in ranges] + + +def affine_from_ranges(global_range, local_range): + g0, g1 = map(float, global_range) + l0, l1 = map(float, local_range) + slope = (g1 - g0) / (l1 - l0) + intercept = g0 - slope * l0 + return [intercept, slope] + + +def assert_close(name, got, expected, atol=1e-12): + got = np.asarray(got, dtype=float) + expected = np.asarray(expected, dtype=float) + print(f"{name}: got = {got}") + print(f"{name}: expected = {expected}") + np.testing.assert_allclose(got, expected, atol=atol, rtol=0.0) + + +def inspect_case(name, selector, expected_global, expected_local, expected_affine): + print(f"\n=== {name} ===") + + got_global = ranges_as_lists(selector.global_ranges) + got_local = ranges_as_lists(selector.local_ranges) + got_affine = [ + affine_from_ranges(global_range, local_range) + for global_range, local_range in zip(got_global, got_local) + ] + + print("Affine form per path: global = intercept + slope * local") + assert_close(f"{name} global_ranges", got_global, expected_global) + assert_close(f"{name} local_ranges", got_local, expected_local) + assert_close(f"{name} affine", got_affine, expected_affine) + + +def make_triangular_filter(): + filt = pyarts.arts.SortedGriddedField1() + filt.grids = (pyarts.arts.AscendingGrid(np.array([5.0, 10.0, 15.0])),) + filt.data = pyarts.arts.Vector(np.array([0.0, 1.0, 0.0])) + filt.gridnames = ("frequency",) + filt.dataname = "triangular" + return filt + + +def make_asymmetric_sideband_filter(): + filt = pyarts.arts.SortedGriddedField1() + filt.grids = ( + pyarts.arts.AscendingGrid(np.array([7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0])), + ) + filt.data = pyarts.arts.Vector( + np.array( + [ + db_to_lin(-10.0), + db_to_lin(-10.0), + db_to_lin(-10.0), + db_to_lin(0.0), + db_to_lin(1.0), + db_to_lin(5.0), + db_to_lin(10.0), + ] + ) + ) + filt.gridnames = ("frequency",) + filt.dataname = "asymmetric_sidebands" + return filt + + +def test_lowpass_then_lo_then_box_spectrometer(): + selector = pyarts.arts.SensorHeterodyneFrequencyRange() + print("test_lowpass_then_lo_then_box_spectrometer") + selector.lowpass(15.0) + print("Apply lowpass keeping 0-15:", selector) + selector.mix(16.0) + print("Apply LO at 16 :", selector) + + inspect_case( + "lowpass -> LO", + selector, + expected_global=[[15.0, 0.0]], + expected_local=[[1.0, 16.0]], + expected_affine=[[16.0, -1.0]], + ) + + channel = pyarts.arts.SensorBoxChannel(2.6, 4.6, 5) + response = selector.channel_response(channel) + + assert_close( + "lowpass -> LO channel points", + np.array(response.grids[0]), + np.array([11.4, 11.9, 12.4, 12.9, 13.4]), + ) + assert_close( + "lowpass -> LO channel weights", + np.array(response.data), + np.full(5, 0.2), + ) + + +def test_default_positive_range_survives_mixes(): + selector = pyarts.arts.SensorHeterodyneFrequencyRange() + selector.mix(10.0) + + assert_close( + "default positive range after one LO global_ranges", + ranges_as_lists(selector.global_ranges), + [[10.0, np.inf], [10.0, 0.0]], + ) + assert_close( + "default positive range after one LO local_ranges", + ranges_as_lists(selector.local_ranges), + [[0.0, np.inf], [0.0, 10.0]], + ) + + selector.mix(1.0) + + assert_close( + "default positive range after two LOs global_ranges", + ranges_as_lists(selector.global_ranges), + [[11.0, np.inf], [11.0, 10.0], [9.0, 0.0], [9.0, 10.0]], + ) + assert_close( + "default positive range after two LOs local_ranges", + ranges_as_lists(selector.local_ranges), + [[0.0, np.inf], [0.0, 1.0], [0.0, 9.0], [0.0, 1.0]], + ) + + response = selector.channel_response(pyarts.arts.SensorDiracChannel(0.25)) + + assert_close( + "default positive range dirac points after two LOs", + np.array(response.grids[0]), + np.array([8.75, 9.25, 10.75, 11.25]), + ) + assert_close( + "default positive range dirac weights after two LOs", + np.array(response.data), + np.full(4, 0.25), + ) + + +def test_bandpass_lo_bandpass_lo_chain(): + selector = pyarts.arts.SensorHeterodyneFrequencyRange() + selector.bandpass(np.array([5.0, 15.0])) + selector.mix(3.0) + selector.bandpass(np.array([3.0, 7.0])) + selector.mix(6.0) + + inspect_case( + "bandpass -> LO -> bandpass -> LO", + selector, + expected_global=[[9.0, 10.0], [9.0, 6.0]], + expected_local=[[0.0, 1.0], [0.0, 3.0]], + expected_affine=[[9.0, 1.0], [9.0, -1.0]], + ) + + channel = pyarts.arts.SensorBoxChannel(0.0, 1.0, 3) + response = selector.channel_response(channel) + + assert_close( + "bandpass -> LO -> bandpass -> LO channel points", + np.array(response.grids[0]), + np.array([8.0, 8.5, 9.0, 9.5, 10.0]), + ) + assert_close( + "bandpass -> LO -> bandpass -> LO channel weights", + np.array(response.data), + np.full(5, 1.0 / 6.0), + ) + + +def test_ideal_split_about_lo(): + selector = pyarts.arts.SensorHeterodyneFrequencyRange() + selector.bandpass(np.array([5.0, 15.0])) + selector.mix(10.0) + + inspect_case( + "ideal split about LO", + selector, + expected_global=[[10.0, 15.0], [10.0, 5.0]], + expected_local=[[0.0, 5.0], [0.0, 5.0]], + expected_affine=[[10.0, 1.0], [10.0, -1.0]], + ) + + assert_close( + "ideal split path 0 global response at 12", + selector.global_response(np.array([12.0]), 0), + np.array([1.0]), + ) + assert_close( + "ideal split path 1 global response at 8", + selector.global_response(np.array([8.0]), 1), + np.array([1.0]), + ) + + +def test_weighted_split_about_lo(): + selector = pyarts.arts.SensorHeterodyneFrequencyRange() + selector.filter(make_triangular_filter()) + selector.mix(10.0) + + inspect_case( + "weighted split about LO", + selector, + expected_global=[[10.0, 15.0], [10.0, 5.0]], + expected_local=[[0.0, 5.0], [0.0, 5.0]], + expected_affine=[[10.0, 1.0], [10.0, -1.0]], + ) + + assert_close( + "weighted split path 0 local response at 2", + selector.local_response(np.array([2.0]), 0), + np.array([0.6]), + ) + assert_close( + "weighted split path 1 local response at 2", + selector.local_response(np.array([2.0]), 1), + np.array([0.6]), + ) + + channel = pyarts.arts.SensorDiracChannel(2.0) + response = selector.channel_response(channel) + + assert_close( + "weighted split channel points", + np.array(response.grids[0]), + np.array([8.0, 12.0]), + ) + assert_close( + "weighted split channel weights", + np.array(response.data), + np.array([0.3, 0.3]), + ) + + +def test_zero_frequency_is_not_double_counted(): + selector = pyarts.arts.SensorHeterodyneFrequencyRange() + selector.bandpass(np.array([5.0, 15.0])) + selector.mix(10.0) + + response = selector.channel_response(pyarts.arts.SensorDiracChannel(0.0)) + + assert_close( + "zero-frequency channel point", + np.array(response.grids[0]), + np.array([10.0]), + ) + assert_close( + "zero-frequency channel weight", + np.array(response.data), + np.array([0.5]), + ) + + +def test_overlapping_sidebands_with_asymmetric_bandpass(): + selector = pyarts.arts.SensorHeterodyneFrequencyRange() + selector.filter(make_asymmetric_sideband_filter()) + selector.mix(10.0) + + inspect_case( + "overlapping sidebands with asymmetric bandpass", + selector, + expected_global=[[10.0, 13.0], [10.0, 7.0]], + expected_local=[[0.0, 3.0], [0.0, 3.0]], + expected_affine=[[10.0, 1.0], [10.0, -1.0]], + ) + + upper_expected = np.array([db_to_lin(1.0), db_to_lin(5.0), db_to_lin(10.0)]) + lower_expected = np.full(3, db_to_lin(-10.0)) + + assert_close( + "overlapping sidebands upper-path local gains", + selector.local_response(np.array([1.0, 2.0, 3.0]), 0), + upper_expected, + ) + assert_close( + "overlapping sidebands lower-path local gains", + selector.local_response(np.array([1.0, 2.0, 3.0]), 1), + lower_expected, + ) + + responses = selector.channel_responses( + [ + pyarts.arts.SensorDiracChannel(1.0), + pyarts.arts.SensorDiracChannel(2.0), + pyarts.arts.SensorDiracChannel(3.0), + ] + ) + + expected_points = [ + np.array([9.0, 11.0]), + np.array([8.0, 12.0]), + np.array([7.0, 13.0]), + ] + expected_weights = [ + np.array([db_to_lin(-10.0) / 2.0, db_to_lin(1.0) / 2.0]), + np.array([db_to_lin(-10.0) / 2.0, db_to_lin(5.0) / 2.0]), + np.array([db_to_lin(-10.0) / 2.0, db_to_lin(10.0) / 2.0]), + ] + + for index, response in enumerate(responses): + assert_close( + f"overlapping sidebands channel {index} points", + np.array(response.grids[0]), + expected_points[index], + ) + assert_close( + f"overlapping sidebands channel {index} weights", + np.array(response.data), + expected_weights[index], + ) + + +test_lowpass_then_lo_then_box_spectrometer() +test_default_positive_range_survives_mixes() +test_bandpass_lo_bandpass_lo_chain() +test_ideal_split_about_lo() +test_weighted_split_about_lo() +test_zero_frequency_is_not_double_counted() +test_overlapping_sidebands_with_asymmetric_bandpass() + +print("\nAll heterodyne frequency-response checks passed.") \ No newline at end of file diff --git a/tests/core/sensor/sensor_builder.py b/tests/core/sensor/sensor_builder.py new file mode 100644 index 0000000000..d03804559c --- /dev/null +++ b/tests/core/sensor/sensor_builder.py @@ -0,0 +1,84 @@ +import numpy as np +import pyarts3 as pyarts + + +def assert_close(name, got, expected, atol=1e-12): + got = np.asarray(got, dtype=float) + expected = np.asarray(expected, dtype=float) + print(f"{name}: got = {got}") + print(f"{name}: expected = {expected}") + np.testing.assert_allclose(got, expected, atol=atol, rtol=0.0) + + +def test_sensor_builder_from_sequences(): + builder = pyarts.arts.SensorBuilder( + [ + pyarts.arts.SensorBoxChannel(pyarts.arts.AscendingGrid([100.0, 101.0])), + pyarts.arts.SensorDiracChannel(200.0), + ], + pyarts.arts.SensorPencilBeamAntenna(), + ) + + pos = [ + pyarts.arts.Vector3([600e3, 10.0, 20.0]), + pyarts.arts.Vector3([601e3, 11.0, 21.0]), + ] + los = [ + pyarts.arts.Vector2([20.0, 30.0]), + pyarts.arts.Vector2([40.0, 50.0]), + ] + + obsels, meta = builder(pos, los) + + assert len(obsels) == 4 + assert len(meta) == 2 + assert meta[0].count == 2 + assert meta[1].count == 2 + assert_close("meta[0] channel grid", meta[0].data.grids[0], [0.0, 1.0]) + assert_close("first pos", obsels[0].poslos[0][:3], pos[0]) + assert_close("first los", obsels[0].poslos[0][3:], los[0]) + assert_close("second geometry pos", obsels[2].poslos[0][:3], pos[1]) + assert_close("second geometry los", obsels[2].poslos[0][3:], los[1]) + assert_close("first channel grid", obsels[0].f_grid, [100.0, 101.0]) + assert_close("second channel grid", obsels[1].f_grid, [200.0]) + assert_close("first channel weights", + obsels[0].weight_matrix.dense[0, 0], + [0.5, 0.0, 0.0, 0.0]) + + +def test_sensor_builder_from_poslos_vector_and_length_guard(): + builder = pyarts.arts.SensorBuilder( + [pyarts.arts.SensorDiracChannel()], + pyarts.arts.SensorPencilBeamAntenna(), + ) + + poslos = pyarts.arts.SensorPosLosVector() + poslos.value = np.array( + [ + [600e3, 10.0, 20.0, 20.0, 30.0], + [601e3, 11.0, 21.0, 40.0, 50.0], + ] + ) + + obsels, meta = builder(poslos) + + assert len(obsels) == 2 + assert len(meta) == 2 + assert meta[0].count == 1 + assert meta[1].count == 1 + assert_close("poslos builder first los", obsels[0].poslos[0][3:], [20.0, 30.0]) + assert_close("poslos builder second los", obsels[1].poslos[0][3:], [40.0, 50.0]) + + try: + builder([pyarts.arts.Vector3([600e3, 10.0, 20.0])], [ + pyarts.arts.Vector2([20.0, 30.0]), + pyarts.arts.Vector2([40.0, 50.0]), + ]) + except ValueError: + pass + else: + raise AssertionError("SensorBuilder must reject mismatching position/LOS lengths") + + +test_sensor_builder_from_sequences() +test_sensor_builder_from_poslos_vector_and_length_guard() \ No newline at end of file