From 0bba1991844b5d4d1ede3622f18d623fdfdca029 Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Tue, 27 Oct 2026 19:00:03 +0100 Subject: [PATCH 1/3] Consolidate radar calculations in and around rtepack --- src/core/CMakeLists.txt | 3 + src/core/options/arts_options.cc | 47 ++ src/core/radar/CMakeLists.txt | 4 + src/core/radar/radar.cc | 268 +++++++ src/core/radar/radar.h | 177 ++++ src/core/rtepack/CMakeLists.txt | 1 + src/core/rtepack/rtepack.h | 1 + src/core/rtepack/rtepack_radar.cc | 174 ++++ src/core/rtepack/rtepack_radar.h | 72 ++ src/core/rtepack/rtepack_scattering.cc | 150 ---- src/core/rtepack/rtepack_scattering.h | 47 -- src/core/rtepack/test/CMakeLists.txt | 5 + src/core/rtepack/test/test_rtepack_radar.cpp | 172 ++++ src/m_montecarlo_radar.cc | 285 +++---- src/m_radar.cc | 799 ++++++------------- src/python_interface/CMakeLists.txt | 1 + src/python_interface/py_module.cpp | 2 + src/python_interface/py_radar.cpp | 27 + src/tests/test_ai_port_regressions.cc | 9 +- src/workspace_groups.cpp | 8 + src/workspace_methods.cpp | 152 ++-- src/workspace_variables.cpp | 10 + tests/core/scat/mc_radar.py | 38 +- tests/core/scat/radar_single_scattering.py | 12 + 24 files changed, 1502 insertions(+), 962 deletions(-) create mode 100644 src/core/radar/CMakeLists.txt create mode 100644 src/core/radar/radar.cc create mode 100644 src/core/radar/radar.h create mode 100644 src/core/rtepack/rtepack_radar.cc create mode 100644 src/core/rtepack/rtepack_radar.h create mode 100644 src/core/rtepack/test/test_rtepack_radar.cpp create mode 100644 src/python_interface/py_radar.cpp diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 383bc9bddc..e66616289e 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -36,6 +36,8 @@ add_subdirectory(predefined) add_subdirectory(quantum) +add_subdirectory(radar) + add_subdirectory(rng) add_subdirectory(rtepack) @@ -81,6 +83,7 @@ target_link_libraries(artscore INTERFACE physics predef quantum + radar rng rtepack sensor diff --git a/src/core/options/arts_options.cc b/src/core/options/arts_options.cc index fb7df0ebd9..ec8eaef301 100644 --- a/src/core/options/arts_options.cc +++ b/src/core/options/arts_options.cc @@ -199,6 +199,53 @@ runs out of damping without an acceptable step. A starting cost above Value{"Lookup", "Lookup table antenna"}}, }); + opts.emplace_back(EnumeratedOption{ + .name = "RadarUnit", + .desc = R"(The unit of a simulated radar return. + +The native return is the gate-averaged, attenuated backscatter coefficient. +The equivalent reflectivity factor is + +.. math:: + + Z_e = \frac{4 \cdot 10^{18} \lambda^4}{\pi^4 |K|^2} \beta, + +where :math:`\lambda` is the wavelength [m], :math:`|K|^2` is the reference +dielectric factor, and :math:`\beta` is the native return. +)", + .values_and_desc = {Value{"unit", "1", "Backscatter coefficient [1 / (m sr)]"}, + Value{"Ze", "ze", "Equivalent reflectivity factor [mm :math:`^{6}` / m :math:`^{3}`]"}, + Value{"dBZe", "dbze", "Equivalent reflectivity factor in decibels, :math:`10 \\log_{10} Z_e`"}}, + }); + + opts.emplace_back(EnumeratedOption{ + .name = "RadarRangeMode", + .desc = R"(The coordinate of radar range-gate limits. +)", + .values_and_desc = {Value{"Altitude", "Altitude [m]"}, + Value{"Distance", "One-way geometric distance from the radar [m]"}, + Value{"RoundTripTime", "Round-trip time from the radar, including path group refractive indices [s]"}, + Value{"Legacy", "The ARTS2 rule: Altitude if the largest limit exceeds 1, otherwise RoundTripTime"}}, + }); + + opts.emplace_back(EnumeratedOption{ + .name = "RadarAuxiliary", + .desc = R"(Auxiliary quantities of a radar simulation. +)", + .values_and_desc = { + Value{"RadiativeBackground", + "Radiative background", + "Type of the path end: 0 for space, 1 for the surface, and 2 otherwise"}, + Value{"Backscattering", "backscattering", "Gate-averaged backscatter without attenuation, in the unit of the measurement"}, + Value{"AbsSpeciesExtinction", + "Abs species extinction", + "Gate-averaged extinction by absorbing species [1 / m]"}, + Value{"ParticleExtinction", + "Particle extinction", + "Gate-averaged particle extinction, after any extinction scaling [1 / m]"}, + }, + }); + opts.emplace_back(EnumeratedOption{ .name = "PredefinedSensor", .desc = diff --git a/src/core/radar/CMakeLists.txt b/src/core/radar/CMakeLists.txt new file mode 100644 index 0000000000..bf218d574b --- /dev/null +++ b/src/core/radar/CMakeLists.txt @@ -0,0 +1,4 @@ +add_library(radar STATIC radar.cc) + +target_link_libraries(radar PUBLIC jacobian path rtepack scattering) +target_include_directories(radar PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/core/radar/radar.cc b/src/core/radar/radar.cc new file mode 100644 index 0000000000..4210bd9314 --- /dev/null +++ b/src/core/radar/radar.cc @@ -0,0 +1,268 @@ +#include "radar.h" + +#include +#include + +#include +#include +#include +#include +#include + +namespace radar { +Numeric liebe93_k2(const Numeric frequency, const Numeric temperature) { + if (temperature < Constant::temperature_at_0c - 40.0 or temperature > Constant::temperature_at_0c + 100.0) + arts::user_error("Automatic radar k2 requires {} K <= ze_tref <= {} K; got {} K", + Constant::temperature_at_0c - 40.0, + Constant::temperature_at_0c + 100.0, + temperature); + if (frequency < 10e9 or frequency > 1000e9) + arts::user_error("Automatic radar k2 uses Liebe-93 and requires 10 GHz <= frequency <= 1000 GHz; got {} Hz", + frequency); + + const Numeric theta = 1.0 - 300.0 / temperature; + const Numeric e0 = 77.66 - 103.3 * theta; + const Numeric e1 = 0.0671 * e0; + const Numeric f1 = 20.2 + 146.0 * theta + 316.0 * theta * theta; + const Numeric e2 = 3.52; + const Numeric f2 = 39.8 * f1; + const Complex ifGHz{0.0, frequency / 1e9}; + const Complex n = std::sqrt(e2 + (e1 - e2) / (1.0 - ifGHz / f2) + (e0 - e1) / (1.0 - ifGHz / f1)); + const Complex n2 = n * n; + return std::norm((n2 - 1.0) / (n2 + 2.0)); +} + +Numeric ze_factor(const Numeric frequency, const Numeric ze_tref, const Numeric k2) { + const Numeric dielectric = k2 < 0.0 ? liebe93_k2(frequency, ze_tref) : k2; + if (dielectric <= 0.0) arts::user_error("Radar k2 must be positive, or negative to request Liebe-93; got {}", k2); + const Numeric wavelength = Constant::speed_of_light / frequency; + return 4e18 * std::pow(wavelength, 4) / (std::pow(Constant::pi, 4) * dielectric); +} + +Numeric unit_factor(const RadarUnit unit, const Numeric frequency, const Numeric ze_tref, const Numeric k2) { + switch (unit) { + case RadarUnit::unit: return 1.0; + case RadarUnit::Ze: + case RadarUnit::dBZe: return ze_factor(frequency, ze_tref, k2); + } + std::unreachable(); +} + +RadarRangeMode resolve(const RadarRangeMode mode, const Numeric largest_limit) { + if (mode != RadarRangeMode::Legacy) return mode; + return largest_limit > 1.0 ? RadarRangeMode::Altitude : RadarRangeMode::RoundTripTime; +} + +PathRange path_range(const ArrayOfPropagationPathPoint& path, const Vector3& radar_pos, const Vector2& ellipsoid) { + const Size np = path.size(); + PathRange out{ + .layer = np > 1 ? path::distance(path, ellipsoid) : Vector{}, .distance = Vector(np), .time = Vector(np)}; + if (np == 0) return out; + + out.distance[0] = path::distance(radar_pos, path.front().pos, ellipsoid); + out.time[0] = 2.0 * out.distance[0] / Constant::speed_of_light; + for (Size ip = 1; ip < np; ++ip) { + out.distance[ip] = out.distance[ip - 1] + out.layer[ip - 1]; + out.time[ip] = + out.time[ip - 1] + out.layer[ip - 1] * (path[ip - 1].ngroup + path[ip].ngroup) / Constant::speed_of_light; + } + return out; +} + +Vector range_coordinate(const PathRange& range, const ArrayOfPropagationPathPoint& path, const RadarRangeMode mode) { + switch (mode) { + case RadarRangeMode::Altitude: { + Vector out(path.size()); + for (Size ip = 0; ip < path.size(); ++ip) out[ip] = path[ip].pos[0]; + return out; + } + case RadarRangeMode::Distance: return range.distance; + case RadarRangeMode::RoundTripTime: return range.time; + case RadarRangeMode::Legacy: arts::user_error("The Legacy radar range mode must be resolved before use"); + } + std::unreachable(); +} + +void add_gate_samples(std::vector& samples, + const ConstVectorView& coordinate, + const Numeric lo, + const Numeric hi, + const Size gate) { + const Numeric width = hi - lo; + for (Size i = 1; i < coordinate.size(); i++) { + const Numeric c0 = coordinate[i - 1]; + const Numeric c1 = coordinate[i]; + const Numeric a = std::max(lo, std::min(c0, c1)); + const Numeric b = std::min(hi, std::max(c0, c1)); + if (not(a < b)) continue; + samples.push_back( + {.layer = i, .x = (std::midpoint(a, b) - c0) / (c1 - c0), .weight = (b - a) / width, .gate = gate}); + } +} + +Size sort_samples(std::vector& samples) { + stdr::stable_sort(samples, {}, &rtepack::radar_sample::layer); + return samples.empty() ? 0 : samples.back().layer + 1; +} + +Vector gate_coverage(const std::span& samples, const Size ngate) { + Vector out(ngate, 0.0); + for (const auto& sample : samples) out[sample.gate] += sample.weight; + return out; +} + +Vector gate_sum(const std::span& samples, const ConstVectorView& level, const Size ngate) { + Vector out(ngate, 0.0); + for (const auto& [layer, x, weight, gate] : samples) + out[gate] += weight * ((1.0 - x) * level[layer - 1] + x * level[layer]); + return out; +} + +namespace { +Numeric normalized_delta_aa(Numeric aa) { + while (aa < -180.0) aa += 360.0; + while (aa > 180.0) aa -= 360.0; + return aa; +} + +bool at_pole(const Numeric za) { return std::abs(std::sin(za * Constant::pi / 180.0)) < 1e-12; } +} // namespace + +Optics optics(const ArrayOfPropagationPathPoint& path, + const Size np, + const AtmField& atm_field, + const AscendingGrid& freq_grid, + const ArrayOfScatteringSpecies& scattering_species, + const JacobianTargets& jac_targets, + const Numeric pext_scaling, + const GasModel& gas) { + const Size nf = freq_grid.size(); + const Size nq = jac_targets.target_count(); + + if (np > path.size()) arts::user_error("Cannot evaluate {} radar levels of a path with {} points", np, path.size()); + + Optics out; + out.K_out.resize(np); + out.K_ret.resize(np); + out.dK_out.resize(np); + out.dK_ret.resize(np); + out.Z.resize(np, MuelmatVector(nf, Muelmat{})); + out.dZ.resize(np, MuelmatMatrix(nq, nf, Muelmat{})); + out.gas_extinction.resize(nf, np); + out.particle_extinction.resize(nf, np); + if (np == 0) return out; + + const ArrayOfAtmPoint atm_path = + forward_atm_path(ArrayOfPropagationPathPoint(path.begin(), path.begin() + np), atm_field); + + for (Size ip = 0; ip < np; ++ip) { + // Propagation-path LOS is the propagation direction toward the receiver; + // its mirror is the transmitted propagation direction. ARO phase + // matrices take the latter as incident and the former as scattered. + const PropagationPathPoint& receive_point = path[ip]; + PropagationPathPoint transmit_point = path[ip]; + transmit_point.los = path::mirror(receive_point.los); + + PropmatVector gas_out, gas_ret; + PropmatMatrix gas_out_jac, gas_ret_jac; + gas(gas_out, gas_out_jac, transmit_point, atm_path[ip]); + gas(gas_ret, gas_ret_jac, receive_point, atm_path[ip]); + if (gas_out.size() != nf or gas_ret.size() != nf or static_cast(gas_out_jac.nrows()) != nq or + static_cast(gas_ret_jac.nrows()) != nq or static_cast(gas_out_jac.ncols()) != nf or + static_cast(gas_ret_jac.ncols()) != nf) + arts::user_error( + "The gas model must give {} propagation matrices and a [{}, {}] Jacobian, got {}, {:B,}, {}, and {:B,}", + nf, + nq, + nf, + gas_out.size(), + gas_out_jac.shape(), + gas_ret.size(), + gas_ret_jac.shape()); + + // Azimuth is undefined at zenith and nadir. Use its canonical value at + // either pole to avoid a 0*NaN rotation without changing the scattering + // angle or the physical polarization basis. + const Vector2& transmit_los = transmit_point.los; + const Vector2& receive_los = receive_point.los; + const Vector delta_aa{at_pole(transmit_los[0]) or at_pole(receive_los[0]) + ? 0.0 + : normalized_delta_aa(receive_los[1] - transmit_los[1])}; + const Vector transmit_za{transmit_los[0]}; + const Vector receive_za{receive_los[0]}; + const auto transmit_za_grid = + std::make_shared(scattering::IrregularZenithAngleGrid(transmit_za)); + const auto receive_za_grid = + std::make_shared(scattering::IrregularZenithAngleGrid(receive_za)); + + const auto bulk_out = scattering_species.get_bulk_scattering_properties_aro_gridded( + atm_path[ip], freq_grid, transmit_za, delta_aa, receive_za_grid); + const auto bulk_ret = scattering_species.get_bulk_scattering_properties_aro_gridded( + atm_path[ip], freq_grid, receive_za, Vector{0.0}, transmit_za_grid); + if (not bulk_out.phase_matrix.has_value()) + arts::user_error("Radar single scattering requires phase matrices from all scattering species"); + + const auto& phase = *bulk_out.phase_matrix; + out.K_out[ip] = gas_out; + out.K_ret[ip] = gas_ret; + for (Size iv = 0; iv < nf; ++iv) { + const Propmat particle_out = bulk_out.extinction_matrix[0, iv, 0].expand(); + out.K_out[ip][iv] += pext_scaling * particle_out; + out.K_ret[ip][iv] += pext_scaling * bulk_ret.extinction_matrix[0, iv, 0].expand(); + out.Z[ip][iv] = phase[0, iv, 0, 0, 0]; + out.gas_extinction[iv, ip] = gas_out[iv].A(); + out.particle_extinction[iv, ip] = pext_scaling * particle_out.A(); + } + + out.dK_out[ip] = std::move(gas_out_jac); + out.dK_ret[ip] = std::move(gas_ret_jac); + for (const auto& target : jac_targets.atm) { + const auto dbulk_out = scattering_species.get_bulk_scattering_properties_aro_gridded_derivative( + atm_path[ip], freq_grid, transmit_za, delta_aa, receive_za_grid, target.type); + const auto dbulk_ret = scattering_species.get_bulk_scattering_properties_aro_gridded_derivative( + atm_path[ip], freq_grid, receive_za, Vector{0.0}, transmit_za_grid, target.type); + if (not dbulk_out.phase_matrix) arts::user_error("Radar derivatives require scattering phase matrices"); + + const auto& dphase = *dbulk_out.phase_matrix; + for (Size iv = 0; iv < nf; ++iv) { + out.dK_out[ip][target.target_pos, iv] += pext_scaling * dbulk_out.extinction_matrix[0, iv, 0].expand(); + out.dK_ret[ip][target.target_pos, iv] += pext_scaling * dbulk_ret.extinction_matrix[0, iv, 0].expand(); + out.dZ[ip][target.target_pos, iv] = dphase[0, iv, 0, 0, 0]; + } + } + } + + return out; +} + +void single_scattering(StokvecMatrix& I, + StokvecTensor4& dI, + const Stokvec& transmitted, + const Optics& optics, + const ConstVectorView& layer, + const std::span& samples, + const Size ngate, + const TransmittanceOption option) { + const Size np = optics.size(); + const Size nf = optics.gas_extinction.nrows(); + const Size nq = np == 0 ? 0 : optics.dK_out.front().nrows(); + + if (np > layer.size() + 1) + arts::user_error("The radar optics have {} levels but there are only {} layers", np, layer.size()); + + I.resize(ngate, nf); + dI.resize(ngate, np, nq, nf); + rtepack::radar_single_scattering(I, + dI, + transmitted, + optics.K_out, + optics.dK_out, + optics.K_ret, + optics.dK_ret, + optics.Z, + optics.dZ, + layer[Range(0, np == 0 ? 0 : np - 1)], + samples, + option); +} +} // namespace radar diff --git a/src/core/radar/radar.h b/src/core/radar/radar.h new file mode 100644 index 0000000000..fd6e822cb0 --- /dev/null +++ b/src/core/radar/radar.h @@ -0,0 +1,177 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +using ArrayOfRadarAuxiliary = Array; + +namespace radar { +/** The liquid-water dielectric factor |K|^2 of the Liebe-93 permittivity model. + * + * @param frequency The frequency [Hz] + * @param temperature The temperature [K] + * @return |K|^2 + */ +[[nodiscard]] Numeric liebe93_k2(Numeric frequency, Numeric temperature); + +/** The factor that converts backscatter [1/(m sr)] to Ze [mm^6/m^3]. + * + * @param frequency The frequency [Hz] + * @param ze_tref The temperature of the Liebe-93 dielectric factor [K] + * @param k2 The dielectric factor |K|^2, negative selects Liebe-93 + * @return The conversion factor + */ +[[nodiscard]] Numeric ze_factor(Numeric frequency, Numeric ze_tref, Numeric k2); + +/** The factor that converts backscatter to the linear part of the unit. + * + * dBZe uses the Ze factor, the logarithm is left to the caller. + */ +[[nodiscard]] Numeric unit_factor(RadarUnit unit, Numeric frequency, Numeric ze_tref, Numeric k2); + +/** Resolve the Legacy mode from the largest range limit, other modes are returned as is */ +[[nodiscard]] RadarRangeMode resolve(RadarRangeMode mode, Numeric largest_limit); + +//! The ranges of the points of a radar path +struct PathRange { + //! The geometric thickness of every layer [np - 1] + Vector layer; + + //! The one-way distance from the radar [np] + Vector distance; + + //! The round-trip time from the radar, including the path group refractive indices [np] + Vector time; +}; + +/** The ranges of a radar path. + * + * The first point of the path need not be at the radar, for instance when + * the radar is outside the atmosphere. The range to it is then the geometric + * distance in vacuum. + * + * @param path The propagation path, starting closest to the radar + * @param radar_pos The position of the radar + * @param ellipsoid The reference ellipsoid + */ +[[nodiscard]] PathRange path_range(const ArrayOfPropagationPathPoint& path, + const Vector3& radar_pos, + const Vector2& ellipsoid); + +/** The coordinate of the range-gate limits at every path point. + * + * @param range The ranges of the path + * @param path The propagation path + * @param mode The coordinate, Legacy must already be resolved + */ +[[nodiscard]] Vector range_coordinate(const PathRange& range, + const ArrayOfPropagationPathPoint& path, + RadarRangeMode mode); + +/** Append the quadrature of a range gate. + * + * Every layer that overlaps the gate gets one sample at the midpoint of the + * overlap. The weights are normalized by the full gate width, so the parts + * of a gate without any path contribute nothing. + * + * @param samples The samples to extend + * @param coordinate The range coordinate at every path point + * @param lo The lower gate limit + * @param hi The upper gate limit + * @param gate The output index of the gate + */ +void add_gate_samples( + std::vector& samples, const ConstVectorView& coordinate, Numeric lo, Numeric hi, Size gate); + +/** Sort the samples by layer. + * + * @return The number of levels the samples need + */ +[[nodiscard]] Size sort_samples(std::vector& samples); + +/** The sum of the sample weights of every gate, which is the covered fraction of the gate */ +[[nodiscard]] Vector gate_coverage(const std::span& samples, Size ngate); + +/** The weighted sum of level values, linearly interpolated to the samples, of every gate */ +[[nodiscard]] Vector gate_sum(const std::span& samples, + const ConstVectorView& level, + Size ngate); + +/** The gas propagation matrices on the frequency grid, and their Jacobian [nq, nf], at a path point */ +using GasModel = std::function; + +//! The optical properties of the levels of a radar path +struct Optics { + //! Propagation matrices away from and towards the radar [np][nf] + std::vector K_out, K_ret; + + //! Their derivatives [np][nq, nf] + std::vector dK_out, dK_ret; + + //! The backscatter matrix [np][nf] + std::vector Z; + + //! Its derivative [np][nq, nf] + std::vector dZ; + + //! The extinction away from the radar by absorbing species and by (scaled) particles [nf, np] + Matrix gas_extinction, particle_extinction; + + [[nodiscard]] Size size() const { return K_out.size(); } +}; + +/** The optical properties of the first np points of a radar path. + * + * The scattering properties are evaluated in the monostatic geometry of + * every path point. Derivatives are with regards to the atmospheric targets + * of jac_targets, the gas model provides all others. + * + * @param path The propagation path, starting closest to the radar + * @param np The number of levels to evaluate + * @param atm_field The atmosphere + * @param freq_grid The frequency grid + * @param scattering_species The scattering species + * @param jac_targets The Jacobian targets + * @param pext_scaling The scaling of particle extinction + * @param gas The gas propagation matrices + */ +[[nodiscard]] Optics optics(const ArrayOfPropagationPathPoint& path, + Size np, + const AtmField& atm_field, + const AscendingGrid& freq_grid, + const ArrayOfScatteringSpecies& scattering_species, + const JacobianTargets& jac_targets, + Numeric pext_scaling, + const GasModel& gas); + +/** The single-scattering return of every gate, and its derivatives. + * + * Wraps rtepack::radar_single_scattering for the levels of optics. + * + * @param[out] I The return of every gate [ngate, nf] + * @param[out] dI The derivatives [ngate, np, nq, nf] + * @param transmitted The transmitted Stokes vector + * @param optics The optical properties of the levels + * @param layer The layer thicknesses of the path, at least optics.size() - 1 + * @param samples The gate samples, sorted by layer + * @param ngate The number of gates + * @param option The transmittance option + */ +void single_scattering(StokvecMatrix& I, + StokvecTensor4& dI, + const Stokvec& transmitted, + const Optics& optics, + const ConstVectorView& layer, + const std::span& samples, + Size ngate, + TransmittanceOption option); +} // namespace radar diff --git a/src/core/rtepack/CMakeLists.txt b/src/core/rtepack/CMakeLists.txt index 30cf67fedd..48c161d6a4 100644 --- a/src/core/rtepack/CMakeLists.txt +++ b/src/core/rtepack/CMakeLists.txt @@ -4,6 +4,7 @@ add_library(rtepack STATIC rtepack_multitype.cc rtepack_polarization.cc rtepack_propagation_matrix.cc + rtepack_radar.cc rtepack_rtestep.cc rtepack_scattering.cc rtepack_source.cc diff --git a/src/core/rtepack/rtepack.h b/src/core/rtepack/rtepack.h index 15da52dae9..0a01541cb0 100644 --- a/src/core/rtepack/rtepack.h +++ b/src/core/rtepack/rtepack.h @@ -9,6 +9,7 @@ #include "rtepack_multitype.h" #include "rtepack_polarization.h" #include "rtepack_propagation_matrix.h" +#include "rtepack_radar.h" #include "rtepack_rtestep.h" #include "rtepack_scattering.h" #include "rtepack_source.h" diff --git a/src/core/rtepack/rtepack_radar.cc b/src/core/rtepack/rtepack_radar.cc new file mode 100644 index 0000000000..a3697d4028 --- /dev/null +++ b/src/core/rtepack/rtepack_radar.cc @@ -0,0 +1,174 @@ +#include "rtepack_radar.h" + +#include + +#include + +#include "rtepack_multitype.h" +#include "rtepack_transmission.h" + +namespace rtepack { +namespace { +constexpr bool uses_magnus(const TransmittanceOption option) { + switch (option) { + case TransmittanceOption::constant: + case TransmittanceOption::linsrc: + case TransmittanceOption::linprop: return false; + case TransmittanceOption::magop: + case TransmittanceOption::magop_linsrc: return true; + } + std::unreachable(); +} + +//! Transmission over a layer from start to end, and its derivatives with regards to either end +struct layer_transmission { + propmat start; + propmat end; + Numeric r; + bool magnus; + tran state; + muelmat T; + + layer_transmission(const propmat &k_start, const propmat &k_end, const Numeric r_, const TransmittanceOption option) + : start(k_start), + end(k_end), + r(r_), + magnus(uses_magnus(option)), + state(magnus ? tran{k_start, k_end, r_, MagnusOperator::magnus} : tran{k_start, k_end, r_}), + T(state()) {} + + [[nodiscard]] muelmat d_start(const propmat &dk) const { + return magnus ? state.magnus_deriv(T, start, end, dk, r, 0.0, true) : state.deriv(T, -0.5 * r * dk); + } + + [[nodiscard]] muelmat d_end(const propmat &dk) const { + return magnus ? state.magnus_deriv(T, start, end, dk, r, 0.0, false) : state.deriv(T, -0.5 * r * dk); + } +}; +} // namespace + +void radar_single_scattering(stokvec_matrix_view I, + stokvec_tensor4_view dI, + const stokvec &transmitted, + const std::span &K_out, + const std::span &dK_out, + const std::span &K_ret, + const std::span &dK_ret, + const std::span &Z, + const std::span &dZ, + const ConstVectorView &r, + const std::span &samples, + const TransmittanceOption option) { + const Size np = K_out.size(); + const Size ngate = I.nrows(); + const Size nf = I.ncols(); + const Size nq = dI.nrows(); + + if (static_cast(dI.nbooks()) != ngate or static_cast(dI.npages()) != np or + static_cast(dI.ncols()) != nf) + arts::user_error("Radar return shapes do not agree: I: {:B,}, dI: {:B,}, np: {}", I.shape(), dI.shape(), np); + if (K_ret.size() != np or Z.size() != np or dK_out.size() != np or dK_ret.size() != np or dZ.size() != np or + r.size() + 1 != std::max(np, 1)) + arts::user_error( + "Radar level sizes do not agree: K_out: {}, K_ret: {}, Z: {}, dK_out: {}, dK_ret: {}, dZ: {}, r: {}", + np, + K_ret.size(), + Z.size(), + dK_out.size(), + dK_ret.size(), + dZ.size(), + r.size()); + for (const auto &sample : samples) + if (sample.gate >= ngate or sample.layer == 0 or sample.layer >= np) + arts::user_error("Radar sample in layer {} for gate {} is outside {} levels and {} gates", + sample.layer, + sample.gate, + np, + ngate); + + I = stokvec{}; + dI = stokvec{}; + if (samples.empty()) return; + + // Tf * s and Tr up to the near level of the current layer, and their + // derivatives with regards to every level. + stokvec_vector Tfs(nf, transmitted); + muelmat_vector Tr(nf, muelmat::id()); + stokvec_tensor3 dTfs(np, nq, nf, stokvec{}); + muelmat_tensor3 dTr(np, nq, nf, muelmat{}); + + Size is = 0; + for (Size i = 1; i < np; i++) { + for (; is < samples.size() and samples[is].layer == i; is++) { + const auto &[layer, x, w, gate] = samples[is]; + const Numeric rx = x * r[i - 1]; + + for (Size iv = 0; iv < nf; iv++) { + const propmat ko = (1.0 - x) * K_out[i - 1][iv] + x * K_out[i][iv]; + const propmat kr = (1.0 - x) * K_ret[i - 1][iv] + x * K_ret[i][iv]; + const muelmat Zx = (1.0 - x) * Z[i - 1][iv] + x * Z[i][iv]; + + const layer_transmission out{K_out[i - 1][iv], ko, rx, option}; + const layer_transmission ret{kr, K_ret[i - 1][iv], rx, option}; + + const stokvec tfs = out.T * Tfs[iv]; + const stokvec zs = Zx * tfs; + const muelmat Trx = Tr[iv] * ret.T; + I[gate, iv] += w * (Trx * zs); + + if (nq == 0) continue; + + // Through the cumulative transmissions of all previous layers + const muelmat M = Trx * Zx * out.T; + const stokvec u = ret.T * zs; + for (Size j = 0; j < i; j++) + for (Size iq = 0; iq < nq; iq++) dI[gate, j, iq, iv] += w * (M * dTfs[j, iq, iv] + dTr[j, iq, iv] * u); + + // Through the partial layer, whose inner end is interpolated + for (Size iq = 0; iq < nq; iq++) { + const propmat &dko0 = dK_out[i - 1][iq, iv]; + const propmat &dko1 = dK_out[i][iq, iv]; + const propmat &dkr0 = dK_ret[i - 1][iq, iv]; + const propmat &dkr1 = dK_ret[i][iq, iv]; + + const muelmat dout0 = out.d_start(dko0) + out.d_end((1.0 - x) * dko0); + const muelmat dout1 = out.d_end(x * dko1); + const muelmat dret0 = ret.d_start((1.0 - x) * dkr0) + ret.d_end(dkr0); + const muelmat dret1 = ret.d_start(x * dkr1); + + dI[gate, i - 1, iq, iv] += w * (Tr[iv] * (dret0 * zs) + Trx * ((1.0 - x) * dZ[i - 1][iq, iv] * tfs) + + Trx * (Zx * (dout0 * Tfs[iv]))); + dI[gate, i, iq, iv] += + w * (Tr[iv] * (dret1 * zs) + Trx * (x * dZ[i][iq, iv] * tfs) + Trx * (Zx * (dout1 * Tfs[iv]))); + } + } + } + + if (is == samples.size()) break; + + for (Size iv = 0; iv < nf; iv++) { + const layer_transmission out{K_out[i - 1][iv], K_out[i][iv], r[i - 1], option}; + const layer_transmission ret{K_ret[i][iv], K_ret[i - 1][iv], r[i - 1], option}; + + for (Size j = 0; j < i; j++) { + for (Size iq = 0; iq < nq; iq++) { + dTfs[j, iq, iv] = out.T * dTfs[j, iq, iv]; + dTr[j, iq, iv] = dTr[j, iq, iv] * ret.T; + } + } + + for (Size iq = 0; iq < nq; iq++) { + dTfs[i - 1, iq, iv] += out.d_start(dK_out[i - 1][iq, iv]) * Tfs[iv]; + dTfs[i, iq, iv] += out.d_end(dK_out[i][iq, iv]) * Tfs[iv]; + dTr[i - 1, iq, iv] += Tr[iv] * ret.d_end(dK_ret[i - 1][iq, iv]); + dTr[i, iq, iv] += Tr[iv] * ret.d_start(dK_ret[i][iq, iv]); + } + + Tfs[iv] = out.T * Tfs[iv]; + Tr[iv] = Tr[iv] * ret.T; + } + } + + if (is != samples.size()) arts::user_error("Radar samples must be sorted by layer"); +} +} // namespace rtepack diff --git a/src/core/rtepack/rtepack_radar.h b/src/core/rtepack/rtepack_radar.h new file mode 100644 index 0000000000..405c2dc6a1 --- /dev/null +++ b/src/core/rtepack/rtepack_radar.h @@ -0,0 +1,72 @@ +#pragma once + +#include +#include + +#include + +#include "rtepack_mueller_matrix.h" +#include "rtepack_propagation_matrix.h" +#include "rtepack_stokes_vector.h" + +namespace rtepack { +/** A quadrature point of a radar range gate inside one layer. + * + * The layer spans level layer - 1, closest to the radar, to level layer. + */ +struct radar_sample { + //! The far level of the layer (>= 1) + Size layer; + + //! Fractional position from level layer - 1 towards level layer, in [0, 1] + Numeric x; + + //! Quadrature weight + Numeric weight; + + //! The output index of the range gate + Size gate; +}; + +/** Single-scattering monostatic radar return of a stack of layers. + * + * Level 0 is closest to the radar. The return at a point is Tr * Z * Tf * s, + * where Tf is the ordered transmission from the radar to the point, Tr the + * ordered transmission back, Z the backscatter matrix, and s the transmitted + * Stokes vector. The outgoing and return propagation matrices need neither + * be equal nor commute. + * + * Inside a layer, the propagation and backscatter matrices are linear in the + * position. The transmission to a point inside a layer uses the same + * exponential as the full layer, so the Magnus operator also applies to + * partial layers. There is no source term, so the linear-source options use + * the transmission of their base option. + * + * The derivatives are with regards to the level quantities only. + * + * @param[out] I The return of every gate [ngate, nf] + * @param[out] dI The derivative of every gate [ngate, np, nq, nf] + * @param[in] transmitted The transmitted Stokes vector + * @param[in] K_out The propagation matrix away from the radar [np][nf] + * @param[in] dK_out The derivative of K_out [np][nq, nf] + * @param[in] K_ret The propagation matrix towards the radar [np][nf] + * @param[in] dK_ret The derivative of K_ret [np][nq, nf] + * @param[in] Z The backscatter matrix [np][nf] + * @param[in] dZ The derivative of Z [np][nq, nf] + * @param[in] r The layer thicknesses [np - 1] + * @param[in] samples The gate samples, sorted by layer + * @param[in] option The transmittance option + */ +void radar_single_scattering(stokvec_matrix_view I, + stokvec_tensor4_view dI, + const stokvec &transmitted, + const std::span &K_out, + const std::span &dK_out, + const std::span &K_ret, + const std::span &dK_ret, + const std::span &Z, + const std::span &dZ, + const ConstVectorView &r, + const std::span &samples, + const TransmittanceOption option); +} // namespace rtepack diff --git a/src/core/rtepack/rtepack_scattering.cc b/src/core/rtepack/rtepack_scattering.cc index 81089f90db..e95a8c4255 100644 --- a/src/core/rtepack/rtepack_scattering.cc +++ b/src/core/rtepack/rtepack_scattering.cc @@ -8,156 +8,6 @@ #include "rtepack_stokes_vector.h" namespace rtepack { -Array bulk_backscatter(const ConstTensor5View &Pe, const ConstMatrixView &pnd) { - assert(Pe.ncols() == 4 and Pe.nrows() == 4); - const Index nv = Pe.npages(); - const Index np = Pe.nbooks(); - const Index ne = Pe.nshelves(); - - Array aotm(np, muelmat_vector(nv)); - if (ne == 0) return aotm; - - for (Index ip = 0; ip < np; ip++) { - for (Index iv = 0; iv < nv; iv++) { - auto m = Pe[0, ip, iv, joker, joker]; - aotm[ip][iv] = pnd[0, ip] * muelmat(m[0, 0], - m[0, 1], - m[0, 2], - m[0, 3], - m[1, 0], - m[1, 1], - m[1, 2], - m[1, 3], - m[2, 0], - m[2, 1], - m[2, 2], - m[2, 3], - m[3, 0], - m[3, 1], - m[3, 2], - m[3, 3]); - } - } - - for (Index ie = 1; ie < ne; ie++) { - for (Index ip = 0; ip < np; ip++) { - for (Index iv = 0; iv < nv; iv++) { - auto m = Pe[ie, ip, iv, joker, joker]; - aotm[ip][iv] += pnd[ie, ip] * muelmat(m[0, 0], - m[0, 1], - m[0, 2], - m[0, 3], - m[1, 0], - m[1, 1], - m[1, 2], - m[1, 3], - m[2, 0], - m[2, 1], - m[2, 2], - m[2, 3], - m[3, 0], - m[3, 1], - m[3, 2], - m[3, 3]); - } - } - } - return aotm; -} - -Array bulk_backscatter_derivative(const ConstTensor5View &Pe, const ArrayOfMatrix &dpnd_dx) { - assert(Pe.ncols() == 4 and Pe.nrows() == 4); - - const Index nv = Pe.npages(); - const Index np = Pe.nbooks(); - const Index ne = Pe.nshelves(); - const Size nq = dpnd_dx.size(); - - Array aoaotm(np, muelmat_matrix(nq, nv)); - - for (Index ip = 0; ip < np; ip++) { - for (Size iq = 0; iq < nq; iq++) { - aoaotm[ip][iq] = muelmat{0.0}; - for (Index iv = 0; iv < nv; iv++) { - for (Index ie = 0; ie < ne; ie++) { - auto m = Pe[ie, ip, iv, joker, joker]; - aoaotm[ip][iq, iv] += dpnd_dx[iq][ie, ip] * muelmat(m[0, 0], - m[0, 1], - m[0, 2], - m[0, 3], - m[1, 0], - m[1, 1], - m[1, 2], - m[1, 3], - m[2, 0], - m[2, 1], - m[2, 2], - m[2, 3], - m[3, 0], - m[3, 1], - m[3, 2], - m[3, 3]); - } - } - } - } - return aoaotm; -} - -namespace { -void setBackscatterTransmission(stokvec_vector &out, - const stokvec_vector &I0, - const muelmat_vector &Tr, - const muelmat_vector &Tf, - const muelmat_vector &Z) { - for (Size i = 0; i < out.size(); i++) out[i] = Tr[i] * Z[i] * Tf[i] * I0[i]; -} - -void setBackscatterTransmissionDerivative(stokvec_matrix &out, - const stokvec_vector &I0, - const muelmat_vector &Tr, - const muelmat_vector &Tf, - const muelmat_matrix &dZ) { - for (Index j = 0; j < dZ.nrows(); j++) - for (Index i = 0; i < dZ.ncols(); i++) out[j, i] += Tr[i] * dZ[j, i] * Tf[i] * I0[i]; -} -} // namespace - -void bulk_backscatter_commutative_transmission_rte(Array &I, - Array> &dI, - const stokvec_vector &I_incoming, - const Array &T, - const Array &PiTf, - const Array &PiTr, - const Array &Z, - const Array &dT1, - const Array &dT2, - const Array &dZ) { - const Size np = dT1.size(); - const Index nv = np ? dT1.front().ncols() : 0; - const Index nq = np ? dT1.front().nrows() : 0; - - // For all transmission, the I-vector is the same - for (Size ip = 0; ip < np; ip++) setBackscatterTransmission(I[ip], I_incoming, PiTr[ip], PiTf[ip], Z[ip]); - - for (Size ip = 0; ip < np; ip++) { - setBackscatterTransmissionDerivative(dI[ip][ip], I_incoming, PiTr[ip], PiTf[ip], dZ[ip]); - } - - for (Size ip = 0; ip < np; ip++) { - for (Size j = ip; j < np; j++) { - for (Index iq = 0; iq < nq; iq++) { - for (Index iv = 0; iv < nv; iv++) { - dI[ip][j][iq, iv] += T[ip][iv] * ((dT1[ip][iq, iv] + dT2[ip][iq, iv]) * I[j][iv]); - - if (j < np - 1 and j > ip) - dI[ip][j][iq, iv] += inv(T[ip + 1][iv]) * ((dT1[ip][iq, iv] + dT2[ip][iq, iv]) * I[j][iv]); - } - } - } - } -} - namespace { Numeric cos_scat_angle(const Vector2 &los_in, const Vector2 &los_out) { using Conversion::cosd; diff --git a/src/core/rtepack/rtepack_scattering.h b/src/core/rtepack/rtepack_scattering.h index 96c0f9104e..4527d284bd 100644 --- a/src/core/rtepack/rtepack_scattering.h +++ b/src/core/rtepack/rtepack_scattering.h @@ -7,52 +7,5 @@ #include "rtepack_stokes_vector.h" namespace rtepack { -/** Polarized monostatic-radar return. - * - * The order is explicit because polarized outgoing and return transmissions - * need neither be equal nor commute. - */ -constexpr stokvec radar_return(const muelmat &return_transmission, - const muelmat &backscatter, - const muelmat &outgoing_transmission, - const stokvec &transmitted) { - return return_transmission * (backscatter * (outgoing_transmission * transmitted)); -} - -/** Directional derivative of radar_return. - * - * Implements the complete product rule for - * ``Tr * Z * Tf * transmitted``. Derivatives of the transmitted Stokes - * vector are deliberately outside this atmospheric-optics helper. - */ -constexpr stokvec radar_return_derivative(const muelmat &return_transmission, - const muelmat &return_transmission_derivative, - const muelmat &backscatter, - const muelmat &backscatter_derivative, - const muelmat &outgoing_transmission, - const muelmat &outgoing_transmission_derivative, - const stokvec &transmitted) { - const stokvec outgoing = outgoing_transmission * transmitted; - const stokvec outgoing_derivative = outgoing_transmission_derivative * transmitted; - const stokvec scattered = backscatter * outgoing; - const stokvec scattered_derivative = backscatter_derivative * outgoing + backscatter * outgoing_derivative; - return return_transmission_derivative * scattered + return_transmission * scattered_derivative; -} - -Array bulk_backscatter(const ConstTensor5View &Pe, const ConstMatrixView &pnd); - -Array bulk_backscatter_derivative(const ConstTensor5View &Pe, const ArrayOfMatrix &dpnd_dx); - -void bulk_backscatter_commutative_transmission_rte(Array &I, - Array> &dI, - const stokvec_vector &I_incoming, - const Array &T, - const Array &PiTf, - const Array &PiTr, - const Array &Z, - const Array &dT1, - const Array &dT2, - const Array &dZ); - muelmat rayleigh_scattering(const Vector2 &los_in, const Vector2 &los_out, const Numeric depolarization_factor); } // namespace rtepack diff --git a/src/core/rtepack/test/CMakeLists.txt b/src/core/rtepack/test/CMakeLists.txt index 2e4fcfa692..56ae926f12 100644 --- a/src/core/rtepack/test/CMakeLists.txt +++ b/src/core/rtepack/test/CMakeLists.txt @@ -16,3 +16,8 @@ add_executable(test_rtepack_polarization test_rtepack_polarization.cpp) target_link_libraries(test_rtepack_polarization PUBLIC rtepack) add_test(NAME "cpp.fast.test_rtepack_polarization" COMMAND test_rtepack_polarization) add_dependencies(check-deps test_rtepack_polarization) + +add_executable(test_rtepack_radar test_rtepack_radar.cpp) +target_link_libraries(test_rtepack_radar PUBLIC rtepack) +add_test(NAME "cpp.fast.test_rtepack_radar" COMMAND test_rtepack_radar) +add_dependencies(check-deps test_rtepack_radar) diff --git a/src/core/rtepack/test/test_rtepack_radar.cpp b/src/core/rtepack/test/test_rtepack_radar.cpp new file mode 100644 index 0000000000..f7a14cb145 --- /dev/null +++ b/src/core/rtepack/test/test_rtepack_radar.cpp @@ -0,0 +1,172 @@ +#include + +#include +#include +#include +#include + +namespace { +using namespace rtepack; + +constexpr Size np = 4; +constexpr Size nf = 1; +constexpr Size nq = 3; + +//! Polarized levels whose neighbouring propagation matrices do not commute +struct Levels { + std::vector K_out, K_ret; + std::vector dK_out, dK_ret; + std::vector Z; + std::vector dZ; + Vector r{300.0, 250.0, 400.0}; + + Levels() + : K_out(np, propmat_vector(nf)), + K_ret(np, propmat_vector(nf)), + dK_out(np, propmat_matrix(nq, nf)), + dK_ret(np, propmat_matrix(nq, nf)), + Z(np, muelmat_vector(nf)), + dZ(np, muelmat_matrix(nq, nf)) { + for (Size i = 0; i < np; i++) { + const Numeric s = static_cast(i + 1); + K_out[i][0] = propmat{1e-3 * s, 2e-4, -1e-4 * s, 5e-5, 1e-4, -2e-4 / s, 3e-5 * s}; + K_ret[i][0] = propmat{1.1e-3 * s, -1e-4 * s, 2e-4, -5e-5, 2e-4 / s, 1e-4, -3e-5}; + Z[i][0] = + 1e-5 * + muelmat{1.0 + s, 0.2, 0.1 * s, 0.0, 0.2, 0.9, -0.1, 0.05 * s, -0.1, 0.1, 0.8 * s, 0.2, 0.0, -0.05, 0.1, 0.7}; + + // The three derivative directions perturb K_out, K_ret and Z in turn + dK_out[i][0, 0] = propmat{1e-4, -3e-5 * s, 2e-5, 1e-5, -2e-5, 4e-5, 1e-5}; + dK_ret[i][1, 0] = propmat{-2e-4, 1e-5, -4e-5 * s, 2e-5, 3e-5, -1e-5, 2e-5}; + dZ[i][2, 0] = 1e-6 * muelmat{1.0, 0.3, 0.0, -0.2, 0.1, 1.0, 0.2, 0.0, 0.0, -0.1, 1.0, 0.1, 0.2, 0.0, 0.1, s}; + } + } + + void perturb(Size level, Size q, Numeric h) { + switch (q) { + case 0: K_out[level][0] += h * dK_out[level][0, 0]; break; + case 1: K_ret[level][0] += h * dK_ret[level][1, 0]; break; + default: Z[level][0] += h * dZ[level][2, 0]; break; + } + } + + void solve(stokvec_matrix& I, + stokvec_tensor4& dI, + const std::vector& samples, + Size ngate, + TransmittanceOption option) const { + const stokvec transmitted{1.0, 0.3, -0.2, 0.1}; + I.resize(ngate, nf); + dI.resize(ngate, np, nq, nf); + radar_single_scattering(I, dI, transmitted, K_out, dK_out, K_ret, dK_ret, Z, dZ, r, samples, option); + } +}; + +const char* name(TransmittanceOption option) { return option == TransmittanceOption::magop ? "magop" : "constant"; } + +bool test_derivatives(TransmittanceOption option) { + const std::vector samples{{.layer = 1, .x = 0.3, .weight = 0.4, .gate = 0}, + {.layer = 2, .x = 0.5, .weight = 0.6, .gate = 0}, + {.layer = 2, .x = 0.9, .weight = 0.2, .gate = 1}, + {.layer = 3, .x = 0.25, .weight = 0.8, .gate = 1}}; + constexpr Size ngate = 2; + + const Levels levels; + stokvec_matrix I; + stokvec_tensor4 dI; + levels.solve(I, dI, samples, ngate, option); + + // The derivatives are linear in the direction, so central differences are + // accurate to the second order in h. + constexpr Numeric h = 1e-3; + stokvec_matrix Ip, Im; + stokvec_tensor4 dummy; + bool ok = true; + for (Size level = 0; level < np; level++) { + for (Size q = 0; q < nq; q++) { + Levels plus = levels, minus = levels; + plus.perturb(level, q, h); + minus.perturb(level, q, -h); + plus.solve(Ip, dummy, samples, ngate, option); + minus.solve(Im, dummy, samples, ngate, option); + + for (Size gate = 0; gate < ngate; gate++) { + for (Size s = 0; s < 4; s++) { + const Numeric numeric = (Ip[gate, 0][s] - Im[gate, 0][s]) / (2.0 * h); + const Numeric analytic = dI[gate, level, q, 0][s]; + const Numeric scale = std::max(std::abs(I[gate, 0][s]), 1e-12); + if (not std::isfinite(analytic) or std::abs(analytic - numeric) > 1e-7 * scale) { + std::cerr << "Radar derivative for " << name(option) << ", level " << level << ", direction " << q + << ", gate " << gate << ", Stokes " << s << ": analytic " << analytic << ", numeric " << numeric + << '\n'; + ok = false; + } + } + } + } + } + return ok; +} + +bool test_layer_continuity(TransmittanceOption option) { + // The far end of a layer and the near end of the next layer are the same + // point, so the partial-layer transmissions must agree with the full ones. + const std::vector samples{{.layer = 1, .x = 1.0, .weight = 1.0, .gate = 0}, + {.layer = 2, .x = 0.0, .weight = 1.0, .gate = 1}, + {.layer = 2, .x = 1.0, .weight = 1.0, .gate = 2}, + {.layer = 3, .x = 0.0, .weight = 1.0, .gate = 3}}; + + const Levels levels; + stokvec_matrix I; + stokvec_tensor4 dI; + levels.solve(I, dI, samples, 4, option); + + bool ok = true; + for (Size gate = 0; gate < 4; gate += 2) { + for (Size s = 0; s < 4; s++) { + const Numeric far = I[gate, 0][s]; + const Numeric near = I[gate + 1, 0][s]; + if (std::abs(far - near) > 1e-12 * std::max(std::abs(far), 1e-12)) { + std::cerr << "Radar return for " << name(option) << " is discontinuous at level " << gate / 2 + 1 << ", Stokes " + << s << ": " << far << " vs " << near << '\n'; + ok = false; + } + } + } + + for (Size gate = 0; gate < 4; gate++) + for (Size level = 0; level < np; level++) + for (Size q = 0; q < nq; q++) + for (Size s = 0; s < 4; s++) + if (not std::isfinite(dI[gate, level, q, 0][s])) { + std::cerr << "Radar derivative for " << name(option) << " is not finite at a level\n"; + return false; + } + return ok; +} + +bool test_magnus_is_used() { + const std::vector samples{{.layer = 3, .x = 0.5, .weight = 1.0, .gate = 0}}; + + const Levels levels; + stokvec_matrix constant, magnus; + stokvec_tensor4 dI; + levels.solve(constant, dI, samples, 1, TransmittanceOption::constant); + levels.solve(magnus, dI, samples, 1, TransmittanceOption::magop); + + if (constant[0, 0] == magnus[0, 0]) { + std::cerr << "The Magnus operator does not change the return of non-commuting layers\n"; + return false; + } + return true; +} +} // namespace + +int main() { + bool ok = test_magnus_is_used(); + for (const auto option : {TransmittanceOption::constant, TransmittanceOption::magop}) { + ok = test_derivatives(option) and ok; + ok = test_layer_continuity(option) and ok; + } + return ok ? 0 : 1; +} diff --git a/src/m_montecarlo_radar.cc b/src/m_montecarlo_radar.cc index a6c9bbc59e..10234f2bbb 100644 --- a/src/m_montecarlo_radar.cc +++ b/src/m_montecarlo_radar.cc @@ -1,7 +1,7 @@ #include -#include #include #include +#include #include #include #include @@ -9,109 +9,19 @@ #include #include #include - -namespace { -Numeric normalized_delta_aa(Numeric aa) { - while (aa < -180.0) aa += 360.0; - while (aa > 180.0) aa -= 360.0; - return aa; -} - -Numeric ze_factor(const Numeric frequency, const Numeric k2) { - const Numeric wavelength = Constant::speed_of_light / frequency; - return 4e18 * std::pow(wavelength, 4) / (std::pow(Constant::pi, 4) * k2); -} - -struct OpticalPoint { - Propmat outgoing_extinction{}; - Propmat return_extinction{}; - Muelmat backscatter{0.0}; -}; - -OpticalPoint optical_point(const Workspace& ws, - const Numeric frequency, - const PropagationPathPoint& ray_point, - const AtmPoint& atm_point, - const ArrayOfScatteringSpecies& scattering_species, - const Agenda& spectral_propmat_agenda) { - const AscendingGrid freq_grid{frequency}; - - PropmatVector gas_return; - StokvecVector src_return; - PropmatMatrix gas_return_jac; - StokvecMatrix src_return_jac; - spectral_propmat_agendaExecute(ws, - gas_return, - src_return, - gas_return_jac, - src_return_jac, - freq_grid, - Vector3{0.0, 0.0, 0.0}, - JacobianTargets{}, - {}, - ray_point, - atm_point, - spectral_propmat_agenda); - - PropagationPathPoint outgoing_point = ray_point; - outgoing_point.los = path::mirror(ray_point.los); - PropmatVector gas_outgoing; - StokvecVector src_outgoing; - PropmatMatrix gas_outgoing_jac; - StokvecMatrix src_outgoing_jac; - spectral_propmat_agendaExecute(ws, - gas_outgoing, - src_outgoing, - gas_outgoing_jac, - src_outgoing_jac, - freq_grid, - Vector3{0.0, 0.0, 0.0}, - JacobianTargets{}, - {}, - outgoing_point, - atm_point, - spectral_propmat_agenda); - - const Vector2 transmit_los = outgoing_point.los; - const Vector2 receive_los = ray_point.los; - const auto at_pole = [](const Numeric za) { return std::abs(std::sin(za * Constant::pi / 180.0)) < 1e-12; }; - const Vector delta_aa{at_pole(transmit_los[0]) or at_pole(receive_los[0]) - ? 0.0 - : normalized_delta_aa(receive_los[1] - transmit_los[1])}; - auto receive_za = - std::make_shared(scattering::IrregularZenithAngleGrid(Vector{receive_los[0]})); - const auto bulk_outgoing = scattering_species.get_bulk_scattering_properties_aro_gridded( - atm_point, freq_grid, Vector{transmit_los[0]}, delta_aa, receive_za); - auto transmit_za = - std::make_shared(scattering::IrregularZenithAngleGrid(Vector{transmit_los[0]})); - const auto bulk_return = scattering_species.get_bulk_scattering_properties_aro_gridded( - atm_point, freq_grid, Vector{receive_los[0]}, Vector{0.0}, transmit_za); - - if (not bulk_outgoing.phase_matrix.has_value()) - arts::user_error("MCRadar requires phase matrices from all scattering species"); - - const auto& phase = *bulk_outgoing.phase_matrix; - - OpticalPoint out; - out.outgoing_extinction = gas_outgoing[0]; - out.outgoing_extinction += bulk_outgoing.extinction_matrix[0, 0, 0].expand(); - out.return_extinction = gas_return[0]; - out.return_extinction += bulk_return.extinction_matrix[0, 0, 0].expand(); - out.backscatter = phase[0, 0, 0, 0, 0]; - return out; -} - -} // namespace +#include void MCRadar(const Workspace& ws, StokvecVector& radar_signal, StokvecVector& radar_error, + Vector& radar_time, const AtmField& atm_field, const SurfaceField& surf_field, const ArrayOfScatteringSpecies& scattering_species, const MCAntenna& mc_antenna, const Agenda& ray_path_observer_agenda, const Agenda& spectral_propmat_agenda, + const TransmittanceOption& rte_option, const Numeric& frequency, const Vector3& sensor_pos, const Vector2& sensor_los, @@ -121,25 +31,63 @@ void MCRadar(const Workspace& ws, const Index& mc_max_iter, const Index& mc_max_scatorder, const Numeric& k2, - const String& unit) try { + const RadarUnit& unit, + const RadarRangeMode& range_mode, + const Numeric& ze_tref, + const Numeric& pext_scaling) try { ARTS_TIME_REPORT + const RadarRangeMode mode = radar::resolve(range_mode, range_bins.empty() ? 0.0 : range_bins.back()); + if (frequency <= 0.0) arts::user_error("frequency must be positive"); if (range_bins.size() < 2) arts::user_error("range_bins must contain at least two edges"); - if (range_bins.front() < 0.0) arts::user_error("range_bins cannot contain negative distances"); + if (mode != RadarRangeMode::Altitude and range_bins.front() < 0.0) + arts::user_error("range_bins cannot contain negative distances or times"); if (mc_max_iter <= 0) arts::user_error("mc_max_iter must be positive"); if (mc_max_scatorder != 1) arts::user_error("This first ARTS3 MCRadar implementation supports mc_max_scatorder=1; got {}", mc_max_scatorder); if (scattering_species.species.empty()) arts::user_error("MCRadar requires at least one scattering species"); - if (unit != "1" and unit != "Ze") arts::user_error(R"(unit must be either "1" or "Ze")"); - if (k2 <= 0.0) arts::user_error("k2 must be positive"); - - const Index nbins = static_cast(range_bins.size()) - 1; - radar_signal.resize(nbins); - radar_error.resize(nbins); - std::fill(radar_signal.begin(), radar_signal.end(), Stokvec{0.0}); - std::fill(radar_error.begin(), radar_error.end(), Stokvec{0.0}); - StokvecVector squared_sum(nbins, Stokvec{0.0}); + if (unit == RadarUnit::dBZe) + arts::user_error("MCRadar returns a mean and its standard error, which are linear quantities; use Ze, not dBZe"); + if (pext_scaling < 0.0 or pext_scaling > 2.0) + arts::user_error("pext_scaling must be between 0 and 2; got {}", pext_scaling); + + const Size nbins = range_bins.size() - 1; + const AscendingGrid freq_grid{frequency}; + const JacobianTargets jac_targets{}; + const Numeric factor = radar::unit_factor(unit, frequency, ze_tref, k2); + + const radar::GasModel gas = [&](PropmatVector& propmat, + PropmatMatrix& propmat_jac, + const PropagationPathPoint& point, + const AtmPoint& atm_point) { + StokvecVector srcvec; + StokvecMatrix srcvec_jac; + spectral_propmat_agendaExecute(ws, + propmat, + srcvec, + propmat_jac, + srcvec_jac, + freq_grid, + Vector3{0.0, 0.0, 0.0}, + jac_targets, + {}, + point, + atm_point, + spectral_propmat_agenda); + }; + + // The mean is weighted by the antenna gain of every sampled direction, so a + // uniformly filled beam returns the same value as a pencil beam. Weighted + // sums of the return, and of its square, give the standard error. + Numeric weight_sum = 0.0; + Numeric squared_weight_sum = 0.0; + StokvecVector sum(nbins, Stokvec{}); + StokvecVector squared_weight_signal_sum(nbins, Stokvec{}); + StokvecVector squared_sum(nbins, Stokvec{}); + Vector time_sum(nbins, 0.0); + Vector coverage_sum(nbins, 0.0); + std::vector covered(nbins, false); RandomNumberGenerator<> rng(mc_seed); // The general RNG avoids reusing seeds process-wide. A forward model with @@ -150,84 +98,73 @@ void MCRadar(const Workspace& ws, for (Index i = 0; i < 3; ++i) for (Index j = 0; j < 3; ++j) enu_to_ant[i, j] = ant_to_enu[j, i]; + StokvecMatrix signal; + StokvecTensor4 signal_jac; for (Index iter = 0; iter < mc_max_iter; ++iter) { - StokvecVector photon(nbins, Stokvec{0.0}); - const auto [sampled_los, ray_rotation] = mc_antenna.draw_los(rng, ant_to_enu, sensor_los); - const Numeric antenna_weight = mc_antenna.return_los(ray_rotation, enu_to_ant); - const Muelmat tx_rotation = rotmat_stokes(1.0, 1.0, ant_to_enu, ray_rotation); - const Muelmat rx_rotation = rotmat_stokes(-1.0, 1.0, ray_rotation, ant_to_enu); - const Stokvec transmitted = tx_rotation * mc_y_tx; + const auto [sampled_los, ray_rotation] = mc_antenna.draw_los(rng, ant_to_enu, sensor_los); + const Numeric antenna_weight = mc_antenna.return_los(ray_rotation, enu_to_ant); + const Muelmat tx_rotation = rotmat_stokes(1.0, 1.0, ant_to_enu, ray_rotation); + const Muelmat rx_rotation = rotmat_stokes(-1.0, 1.0, ray_rotation, ant_to_enu); + weight_sum += antenna_weight; + squared_weight_sum += antenna_weight * antenna_weight; ArrayOfPropagationPathPoint ray_path; ray_path_observer_agendaExecute(ws, ray_path, sensor_pos, sampled_los, ray_path_observer_agenda); - if (ray_path.size() < 2) continue; - - const ArrayOfAtmPoint atm_path = forward_atm_path(ray_path, atm_field); - const Vector distances = path::distance(ray_path, surf_field.ellipsoid); - Muelmat outbound = Muelmat::id(); - Muelmat returning = Muelmat::id(); - Numeric travelled = 0.0; - OpticalPoint previous = - optical_point(ws, frequency, ray_path.front(), atm_path.front(), scattering_species, spectral_propmat_agenda); - - for (Size ip = 1; ip < ray_path.size(); ++ip) { - const Numeric ds = distances[ip - 1]; - if (ds <= 0.0) continue; - - const OpticalPoint current = - optical_point(ws, frequency, ray_path[ip], atm_path[ip], scattering_species, spectral_propmat_agenda); - const Muelmat outgoing_segment = rtepack::tran(previous.outgoing_extinction, current.outgoing_extinction, ds)(); - const Muelmat return_segment = rtepack::tran(current.return_extinction, previous.return_extinction, ds)(); - // Integrate only the overlap with each gate. Interpolate optical - // properties to its midpoint, including the extinction along the - // partial segment between the previous path point and that midpoint. - const auto first_edge = stdr::upper_bound(range_bins, travelled); - const Index first_bin = std::max(0, first_edge - range_bins.begin() - 1); - for (Index ibin = first_bin; ibin < nbins and range_bins[ibin] < travelled + ds; ++ibin) { - const Numeric begin = std::max(travelled, range_bins[ibin]); - const Numeric end = std::min(travelled + ds, range_bins[ibin + 1]); - if (end <= begin) continue; - const Numeric midpoint = 0.5 * (begin + end) - travelled; - const Numeric fraction = midpoint / ds; - const Propmat outgoing_extinction = - (1.0 - fraction) * previous.outgoing_extinction + fraction * current.outgoing_extinction; - const Propmat return_extinction = - (1.0 - fraction) * previous.return_extinction + fraction * current.return_extinction; - const Muelmat outgoing_midpoint = - rtepack::tran(previous.outgoing_extinction, outgoing_extinction, midpoint)() * outbound; - const Muelmat return_midpoint = - returning * rtepack::tran(return_extinction, previous.return_extinction, midpoint)(); - const Muelmat backscatter = (1.0 - fraction) * previous.backscatter + fraction * current.backscatter; - const Stokvec contribution = - antenna_weight * (end - begin) * - (rx_rotation * rtepack::radar_return(return_midpoint, backscatter, outgoing_midpoint, transmitted)); - photon[ibin] += contribution; - } - travelled += ds; - if (travelled >= range_bins.back()) break; - outbound = outgoing_segment * outbound; - returning = returning * return_segment; - previous = current; - } - - for (Index b = 0; b < nbins; ++b) { - const Numeric width = range_bins[b + 1] - range_bins[b]; - for (Index s = 0; s < 4; ++s) { - const Numeric value = photon[b][s] / width; - radar_signal[b][s] += value; - squared_sum[b][s] += value * value; + const radar::PathRange range = radar::path_range(ray_path, sensor_pos, surf_field.ellipsoid); + const Vector coordinate = radar::range_coordinate(range, ray_path, mode); + std::vector samples; + for (Size ib = 0; ib < nbins; ++ib) + radar::add_gate_samples(samples, coordinate, range_bins[ib], range_bins[ib + 1], ib); + const Size np = radar::sort_samples(samples); + if (samples.empty()) continue; + + const radar::Optics optics = + radar::optics(ray_path, np, atm_field, freq_grid, scattering_species, jac_targets, pext_scaling, gas); + radar::single_scattering( + signal, signal_jac, tx_rotation * mc_y_tx, optics, range.layer, samples, nbins, rte_option); + + const Vector coverage = radar::gate_coverage(samples, nbins); + const Vector time = radar::gate_sum(samples, range.time, nbins); + for (Size ib = 0; ib < nbins; ++ib) { + if (coverage[ib] == 0.0) continue; + + covered[ib] = true; + time_sum[ib] += antenna_weight * time[ib]; + coverage_sum[ib] += antenna_weight * coverage[ib]; + + const Stokvec value = rx_rotation * signal[ib, 0]; + for (Size s = 0; s < 4; ++s) { + sum[ib][s] += antenna_weight * value[s]; + squared_weight_signal_sum[ib][s] += antenna_weight * antenna_weight * value[s]; + squared_sum[ib][s] += antenna_weight * antenna_weight * value[s] * value[s]; } } } - const Numeric factor = unit == "Ze" ? ze_factor(frequency, k2) / (2.0 * Constant::pi) : 1.0; - for (Index b = 0; b < nbins; ++b) { - for (Index s = 0; s < 4; ++s) { - const Numeric mean = radar_signal[b][s] / static_cast(mc_max_iter); - const Numeric variance = std::max(0.0, squared_sum[b][s] / static_cast(mc_max_iter) - mean * mean); - radar_signal[b][s] = factor * mean; - radar_error[b][s] = factor * std::sqrt(variance / static_cast(mc_max_iter)); + if (not(weight_sum > 0.0)) arts::user_error("No sampled antenna direction has a positive antenna weight"); + + constexpr Numeric nan = std::numeric_limits::quiet_NaN(); + radar_signal.resize(nbins); + radar_error.resize(nbins); + radar_time.resize(nbins); + for (Size ib = 0; ib < nbins; ++ib) { + if (not covered[ib]) { + radar_signal[ib] = Stokvec{nan, nan, nan, nan}; + radar_error[ib] = Stokvec{nan, nan, nan, nan}; + radar_time[ib] = nan; + continue; + } + + radar_time[ib] = coverage_sum[ib] > 0.0 ? time_sum[ib] / coverage_sum[ib] : nan; + for (Size s = 0; s < 4; ++s) { + const Numeric mean = sum[ib][s] / weight_sum; + const Numeric variance = std::max(0.0, + squared_sum[ib][s] - 2.0 * mean * squared_weight_signal_sum[ib][s] + + mean * mean * squared_weight_sum) / + (weight_sum * weight_sum); + radar_signal[ib][s] = factor * mean; + radar_error[ib][s] = factor * std::sqrt(variance); } } } diff --git a/src/m_radar.cc b/src/m_radar.cc index a16e498929..fc18ecdd20 100644 --- a/src/m_radar.cc +++ b/src/m_radar.cc @@ -1,7 +1,7 @@ #include -#include #include #include +#include #include #include #include @@ -13,503 +13,262 @@ #include #include #include -#include -#include #include namespace { - -Numeric liebe93_k2(const Numeric frequency, const Numeric temperature) { - if (temperature < Constant::temperature_at_0c - 40.0 or temperature > Constant::temperature_at_0c + 100.0) - arts::user_error("Automatic radar k2 requires {} K <= ze_tref <= {} K; got {} K", - Constant::temperature_at_0c - 40.0, - Constant::temperature_at_0c + 100.0, - temperature); - if (frequency < 10e9 or frequency > 1000e9) - arts::user_error("Automatic radar k2 uses Liebe-93 and requires 10 GHz <= frequency <= 1000 GHz; got {} Hz", - frequency); - - const Numeric theta = 1.0 - 300.0 / temperature; - const Numeric e0 = 77.66 - 103.3 * theta; - const Numeric e1 = 0.0671 * e0; - const Numeric f1 = 20.2 + 146.0 * theta + 316.0 * theta * theta; - const Numeric e2 = 3.52; - const Numeric f2 = 39.8 * f1; - const Complex ifGHz{0.0, frequency / 1e9}; - const Complex n = std::sqrt(e2 + (e1 - e2) / (1.0 - ifGHz / f2) + (e0 - e1) / (1.0 - ifGHz / f1)); - const Complex n2 = n * n; - return std::norm((n2 - 1.0) / (n2 + 2.0)); -} - -Numeric ze_factor(const Numeric frequency, const Numeric ze_tref, const Numeric k2) { - const Numeric dielectric = k2 < 0.0 ? liebe93_k2(frequency, ze_tref) : k2; - if (dielectric <= 0.0) arts::user_error("Radar k2 must be positive, or negative to request Liebe-93; got {}", k2); - const Numeric wavelength = Constant::speed_of_light / frequency; - return 4e18 * std::pow(wavelength, 4) / (std::pow(Constant::pi, 4) * dielectric); -} - -struct OpticalProfile { - std::vector outgoing_extinction; - std::vector return_extinction; - std::vector gas_extinction; - std::vector particle_extinction; - std::vector backscatter; - std::vector outgoing_extinction_jac; - std::vector return_extinction_jac; - std::vector backscatter_jac; +struct Settings { + Stokvec transmitted; + RadarRangeMode range_mode; + RadarUnit unit; + Numeric ze_tref; + Numeric k2; + Numeric dbze_min; + Numeric pext_scaling; + TransmittanceOption rte_option; }; -Numeric normalized_delta_aa(Numeric aa) { - while (aa < -180.0) aa += 360.0; - while (aa > 180.0) aa -= 360.0; - return aa; -} - -OpticalProfile optical_profile(const Workspace& ws, - const AscendingGrid& frequency, - const ArrayOfPropagationPathPoint& ray_path, - const AtmField& atm_field, - const ArrayOfScatteringSpecies& scattering_species, - const Agenda& spectral_propmat_agenda, - const JacobianTargets& jac_targets, - const Numeric pext_scaling) { - const Size np = ray_path.size(); - const Size nf = frequency.size(); - const Size nq = jac_targets.target_count(); - - OpticalProfile out; - out.outgoing_extinction.resize(np, PropmatVector(nf)); - out.return_extinction.resize(np, PropmatVector(nf)); - out.gas_extinction.resize(np, PropmatVector(nf)); - out.particle_extinction.resize(np, PropmatVector(nf)); - out.backscatter.resize(np, MuelmatVector(nf, Muelmat{0.0})); - out.outgoing_extinction_jac.resize(np, PropmatMatrix(nq, nf)); - out.return_extinction_jac.resize(np, PropmatMatrix(nq, nf)); - out.backscatter_jac.resize(np, MuelmatMatrix(nq, nf, Muelmat{0.0})); - - const ArrayOfAtmPoint atm_path = forward_atm_path(ray_path, atm_field); - - for (Size ip = 0; ip < np; ++ip) { - PropmatVector gas_return; - StokvecVector source_return; - PropmatMatrix gas_return_jac; - StokvecMatrix source_return_jac; - spectral_propmat_agendaExecute(ws, - gas_return, - source_return, - gas_return_jac, - source_return_jac, - frequency, - Vector3{0.0, 0.0, 0.0}, - jac_targets, - {}, - ray_path[ip], - atm_path[ip], - spectral_propmat_agenda); - PropagationPathPoint outgoing_point = ray_path[ip]; - outgoing_point.los = path::mirror(outgoing_point.los); - PropmatVector gas_outgoing; - StokvecVector source_outgoing; - PropmatMatrix gas_outgoing_jac; - StokvecMatrix source_outgoing_jac; - spectral_propmat_agendaExecute(ws, - gas_outgoing, - source_outgoing, - gas_outgoing_jac, - source_outgoing_jac, - frequency, - Vector3{0.0, 0.0, 0.0}, - jac_targets, - {}, - outgoing_point, - atm_path[ip], - spectral_propmat_agenda); - - // Propagation-path LOS is the propagation direction toward the receiver; - // its mirror is the transmitted propagation direction. ARO phase - // matrices take the latter as incident and the former as scattered. - const Vector2 receive_los = ray_path[ip].los; - const Vector2 transmit_los = outgoing_point.los; - // Azimuth is undefined at zenith and nadir. Use its canonical value at - // either pole to avoid a 0*NaN rotation without changing the scattering - // angle or the physical polarization basis. - const auto at_pole = [](const Numeric za) { return std::abs(std::sin(za * Constant::pi / 180.0)) < 1e-12; }; - const Vector delta_aa{at_pole(transmit_los[0]) or at_pole(receive_los[0]) - ? 0.0 - : normalized_delta_aa(receive_los[1] - transmit_los[1])}; - auto receive_za = - std::make_shared(scattering::IrregularZenithAngleGrid(Vector{receive_los[0]})); - const auto bulk_outgoing = scattering_species.get_bulk_scattering_properties_aro_gridded( - atm_path[ip], frequency, Vector{transmit_los[0]}, delta_aa, receive_za); - auto transmit_za = - std::make_shared(scattering::IrregularZenithAngleGrid(Vector{transmit_los[0]})); - const auto bulk_return = scattering_species.get_bulk_scattering_properties_aro_gridded( - atm_path[ip], frequency, Vector{receive_los[0]}, Vector{0.0}, transmit_za); - if (not bulk_outgoing.phase_matrix.has_value()) - arts::user_error("Radar single scattering requires phase matrices from all scattering species"); - - const auto& phase = *bulk_outgoing.phase_matrix; - for (Size iv = 0; iv < nf; ++iv) { - out.gas_extinction[ip][iv] = gas_outgoing[iv]; - out.particle_extinction[ip][iv] = bulk_outgoing.extinction_matrix[0, iv, 0].expand(); - out.outgoing_extinction[ip][iv] = gas_outgoing[iv]; - out.outgoing_extinction[ip][iv] += pext_scaling * out.particle_extinction[ip][iv]; - out.return_extinction[ip][iv] = gas_return[iv]; - out.return_extinction[ip][iv] += pext_scaling * bulk_return.extinction_matrix[0, iv, 0].expand(); - out.backscatter[ip][iv] = phase[0, iv, 0, 0, 0]; - } - for (const auto& target : jac_targets.atm) { - auto receive_za_jac = - std::make_shared(scattering::IrregularZenithAngleGrid(Vector{receive_los[0]})); - const auto dbulk_outgoing = scattering_species.get_bulk_scattering_properties_aro_gridded_derivative( - atm_path[ip], frequency, Vector{transmit_los[0]}, delta_aa, receive_za_jac, target.type); - auto transmit_za_jac = - std::make_shared(scattering::IrregularZenithAngleGrid(Vector{transmit_los[0]})); - const auto dbulk_return = scattering_species.get_bulk_scattering_properties_aro_gridded_derivative( - atm_path[ip], frequency, Vector{receive_los[0]}, Vector{0.0}, transmit_za_jac, target.type); - if (not dbulk_outgoing.phase_matrix) arts::user_error("Radar derivatives require scattering phase matrices"); - const auto& dphase = *dbulk_outgoing.phase_matrix; - for (Size iv = 0; iv < nf; ++iv) { - out.outgoing_extinction_jac[ip][target.target_pos, iv] = gas_outgoing_jac[target.target_pos, iv]; - out.outgoing_extinction_jac[ip][target.target_pos, iv] += - pext_scaling * dbulk_outgoing.extinction_matrix[0, iv, 0].expand(); - out.return_extinction_jac[ip][target.target_pos, iv] = gas_return_jac[target.target_pos, iv]; - out.return_extinction_jac[ip][target.target_pos, iv] += - pext_scaling * dbulk_return.extinction_matrix[0, iv, 0].expand(); - out.backscatter_jac[ip][target.target_pos, iv] = dphase[0, iv, 0, 0, 0]; - } - } +void check_settings(const ArrayOfSensorObsel& sensor, + const Matrix& limits, + const JacobianTargets& jac_targets, + const ArrayOfScatteringSpecies& scattering_species, + const Settings& settings) { + const Size ny = sensor.size(); + if (ny == 0) arts::user_error("measurement_sensor cannot be empty for a radar calculation"); + if (limits.nrows() != static_cast(ny) or limits.ncols() != 2) + arts::user_error("radar_range_limits must have shape [{}, 2], got {:B,}", ny, limits.shape()); + for (Size iy = 0; iy < ny; ++iy) { + sensor[iy].check(); + if (not(limits[iy, 0] < limits[iy, 1])) + arts::user_error( + "Radar range limits in row {} are not strictly increasing: [{}, {}]", iy, limits[iy, 0], limits[iy, 1]); } - - return out; + if (settings.transmitted.I() != 1.0) + arts::user_error("The transmitted Stokes I component must equal one, got {}", settings.transmitted.I()); + if (settings.pext_scaling < 0.0 or settings.pext_scaling > 2.0) + arts::user_error("pext_scaling must be between 0 and 2; got {}", settings.pext_scaling); + if (scattering_species.species.empty()) arts::user_error("Radar single scattering requires scattering species"); + if (not jac_targets.surf.empty() or not jac_targets.subsurf.empty() or not jac_targets.line.empty() or + not jac_targets.sensor.empty() or not jac_targets.error.empty()) + arts::user_error("The radar forward model currently supports atmospheric Jacobian targets only"); } -enum class RadarRangeMode : char { Altitude, Distance, RoundTripTime }; - -RadarRangeMode range_mode(const String& mode, const Matrix& limits) { - if (mode == "Altitude") return RadarRangeMode::Altitude; - if (mode == "Distance") return RadarRangeMode::Distance; - if (mode == "RoundTripTime") return RadarRangeMode::RoundTripTime; - if (mode == "Legacy") { - Numeric largest = -std::numeric_limits::infinity(); - for (Index i = 0; i < limits.nrows(); ++i) - for (Index j = 0; j < limits.ncols(); ++j) largest = std::max(largest, limits[i, j]); - return largest > 1.0 ? RadarRangeMode::Altitude : RadarRangeMode::RoundTripTime; - } - arts::user_error("Unknown radar range_mode \"{}\". Expected Altitude, Distance, RoundTripTime, or Legacy.", mode); +RadarRangeMode resolved_range_mode(const RadarRangeMode mode, const Matrix& limits) { + Numeric largest = -std::numeric_limits::infinity(); + for (Index i = 0; i < limits.nrows(); ++i) + for (Index j = 0; j < limits.ncols(); ++j) largest = std::max(largest, limits[i, j]); + return radar::resolve(mode, largest); } -Vector range_coordinate(const ArrayOfPropagationPathPoint& ray_path, - const Vector3& observer, - const Vector2& ellipsoid, - const RadarRangeMode mode) { - const Size np = ray_path.size(); - Vector out(np); - if (mode == RadarRangeMode::Altitude) { - for (Size ip = 0; ip < np; ++ip) out[ip] = ray_path[ip].pos[0]; - return out; - } - - const Numeric initial = path::distance(observer, ray_path.front().pos, ellipsoid); - if (mode == RadarRangeMode::Distance) { - out[0] = initial; - const Vector segment = path::distance(ray_path, ellipsoid); - for (Size ip = 1; ip < np; ++ip) out[ip] = out[ip - 1] + segment[ip - 1]; - return out; - } - - out[0] = 2.0 * initial / Constant::speed_of_light; - const Vector segment = path::distance(ray_path, ellipsoid); - for (Size ip = 1; ip < np; ++ip) { - out[ip] = - out[ip - 1] + segment[ip - 1] * (ray_path[ip - 1].ngroup + ray_path[ip].ngroup) / Constant::speed_of_light; - } - return out; -} +//! A sensor geometry, its path, and the samples of the gates of all sensor elements that use it +struct PreparedSimulation { + const AscendingGrid* frequency{}; + Size iposlos{}; + ArrayOfPropagationPathPoint ray_path{}; + radar::PathRange range{}; + std::vector gates{}; + std::vector samples{}; + Size np{}; +}; -template std::optional bin_average(const Vector& x, const std::vector& y, Numeric lo, Numeric hi) { - if (x.size() != y.size()) arts::user_error("Internal radar range/profile size mismatch"); - if (x.size() < 2) return std::nullopt; - - const bool increasing = x.back() > x.front(); - const Numeric xmin = std::min(x.front(), x.back()); - const Numeric xmax = std::max(x.front(), x.back()); - lo = std::max(lo, xmin); - hi = std::min(hi, xmax); - if (not(lo < hi)) return std::nullopt; - - auto interpolate = [&](const Numeric q) { - Size i = 0; - if (increasing) { - while (i + 1 < x.size() and x[i + 1] < q) ++i; - } else { - while (i + 1 < x.size() and x[i + 1] > q) ++i; +std::vector prepare_simulations(const Workspace& ws, + const ArrayOfSensorObsel& sensor, + const Matrix& limits, + const RadarRangeMode range_mode, + const SurfaceField& surf_field, + const Agenda& ray_path_observer_agenda) { + std::vector out; + for (const auto& [frequency, poslos, iposlos] : collect_simulations(sensor)) { + PreparedSimulation simulation{.frequency = &frequency, .iposlos = iposlos}; + ray_path_observer_agendaExecute( + ws, simulation.ray_path, poslos[iposlos].pos, poslos[iposlos].los, ray_path_observer_agenda); + simulation.range = radar::path_range(simulation.ray_path, poslos[iposlos].pos, surf_field.ellipsoid); + + const Vector coordinate = radar::range_coordinate(simulation.range, simulation.ray_path, range_mode); + for (Size iy = 0; iy < sensor.size(); ++iy) { + if (not sensor[iy].same_freqs(frequency) or not sensor[iy].same_poslos(poslos)) continue; + radar::add_gate_samples(simulation.samples, coordinate, limits[iy, 0], limits[iy, 1], simulation.gates.size()); + simulation.gates.push_back(iy); } - if (i + 1 >= x.size()) return y.back(); - const Numeric w = (q - x[i]) / (x[i + 1] - x[i]); - return (1.0 - w) * y[i] + w * y[i + 1]; - }; + simulation.np = radar::sort_samples(simulation.samples); - std::vector q{lo}; - q.reserve(x.size() + 2); - for (auto v : x) - if (v > lo and v < hi) q.push_back(v); - q.push_back(hi); - std::ranges::sort(q); - - T integral{0.0}; - T past = interpolate(q.front()); - for (Size i = 1; i < q.size(); ++i) { - const T now = interpolate(q[i]); - integral += 0.5 * (q[i] - q[i - 1]) * (past + now); - past = now; + out.push_back(std::move(simulation)); } - return integral * (1.0 / (hi - lo)); + return out; } -struct PreparedSimulation { - const AscendingGrid* frequency{}; - const SensorPosLosVector* poslos{}; - Size iposlos{}; - ArrayOfPropagationPathPoint ray_path{}; -}; - struct RadarResult { Vector measurement; Matrix jacobian; Matrix auxiliary; + Vector time; }; RadarResult radar_forward(const Workspace& ws, const ArrayOfSensorObsel& sensor, - const Matrix& limits, const std::vector& simulations, const AtmField& atm_field, - const SurfaceField& surf_field, const ArrayOfScatteringSpecies& scattering_species, const Agenda& spectral_propmat_agenda, - const Stokvec& transmitted, - const RadarRangeMode coordinate_mode, - const String& unit, - const Numeric ze_tref, - const Numeric k2, - const Numeric dbze_min, - const Numeric pext_scaling, - const ArrayOfString& aux_vars, + const Settings& settings, + const ArrayOfRadarAuxiliary& aux_vars, const JacobianTargets& jac_targets) { const Size ny = sensor.size(); - const Size nq = jac_targets.target_count(); RadarResult out{.measurement = Vector(ny, 0.0), .jacobian = Matrix(ny, jac_targets.x_size(), 0.0), - .auxiliary = Matrix(aux_vars.size(), ny, 0.0)}; - std::vector seen(ny, false); + .auxiliary = Matrix(aux_vars.size(), ny, 0.0), + .time = Vector(ny, 0.0)}; + std::vector covered(ny, false); Vector scalar_weight_sum(ny, 0.0); for (const auto& simulation : simulations) { + if (simulation.samples.empty()) continue; + const auto& frequency = *simulation.frequency; - const auto& poslos = *simulation.poslos; const auto& ray_path = simulation.ray_path; + const auto& samples = simulation.samples; + const Size nf = frequency.size(); + const Size ngate = simulation.gates.size(); const Size iposlos = simulation.iposlos; - if (ray_path.size() < 2) continue; - - const OpticalProfile optical = optical_profile( - ws, frequency, ray_path, atm_field, scattering_species, spectral_propmat_agenda, jac_targets, pext_scaling); - const Size np = ray_path.size(); - const Size nf = frequency.size(); - - std::vector attenuated(np, StokvecVector(nf)); - std::vector unattenuated(np, StokvecVector(nf)); - std::vector gas_ext(np, StokvecVector(nf)); - std::vector particle_ext(np, StokvecVector(nf)); - std::vector outgoing_transmission(nf, Muelmat::id()); - std::vector return_transmission(nf, Muelmat::id()); - std::vector outgoing_transmission_jac(nf, MuelmatMatrix(np, nq, Muelmat{0.0})); - std::vector return_transmission_jac(nf, MuelmatMatrix(np, nq, Muelmat{0.0})); - std::vector> attenuated_jac( - np * nq, std::vector(np, StokvecVector(nf, Stokvec{0.0}))); - - for (Size ip = 0; ip < np; ++ip) { - if (ip > 0) { - const Numeric ds = path::distance(ray_path[ip - 1].pos, ray_path[ip].pos, surf_field.ellipsoid); - for (Size iv = 0; iv < nf; ++iv) { - const Muelmat previous_outgoing = outgoing_transmission[iv]; - const Muelmat previous_return = return_transmission[iv]; - const rtepack::tran outgoing_segment( - optical.outgoing_extinction[ip - 1][iv], optical.outgoing_extinction[ip][iv], ds); - const rtepack::tran return_segment( - optical.return_extinction[ip][iv], optical.return_extinction[ip - 1][iv], ds); - const Muelmat outgoing_segment_transmission = outgoing_segment(); - const Muelmat return_segment_transmission = return_segment(); - for (Size jp = 0; jp < np; ++jp) - for (Size iq = 0; iq < nq; ++iq) { - outgoing_transmission_jac[iv][jp, iq] = - outgoing_segment_transmission * outgoing_transmission_jac[iv][jp, iq]; - return_transmission_jac[iv][jp, iq] = return_transmission_jac[iv][jp, iq] * return_segment_transmission; - } - for (Size iq = 0; iq < nq; ++iq) { - outgoing_transmission_jac[iv][ip - 1, iq] += - outgoing_segment.deriv(outgoing_segment_transmission, - optical.outgoing_extinction[ip - 1][iv], - optical.outgoing_extinction[ip][iv], - optical.outgoing_extinction_jac[ip - 1][iq, iv], - ds, - 0.0) * - previous_outgoing; - outgoing_transmission_jac[iv][ip, iq] += outgoing_segment.deriv(outgoing_segment_transmission, - optical.outgoing_extinction[ip - 1][iv], - optical.outgoing_extinction[ip][iv], - optical.outgoing_extinction_jac[ip][iq, iv], - ds, - 0.0) * - previous_outgoing; - return_transmission_jac[iv][ip - 1, iq] += - previous_return * return_segment.deriv(return_segment_transmission, - optical.return_extinction[ip][iv], - optical.return_extinction[ip - 1][iv], - optical.return_extinction_jac[ip - 1][iq, iv], - ds, - 0.0); - return_transmission_jac[iv][ip, iq] += - previous_return * return_segment.deriv(return_segment_transmission, - optical.return_extinction[ip][iv], - optical.return_extinction[ip - 1][iv], - optical.return_extinction_jac[ip][iq, iv], - ds, - 0.0); - } - outgoing_transmission[iv] = outgoing_segment_transmission * previous_outgoing; - return_transmission[iv] = previous_return * return_segment_transmission; - } - } - for (Size iv = 0; iv < nf; ++iv) { - unattenuated[ip][iv] = optical.backscatter[ip][iv] * transmitted; - attenuated[ip][iv] = rtepack::radar_return( - return_transmission[iv], optical.backscatter[ip][iv], outgoing_transmission[iv], transmitted); - gas_ext[ip][iv] = Stokvec{optical.gas_extinction[ip][iv].A(), 0.0, 0.0, 0.0}; - particle_ext[ip][iv] = Stokvec{pext_scaling * optical.particle_extinction[ip][iv].A(), 0.0, 0.0, 0.0}; - for (Size jp = 0; jp < np; ++jp) { - for (Size iq = 0; iq < nq; ++iq) { - const Muelmat dbackscatter = jp == ip ? optical.backscatter_jac[ip][iq, iv] : Muelmat{0.0}; - const Stokvec derivative = rtepack::radar_return_derivative(return_transmission[iv], - return_transmission_jac[iv][jp, iq], - optical.backscatter[ip][iv], - dbackscatter, - outgoing_transmission[iv], - outgoing_transmission_jac[iv][jp, iq], - transmitted); - attenuated_jac[jp * nq + iq][ip][iv] = derivative; - } - } - } - } - const Vector coordinate = range_coordinate(ray_path, poslos[iposlos].pos, surf_field.ellipsoid, coordinate_mode); + const radar::GasModel gas = [&](PropmatVector& propmat, + PropmatMatrix& propmat_jac, + const PropagationPathPoint& point, + const AtmPoint& atm_point) { + StokvecVector srcvec; + StokvecMatrix srcvec_jac; + spectral_propmat_agendaExecute(ws, + propmat, + srcvec, + propmat_jac, + srcvec_jac, + frequency, + Vector3{0.0, 0.0, 0.0}, + jac_targets, + {}, + point, + atm_point, + spectral_propmat_agenda); + }; + const radar::Optics optics = radar::optics( + ray_path, simulation.np, atm_field, frequency, scattering_species, jac_targets, settings.pext_scaling, gas); + + StokvecMatrix signal; + StokvecTensor4 signal_jac; + radar::single_scattering( + signal, signal_jac, settings.transmitted, optics, simulation.range.layer, samples, ngate, settings.rte_option); + + Vector factor(nf); + for (Size iv = 0; iv < nf; ++iv) + factor[iv] = radar::unit_factor(settings.unit, frequency[iv], settings.ze_tref, settings.k2); + + const Vector coverage = radar::gate_coverage(samples, ngate); + const Vector time = radar::gate_sum(samples, simulation.range.time, ngate); const Numeric background = ray_path.back().los_type == PathPositionType::space ? 0.0 : (ray_path.back().los_type == PathPositionType::surface ? 1.0 : 2.0); - for (Size iy = 0; iy < ny; ++iy) { + const auto contains = [&](const RadarAuxiliary key) { return stdr::contains(aux_vars, key); }; + + StokvecMatrix backscatter(ngate, nf, Stokvec{}); + if (contains(RadarAuxiliary::Backscattering)) { + for (const auto& [layer, x, weight, gate] : samples) + for (Size iv = 0; iv < nf; ++iv) + backscatter[gate, iv] += + weight * (((1.0 - x) * optics.Z[layer - 1][iv] + x * optics.Z[layer][iv]) * settings.transmitted); + } + + Matrix gas_extinction(ngate, nf, 0.0); + if (contains(RadarAuxiliary::AbsSpeciesExtinction)) + for (Size iv = 0; iv < nf; ++iv) + gas_extinction[joker, iv] = radar::gate_sum(samples, optics.gas_extinction[iv], ngate); + + Matrix particle_extinction(ngate, nf, 0.0); + if (contains(RadarAuxiliary::ParticleExtinction)) + for (Size iv = 0; iv < nf; ++iv) + particle_extinction[joker, iv] = radar::gate_sum(samples, optics.particle_extinction[iv], ngate); + + StokvecVector bin(nf); + for (Size ig = 0; ig < ngate; ++ig) { + if (coverage[ig] == 0.0) continue; + + const Size iy = simulation.gates[ig]; const auto& obsel = sensor[iy]; - if (not obsel.same_freqs(frequency)) continue; - - StokvecVector bin_signal(nf, Stokvec{0.0}); - StokvecVector bin_backscatter(nf, Stokvec{0.0}); - StokvecVector bin_gas(nf, Stokvec{0.0}); - StokvecVector bin_particle(nf, Stokvec{0.0}); - std::vector bin_jac(np * nq, StokvecVector(nf, Stokvec{0.0})); - bool inside = false; - for (Size iv = 0; iv < nf; ++iv) { - std::vector a(np), b(np), g(np), p(np); - for (Size ip = 0; ip < np; ++ip) { - a[ip] = attenuated[ip][iv]; - b[ip] = unattenuated[ip][iv]; - g[ip] = gas_ext[ip][iv]; - p[ip] = particle_ext[ip][iv]; - } - if (auto avg = bin_average(coordinate, a, limits[iy, 0], limits[iy, 1])) { - inside = true; - bin_signal[iv] = *avg; - bin_backscatter[iv] = *bin_average(coordinate, b, limits[iy, 0], limits[iy, 1]); - bin_gas[iv] = *bin_average(coordinate, g, limits[iy, 0], limits[iy, 1]); - bin_particle[iv] = *bin_average(coordinate, p, limits[iy, 0], limits[iy, 1]); - - const Numeric factor = unit == "1" ? 1.0 : ze_factor(frequency[iv], ze_tref, k2); - bin_signal[iv] *= factor; - bin_backscatter[iv] *= factor; - for (Size jp = 0; jp < np; ++jp) { - for (Size iq = 0; iq < nq; ++iq) { - std::vector derivative(np); - for (Size kp = 0; kp < np; ++kp) derivative[kp] = attenuated_jac[jp * nq + iq][kp][iv]; - bin_jac[jp * nq + iq][iv] = factor * *bin_average(coordinate, derivative, limits[iy, 0], limits[iy, 1]); - } - } - } - } - if (not inside) { - if (not seen[iy]) { - out.measurement[iy] = std::numeric_limits::quiet_NaN(); - for (Size ia = 0; ia < aux_vars.size(); ++ia) - out.auxiliary[ia, iy] = std::numeric_limits::quiet_NaN(); - } - continue; - } + covered[iy] = true; - if (not seen[iy] or std::isnan(out.measurement[iy])) { - out.measurement[iy] = 0.0; - for (Size ia = 0; ia < aux_vars.size(); ++ia) out.auxiliary[ia, iy] = 0.0; - } - const Numeric signal = obsel.sumup(bin_signal, iposlos); - out.measurement[iy] += signal; - for (const auto& target : jac_targets.atm) { - const auto& field = atm_field[target.type]; - for (Size jp = 0; jp < np; ++jp) { - const Numeric local = obsel.sumup(bin_jac[jp * nq + target.target_pos], iposlos); - for (const auto& [index, weight] : field.flat_weight(ray_path[jp].pos)) - out.jacobian[iy, target.x_start + index] += weight * local; - } - } + for (Size iv = 0; iv < nf; ++iv) bin[iv] = factor[iv] * signal[ig, iv]; + out.measurement[iy] += obsel.sumup(bin, iposlos); + // Extinction, background and time are means over the covered part of the gate const Numeric scalar_weight = obsel.sumup(Stokvec{1.0, 0.0, 0.0, 0.0}, iposlos); scalar_weight_sum[iy] += scalar_weight; + out.time[iy] += scalar_weight * time[ig] / coverage[ig]; for (Size ia = 0; ia < aux_vars.size(); ++ia) { - if (aux_vars[ia] == "Radiative background") { - out.auxiliary[ia, iy] += scalar_weight * background; - } else if (aux_vars[ia] == "Backscattering") { - out.auxiliary[ia, iy] += obsel.sumup(bin_backscatter, iposlos); - } else if (aux_vars[ia] == "Abs species extinction") { - out.auxiliary[ia, iy] += obsel.sumup(bin_gas, iposlos); - } else if (aux_vars[ia] == "Particle extinction") { - out.auxiliary[ia, iy] += obsel.sumup(bin_particle, iposlos); + switch (aux_vars[ia]) { + case RadarAuxiliary::RadiativeBackground: out.auxiliary[ia, iy] += scalar_weight * background; break; + case RadarAuxiliary::Backscattering: + for (Size iv = 0; iv < nf; ++iv) bin[iv] = factor[iv] * backscatter[ig, iv]; + out.auxiliary[ia, iy] += obsel.sumup(bin, iposlos); + break; + case RadarAuxiliary::AbsSpeciesExtinction: + for (Size iv = 0; iv < nf; ++iv) bin[iv] = Stokvec{gas_extinction[ig, iv] / coverage[ig]}; + out.auxiliary[ia, iy] += obsel.sumup(bin, iposlos); + break; + case RadarAuxiliary::ParticleExtinction: + for (Size iv = 0; iv < nf; ++iv) bin[iv] = Stokvec{particle_extinction[ig, iv] / coverage[ig]}; + out.auxiliary[ia, iy] += obsel.sumup(bin, iposlos); + break; + } + } + } + + for (const auto& target : jac_targets.atm) { + const auto& field = atm_field[target.type]; + for (Size jp = 0; jp < simulation.np; ++jp) { + const auto weights = field.flat_weight(ray_path[jp].pos); + for (Size ig = 0; ig < ngate; ++ig) { + if (coverage[ig] == 0.0) continue; + + const Size iy = simulation.gates[ig]; + for (Size iv = 0; iv < nf; ++iv) bin[iv] = factor[iv] * signal_jac[ig, jp, target.target_pos, iv]; + const Numeric local = sensor[iy].sumup(bin, iposlos); + for (const auto& [index, weight] : weights) out.jacobian[iy, target.x_start + index] += weight * local; } } - seen[iy] = true; } } - for (Size ia = 0; ia < aux_vars.size(); ++ia) { - if (aux_vars[ia] == "Backscattering") continue; - for (Size iy = 0; iy < ny; ++iy) - if (seen[iy] and scalar_weight_sum[iy] != 0.0) out.auxiliary[ia, iy] /= scalar_weight_sum[iy]; + constexpr Numeric nan = std::numeric_limits::quiet_NaN(); + for (Size iy = 0; iy < ny; ++iy) { + if (not covered[iy]) { + out.measurement[iy] = nan; + out.time[iy] = nan; + out.auxiliary[joker, iy] = nan; + continue; + } + + if (scalar_weight_sum[iy] == 0.0) continue; + out.time[iy] /= scalar_weight_sum[iy]; + for (Size ia = 0; ia < aux_vars.size(); ++ia) + if (aux_vars[ia] != RadarAuxiliary::Backscattering) out.auxiliary[ia, iy] /= scalar_weight_sum[iy]; } - if (unit == "dBZe") { - const Numeric ze_min = std::pow(10.0, dbze_min / 10.0); + if (settings.unit == RadarUnit::dBZe) { + const Numeric ze_min = std::pow(10.0, settings.dbze_min / 10.0); for (Size iy = 0; iy < ny; ++iy) { if (std::isnan(out.measurement[iy])) continue; if (out.measurement[iy] <= ze_min) { - out.measurement[iy] = dbze_min; + out.measurement[iy] = settings.dbze_min; out.jacobian[iy, joker] = 0.0; } else { - const Numeric factor = 10.0 / (std::log(10.0) * out.measurement[iy]); - out.jacobian[iy, joker] *= factor; + out.jacobian[iy, joker] *= 10.0 / (std::log(10.0) * out.measurement[iy]); out.measurement[iy] = 10.0 * std::log10(out.measurement[iy]); } } for (Size ia = 0; ia < aux_vars.size(); ++ia) { - if (aux_vars[ia] != "Backscattering") continue; + if (aux_vars[ia] != RadarAuxiliary::Backscattering) continue; for (Size iy = 0; iy < ny; ++iy) { auto& value = out.auxiliary[ia, iy]; - if (not std::isnan(value)) value = value <= ze_min ? dbze_min : 10.0 * std::log10(value); + if (not std::isnan(value)) value = value <= ze_min ? settings.dbze_min : 10.0 * std::log10(value); } } } @@ -517,6 +276,15 @@ RadarResult radar_forward(const Workspace& ws, return out; } +//! Apply the target transformations to a Jacobian with regards to the untransformed model state +void transform_jacobian(Matrix& jacobian, const JacobianTargets& jac_targets, const AtmField& atm_field) { + const Size nx = jac_targets.x_size(); + if (nx == 0) return; + + Vector state(nx, 0.0); + for (const auto& target : jac_targets.atm) target.update_state(state, atm_field); + for (const auto& target : jac_targets.atm) target.update_jac(jacobian, state, atm_field); +} } // namespace void measurement_sensorAddSimpleRadar(ArrayOfSensorObsel& measurement_sensor, @@ -577,6 +345,7 @@ void measurement_vecFromRadarSingleScattering(const Workspace& ws Vector& measurement_vec, Matrix& measurement_jac, Matrix& radar_aux, + Vector& radar_time, const ArrayOfSensorObsel& measurement_sensor, const Matrix& radar_range_limits, const JacobianTargets& jac_targets, @@ -585,86 +354,44 @@ void measurement_vecFromRadarSingleScattering(const Workspace& ws const ArrayOfScatteringSpecies& scattering_species, const Agenda& ray_path_observer_agenda, const Agenda& spectral_propmat_agenda, + const TransmittanceOption& rte_option, const Stokvec& transmitted_stokes, - const String& range_mode_string, - const String& unit, + const RadarRangeMode& range_mode, + const RadarUnit& unit, const Numeric& ze_tref, const Numeric& k2, const Numeric& dbze_min, const Numeric& pext_scaling, - const ArrayOfString& aux_vars) try { + const ArrayOfRadarAuxiliary& aux_vars) try { ARTS_TIME_REPORT - const Size ny = measurement_sensor.size(); - if (ny == 0) arts::user_error("measurement_sensor cannot be empty for a radar calculation"); - if (radar_range_limits.nrows() != static_cast(ny) or radar_range_limits.ncols() != 2) - arts::user_error("radar_range_limits must have shape [{}, 2], got {:B,}", ny, radar_range_limits.shape()); - for (Size iy = 0; iy < ny; ++iy) { - measurement_sensor[iy].check(); - if (not(radar_range_limits[iy, 0] < radar_range_limits[iy, 1])) - arts::user_error("Radar range limits in row {} are not strictly increasing: [{}, {}]", - iy, - radar_range_limits[iy, 0], - radar_range_limits[iy, 1]); - } - if (transmitted_stokes.I() != 1.0) - arts::user_error("The transmitted Stokes I component must equal one, got {}", transmitted_stokes.I()); - if (unit != "1" and unit != "Ze" and unit != "dBZe") - arts::user_error("Radar unit must be \"1\", \"Ze\", or \"dBZe\"; got \"{}\"", unit); - if (pext_scaling < 0.0 or pext_scaling > 2.0) - arts::user_error("pext_scaling must be between 0 and 2; got {}", pext_scaling); - if (scattering_species.species.empty()) arts::user_error("Radar single scattering requires scattering species"); - if (not jac_targets.surf.empty() or not jac_targets.subsurf.empty() or not jac_targets.line.empty() or - not jac_targets.sensor.empty() or not jac_targets.error.empty()) - arts::user_error("The radar forward model currently supports atmospheric Jacobian targets only"); - for (const auto& name : aux_vars) { - if (name != "Radiative background" and name != "Backscattering" and name != "Abs species extinction" and - name != "Particle extinction") - arts::user_error("Unknown radar auxiliary variable \"{}\"", name); - } - - const RadarRangeMode coordinate_mode = range_mode(range_mode_string, radar_range_limits); - const SensorSimulations collected = collect_simulations(measurement_sensor); - std::vector simulations; - simulations.reserve(collected.size()); - for (const auto& simulation : collected) { - const auto& [frequency, poslos, iposlos] = simulation; - PreparedSimulation prepared{.frequency = &frequency, .poslos = &poslos, .iposlos = iposlos}; - ray_path_observer_agendaExecute( - ws, prepared.ray_path, poslos[iposlos].pos, poslos[iposlos].los, ray_path_observer_agenda); - simulations.push_back(std::move(prepared)); - } - const auto calculate = [&](const AtmField& field) { - return radar_forward(ws, - measurement_sensor, - radar_range_limits, - simulations, - field, - surf_field, - scattering_species, - spectral_propmat_agenda, - transmitted_stokes, - coordinate_mode, - unit, - ze_tref, - k2, - dbze_min, - pext_scaling, - aux_vars, - jac_targets); - }; - - RadarResult base = calculate(atm_field); - measurement_vec = std::move(base.measurement); - measurement_jac = std::move(base.jacobian); - radar_aux = std::move(base.auxiliary); - - const Size nx = jac_targets.x_size(); - if (nx == 0) return; - - Vector state(nx, 0.0); - for (const auto& target : jac_targets.atm) target.update_state(state, atm_field); - for (const auto& target : jac_targets.atm) target.update_jac(measurement_jac, state, atm_field); + const Settings settings{.transmitted = transmitted_stokes, + .range_mode = resolved_range_mode(range_mode, radar_range_limits), + .unit = unit, + .ze_tref = ze_tref, + .k2 = k2, + .dbze_min = dbze_min, + .pext_scaling = pext_scaling, + .rte_option = rte_option}; + check_settings(measurement_sensor, radar_range_limits, jac_targets, scattering_species, settings); + + const auto simulations = prepare_simulations( + ws, measurement_sensor, radar_range_limits, settings.range_mode, surf_field, ray_path_observer_agenda); + RadarResult result = radar_forward(ws, + measurement_sensor, + simulations, + atm_field, + scattering_species, + spectral_propmat_agenda, + settings, + aux_vars, + jac_targets); + + measurement_vec = std::move(result.measurement); + measurement_jac = std::move(result.jacobian); + radar_aux = std::move(result.auxiliary); + radar_time = std::move(result.time); + transform_jacobian(measurement_jac, jac_targets, atm_field); } ARTS_METHOD_ERROR_CATCH @@ -681,9 +408,10 @@ void model_state_vecFromRadarOnionPeeling(const Workspace& ws, const ArrayOfScatteringSpecies& scattering_species, const Agenda& ray_path_observer_agenda, const Agenda& spectral_propmat_agenda, + const TransmittanceOption& rte_option, const Stokvec& transmitted_stokes, - const String& range_mode_string, - const String& unit, + const RadarRangeMode& range_mode, + const RadarUnit& unit, const Numeric& ze_tref, const Numeric& k2, const Numeric& dbze_min, @@ -698,13 +426,20 @@ void model_state_vecFromRadarOnionPeeling(const Workspace& ws, ARTS_TIME_REPORT const Size ny = measurement_sensor.size(); const Size nx = jac_targets.x_size(); - if (radar_range_limits.nrows() != static_cast(ny) or radar_range_limits.ncols() != 2) - arts::user_error("radar_range_limits must have shape [{}, 2], got {:B,}", ny, radar_range_limits.shape()); + + const Settings settings{.transmitted = transmitted_stokes, + .range_mode = resolved_range_mode(range_mode, radar_range_limits), + .unit = unit, + .ze_tref = ze_tref, + .k2 = k2, + .dbze_min = dbze_min, + .pext_scaling = pext_scaling, + .rte_option = rte_option}; + check_settings(measurement_sensor, radar_range_limits, jac_targets, scattering_species, settings); + if (measurement_vec.size() != ny) + arts::user_error("measurement_vec has {} elements but there are {} radar observations", measurement_vec.size(), ny); if (nx == 0 or jac_targets.atm.empty()) arts::user_error("Radar onion peeling requires at least one finalized atmospheric Jacobian target"); - if (not jac_targets.surf.empty() or not jac_targets.subsurf.empty() or not jac_targets.line.empty() or - not jac_targets.sensor.empty() or not jac_targets.error.empty()) - arts::user_error("Radar onion peeling currently changes atmospheric state coordinates only"); if (state_min > state_max) arts::user_error("state_min ({}) exceeds state_max ({})", state_min, state_max); if (max_step <= 0.0) arts::user_error("max_step must be positive, got {}", max_step); if (tolerance <= 0.0) arts::user_error("tolerance must be positive, got {}", tolerance); @@ -718,28 +453,23 @@ void model_state_vecFromRadarOnionPeeling(const Workspace& ws, model_state_vec = 0.0; for (const auto& target : jac_targets.atm) target.update_state(model_state_vec, atm_field); + // The paths are traced for every update, since the path agenda may depend + // on the atmosphere that is being peeled. const auto calculate = [&]() { - Matrix radar_aux; - measurement_vecFromRadarSingleScattering(ws, - measurement_vec_fit, - measurement_jac, - radar_aux, - measurement_sensor, - radar_range_limits, - jac_targets, - atm_field, - surf_field, - scattering_species, - ray_path_observer_agenda, - spectral_propmat_agenda, - transmitted_stokes, - range_mode_string, - unit, - ze_tref, - k2, - dbze_min, - pext_scaling, - {}); + const auto simulations = prepare_simulations( + ws, measurement_sensor, radar_range_limits, settings.range_mode, surf_field, ray_path_observer_agenda); + RadarResult result = radar_forward(ws, + measurement_sensor, + simulations, + atm_field, + scattering_species, + spectral_propmat_agenda, + settings, + {}, + jac_targets); + measurement_vec_fit = std::move(result.measurement); + measurement_jac = std::move(result.jacobian); + transform_jacobian(measurement_jac, jac_targets, atm_field); }; calculate(); @@ -752,10 +482,9 @@ void model_state_vecFromRadarOnionPeeling(const Workspace& ws, // blocks have zero cross-sensitivity. std::vector order(ny); stdr::iota(order, Size{0}); - const RadarRangeMode coordinate_mode = range_mode(range_mode_string, radar_range_limits); - const auto distance_key = [&](const Size iy) { + const auto distance_key = [&](const Size iy) { const Numeric center = std::midpoint(radar_range_limits[iy, 0], radar_range_limits[iy, 1]); - if (coordinate_mode != RadarRangeMode::Altitude) return center; + if (settings.range_mode != RadarRangeMode::Altitude) return center; const auto& poslos = measurement_sensor[iy].poslos_grid(); if (poslos.empty()) arts::user_error("Radar observation {} has no position/LOS geometry", iy); return std::abs(center - poslos.front().pos[0]); diff --git a/src/python_interface/CMakeLists.txt b/src/python_interface/CMakeLists.txt index 9339736b10..9f23013131 100644 --- a/src/python_interface/CMakeLists.txt +++ b/src/python_interface/CMakeLists.txt @@ -184,6 +184,7 @@ nanobind_add_module(pyarts_cpp NB_STATIC NB_DOMAIN ARTS ${NANOBIND_LTO} STABLE_A py_geodesy.cpp py_planets.cpp py_montecarlo.cpp + py_radar.cpp py_rng.cpp py_tmatrix.cpp py_polradtran.cpp diff --git a/src/python_interface/py_module.cpp b/src/python_interface/py_module.cpp index 924fb0721f..2c98e83c2f 100644 --- a/src/python_interface/py_module.cpp +++ b/src/python_interface/py_module.cpp @@ -55,6 +55,7 @@ void py_planets(py::module_& m); void py_auto_agenda_operators(py::module_& m); void py_rng(py::module_& m); void py_montecarlo(py::module_& m); +void py_radar(py::module_& m); void py_tmatrix(py::module_& m); void py_polradtran(py::module_& m); void py_rt3(py::module_& m); @@ -129,6 +130,7 @@ NB_MODULE(arts, m) try { py_auto_agenda_operators(m); py_rng(m); py_montecarlo(m); + py_radar(m); py_tmatrix(m); py_polradtran(m); py_rt3(m); diff --git a/src/python_interface/py_radar.cpp b/src/python_interface/py_radar.cpp new file mode 100644 index 0000000000..687dc52b72 --- /dev/null +++ b/src/python_interface/py_radar.cpp @@ -0,0 +1,27 @@ +#include +#include +#include +#include +#include + +#include + +#include "hpy_arts.h" +#include "hpy_vector.h" + +namespace Python { +void py_radar(py::module_& m) try { + auto aora = py::bind_vector(m, "ArrayOfRadarAuxiliary"); + generic_interface(aora); + vector_interface(aora); + aora.def("__init__", [](ArrayOfRadarAuxiliary* self, const std::vector& x) { + new (self) ArrayOfRadarAuxiliary(); + self->reserve(x.size()); + std::transform( + x.begin(), x.end(), std::back_inserter(*self), [](const std::string& s) { return to(s); }); + }); + py::implicitly_convertible, ArrayOfRadarAuxiliary>(); +} catch (std::exception& e) { + throw std::runtime_error(std::format("DEV ERROR:\nCannot initialize radar\n{}", e.what())); +} +} // namespace Python diff --git a/src/tests/test_ai_port_regressions.cc b/src/tests/test_ai_port_regressions.cc index 03fb5d780f..b4c554db9c 100644 --- a/src/tests/test_ai_port_regressions.cc +++ b/src/tests/test_ai_port_regressions.cc @@ -216,16 +216,19 @@ void radar_gates() { antenna.set_pencil_beam(); const AscendingGrid bins{100.0, 200.0, 400.0}; StokvecVector signal, error; + Vector time; const auto run = [&](const Vector& altitudes) { MCRadar(ws, signal, error, + time, atmosphere, surface, species, antenna, path_agenda(altitudes), gas, + TransmittanceOption::constant, 94e9, Vector3{1000.0, 0.0, 0.0}, Vector2{180.0, 0.0}, @@ -235,7 +238,10 @@ void radar_gates() { 1, 1, 0.93, - "1"); + RadarUnit::unit, + RadarRangeMode::Distance, + 273.15, + 1.0); return signal; }; // The coarse segment midpoint (500 m) is beyond the final gate (400 m). @@ -248,6 +254,7 @@ void radar_gates() { close(coarse[i].I() / expected, 1.0, "Radar must integrate every overlapped gate"); close(split[i].I() / coarse[i].I(), 1.0, "Gate splitting must preserve the homogeneous return"); for (Numeric x : error[i]) close(x, 0.0, "Pencil beam is deterministic"); + close(time[i] * Constant::speed_of_light / (2.0 * midpoint), 1.0, "Radar time must be the round trip to the gate"); } } } // namespace diff --git a/src/workspace_groups.cpp b/src/workspace_groups.cpp index 15bbfd6de7..9c4f8d8fe6 100644 --- a/src/workspace_groups.cpp +++ b/src/workspace_groups.cpp @@ -1114,6 +1114,12 @@ of this term multiplied by a negative distance. .invariant_printables = {"{}"}, }; + wsg_data["ArrayOfRadarAuxiliary"] = { + .file = "radar.h", + .desc = "A list of *RadarAuxiliary*\n", + .dim_size = {"{}.size()"}, + }; + wsg_data["XsecRecords"] = { .file = "xsec_fit.h", .desc = R"--(A map from *SpeciesEnum* to *XsecRecord* @@ -1139,6 +1145,8 @@ of this term multiplied by a negative distance. "SurfaceKey", "SubsurfaceKey", "SpectralRadianceUnitType", + "RadarUnit", + "RadarRangeMode", }); agenda_operators(wsg_data); diff --git a/src/workspace_methods.cpp b/src/workspace_methods.cpp index 6de2c76f9b..4150b2ccf3 100644 --- a/src/workspace_methods.cpp +++ b/src/workspace_methods.cpp @@ -6824,36 +6824,70 @@ Neither polarization nor wind calculations are possible with this method. }; wsm_data["MCRadar"] = { - .desc = R"--(Simulates active-radar returns using Monte Carlo antenna sampling. + .desc = R"--(Simulates active-radar returns using Monte Carlo antenna sampling. This ARTS3-native implementation uses *AtmField*, *SurfaceField*, -*ArrayOfScatteringSpecies*, and *ray_path_observer_agenda*. Range-bin edges -are one-way geometric distances in metres. The outputs contain one -*Stokvec* per range bin, and the error contains the component-wise standard -error of the Monte Carlo mean. The implementation currently supports single -scattering (``mc_max_scatorder=1``); the argument is retained so higher orders -can be added without changing the interface. -)--", - .author = {"Richard Larsson", "OpenAI Codex"}, - .out = {"radar_signal", "radar_error"}, - .in = {"atm_field", - "surf_field", - "scat_species", - "mc_antenna", - "ray_path_observer_agenda", - "spectral_propmat_agenda"}, - .gin = {"frequency", - "sensor_pos", - "sensor_los", - "mc_y_tx", - "range_bins", - "mc_seed", - "mc_max_iter", - "mc_max_scatorder", - "k2", - "unit"}, - .gin_type = - {"Numeric", "Vector3", "Vector2", "Stokvec", "AscendingGrid", "Index", "Index", "Index", "Numeric", "String"}, +*ArrayOfScatteringSpecies*, and *ray_path_observer_agenda*. Every sampled +antenna direction is solved with the same range gates, optical properties, +and two-way rtepack transmission as *measurement_vecFromRadarSingleScattering*, +so a pencil beam gives the same result. + +``range_bins`` are the gate edges in the coordinate of ``range_mode``. The +range includes the distance from the radar to the first path point, as for a +radar outside the atmosphere. Gates are averaged over their full width. + +The mean over the sampled directions is weighted by the antenna gain, so a +uniformly filled beam returns the same value as a pencil beam. The outputs +contain one *Stokvec* per range bin, and the error contains the +component-wise standard error of the weighted mean. *radar_time* is the mean +round-trip time of the part of each bin that the paths cover. Bins that no +path reaches are NaN. + +``unit`` may be ``"1"`` or ``"Ze"``. A negative ``k2`` derives the +liquid-water dielectric factor at ``ze_tref`` using the Liebe-93 +parameterization. ``pext_scaling`` scales only particulate extinction, not +particle backscatter. *rte_option* selects the transmission of every layer, +and of the part of a layer inside a gate. + +The implementation currently supports single scattering +(``mc_max_scatorder=1``); the argument is retained so higher orders can be +added without changing the interface. +)--", + .author = {"Richard Larsson", "OpenAI Codex"}, + .out = {"radar_signal", "radar_error", "radar_time"}, + .in = {"atm_field", + "surf_field", + "scat_species", + "mc_antenna", + "ray_path_observer_agenda", + "spectral_propmat_agenda", + "rte_option"}, + .gin = {"frequency", + "sensor_pos", + "sensor_los", + "mc_y_tx", + "range_bins", + "mc_seed", + "mc_max_iter", + "mc_max_scatorder", + "k2", + "unit", + "range_mode", + "ze_tref", + "pext_scaling"}, + .gin_type = {"Numeric", + "Vector3", + "Vector2", + "Stokvec", + "AscendingGrid", + "Index", + "Index", + "Index", + "Numeric", + "RadarUnit", + "RadarRangeMode", + "Numeric", + "Numeric"}, .gin_value = {std::nullopt, std::nullopt, std::nullopt, @@ -6863,17 +6897,23 @@ can be added without changing the interface. Index{1000}, Index{1}, Numeric{0.93}, - String{"1"}}, + RadarUnit::unit, + RadarRangeMode::Distance, + Numeric{273.15}, + Numeric{1.0}}, .gin_desc = {"Radar frequency [Hz].", "Sensor position [alt, lat, lon].", "Sensor boresight [zenith, azimuth].", "Transmitted Stokes vector.", - "One-way range-bin edges [m].", + "Range-bin edges in the coordinate of range_mode.", "Random seed.", "Number of sampled antenna rays.", "Maximum scattering order (currently must be 1).", - "Reference dielectric factor squared for Ze conversion.", - "Output unit: '1' or 'Ze'."}, + "Reference dielectric factor squared; negative selects Liebe-93.", + "Output unit: 1 or Ze.", + "Range coordinate: Altitude, Distance, RoundTripTime, or Legacy.", + "Liquid-water reference temperature [K] for automatic k2.", + "Multiplicative factor for particulate extinction (0 to 2)."}, .pass_workspace = true, }; @@ -6920,7 +6960,15 @@ frequency/polarization/range-bin arrangement. ``range_mode`` may be ``"Altitude"`` [m], ``"Distance"`` [m one way], ``"RoundTripTime"`` [s, including path group refractive indices], or ``"Legacy"`` (the ARTS2 rule: altitude if the largest edge exceeds 1, -otherwise round-trip time). +otherwise round-trip time). Distance and round-trip time include the leg from +the radar to the first path point, as for a radar outside the atmosphere. + +Every gate is integrated layer by layer. The return at the midpoint of the +overlap of a layer and a gate uses the transmission of *rte_option* from the +layer edge, so the Magnus operator also applies to partial layers. Gates are +averaged over their full width, so the parts of a gate without any path +contribute nothing. Gates that no path reaches are NaN. *radar_time* is the +mean round-trip time of the covered part of each gate. The output unit may be ``"1"`` (backscatter coefficient, 1/(m sr)), ``"Ze"`` (mm6/m3), or ``"dBZe"``. A negative ``k2`` derives the liquid-water @@ -6936,13 +6984,12 @@ come from *spectral_propmat_agenda*, scattering derivatives come from *scat_species*, and their ordered two-way transmission product is propagated analytically through rtepack. There is no production finite-difference loop. -Allowed auxiliary quantities are ``"Radiative background"``, -``"Backscattering"``, ``"Abs species extinction"``, and -``"Particle extinction"``. They are returned as rows of *radar_aux* in the -requested order. +The auxiliary quantities are listed by *RadarAuxiliary*. They are returned +as rows of *radar_aux* in the requested order. Extinction and background are +means over the covered part of each gate. )--", .author = {"Patrick Eriksson", "OpenAI Codex"}, - .out = {"measurement_vec", "measurement_jac", "radar_aux"}, + .out = {"measurement_vec", "measurement_jac", "radar_aux", "radar_time"}, .in = {"measurement_sensor", "radar_range_limits", "jac_targets", @@ -6950,17 +6997,25 @@ requested order. "surf_field", "scat_species", "ray_path_observer_agenda", - "spectral_propmat_agenda"}, + "spectral_propmat_agenda", + "rte_option"}, .gin = {"transmitted_stokes", "range_mode", "unit", "ze_tref", "k2", "dbze_min", "pext_scaling", "aux_vars"}, - .gin_type = {"Stokvec", "String", "String", "Numeric", "Numeric", "Numeric", "Numeric", "ArrayOfString"}, + .gin_type = {"Stokvec", + "RadarRangeMode", + "RadarUnit", + "Numeric", + "Numeric", + "Numeric", + "Numeric", + "ArrayOfRadarAuxiliary"}, .gin_value = {Stokvec{1.0, 1.0, 0.0, 0.0}, - String{"Legacy"}, - String{"1"}, + RadarRangeMode::Legacy, + RadarUnit::unit, Numeric{273.15}, Numeric{-1.0}, Numeric{-99.0}, Numeric{1.0}, - ArrayOfString{}}, + ArrayOfRadarAuxiliary{}}, .gin_desc = {"Transmitted Stokes vector; its I component must equal one.", "Range coordinate: Altitude, Distance, RoundTripTime, or Legacy.", "Output unit: 1, Ze, or dBZe.", @@ -7008,7 +7063,8 @@ and atmospheric state chosen by the caller. "surf_field", "scat_species", "ray_path_observer_agenda", - "spectral_propmat_agenda"}, + "spectral_propmat_agenda", + "rte_option"}, .gin = {"transmitted_stokes", "range_mode", "unit", @@ -7024,8 +7080,8 @@ and atmospheric state chosen by the caller. "max_iterations", "max_sweeps"}, .gin_type = {"Stokvec", - "String", - "String", + "RadarRangeMode", + "RadarUnit", "Numeric", "Numeric", "Numeric", @@ -7038,8 +7094,8 @@ and atmospheric state chosen by the caller. "Index", "Index"}, .gin_value = {Stokvec{1.0, 1.0, 0.0, 0.0}, - String{"Legacy"}, - String{"1"}, + RadarRangeMode::Legacy, + RadarUnit::unit, Numeric{273.15}, Numeric{-1.0}, Numeric{-99.0}, diff --git a/src/workspace_variables.cpp b/src/workspace_variables.cpp index f43238f900..4227a6ccdb 100644 --- a/src/workspace_variables.cpp +++ b/src/workspace_variables.cpp @@ -221,6 +221,16 @@ When Bath is selected, all species are used. Otherwise, this variable should co .type = "StokvecVector", }; + wsv_data["radar_time"] = { + .desc = R"--(Mean round-trip time of every radar range gate [s]. + +The time includes the path group refractive indices and the leg from the +radar to the first path point. It is the mean over the part of the gate that +the paths cover, and NaN for gates that no path reaches. +)--", + .type = "Vector", + }; + wsv_data["radar_range_limits"] = { .desc = R"--(Range-gate limits for active-radar measurements. diff --git a/tests/core/scat/mc_radar.py b/tests/core/scat/mc_radar.py index dc669b2ded..8714f8e8e6 100644 --- a/tests/core/scat/mc_radar.py +++ b/tests/core/scat/mc_radar.py @@ -56,6 +56,7 @@ ws.MCRadar(**kwargs) signal = np.array(ws.radar_signal) error = np.array(ws.radar_error) +gaussian_time = np.array(ws.radar_time) assert isinstance(ws.radar_signal, A.StokvecVector) assert isinstance(ws.radar_error, A.StokvecVector) @@ -65,6 +66,8 @@ assert np.all(np.isfinite(error)) assert np.all(error >= 0.0) assert np.any(signal[:, 0] > 0.0) +assert gaussian_time.shape == (20,) +assert np.all(np.diff(gaussian_time) > 0.0) # The same seed must give bitwise-identical Monte Carlo output. ws.MCRadar(**kwargs) @@ -83,13 +86,16 @@ # In the pencil-beam, single-scattering limit, MCRadar must reproduce the # ARTS3 deterministic rtepack radar result. One history is sufficient because -# there is then no random angular sampling. The small tolerance reflects the -# midpoint quadrature used by MCRadar versus path interpolation in the -# deterministic method. +# there is then no random angular sampling. Both methods share the gate +# quadrature and the two-way transmission, so they agree to rounding, also in +# Ze. ws.mc_antenna.set_pencil_beam() ws.max_stepsize = 100.0 ws.MCRadar(**(kwargs | {"mc_max_iter": 1})) mc_pencil = np.asarray(ws.radar_signal)[:, 0].copy() +mc_time = np.asarray(ws.radar_time).copy() +ws.MCRadar(**(kwargs | {"mc_max_iter": 1, "unit": "Ze"})) +mc_pencil_ze = np.asarray(ws.radar_signal)[:, 0].copy() ws.freq_grid = [kwargs["frequency"]] ws.measurement_sensorInit() @@ -101,14 +107,32 @@ range_bins=kwargs["range_bins"], ) ws.jac_targetsOff() -ws.measurement_vecFromRadarSingleScattering( +deterministic_kwargs = dict( transmitted_stokes=kwargs["mc_y_tx"], range_mode="Distance", - unit="1", pext_scaling=1.0, + k2=0.93, ) -deterministic = np.asarray(ws.measurement_vec) -np.testing.assert_allclose(mc_pencil, deterministic, rtol=2e-3, atol=1e-15) +ws.measurement_vecFromRadarSingleScattering(**deterministic_kwargs, unit="1") +deterministic = np.asarray(ws.measurement_vec).copy() +np.testing.assert_allclose(mc_pencil, deterministic, rtol=1e-10, atol=1e-30) +np.testing.assert_allclose(mc_time, ws.radar_time, rtol=1e-12) + +ws.measurement_vecFromRadarSingleScattering(**deterministic_kwargs, unit="Ze") +np.testing.assert_allclose(mc_pencil_ze, ws.measurement_vec, rtol=1e-10, atol=1e-30) + +# The round-trip time of a vertical path in vacuum is twice the distance to +# the gate centre divided by the speed of light. +centres = 0.5 * (kwargs["range_bins"][1:] + kwargs["range_bins"][:-1]) +np.testing.assert_allclose(mc_time, 2.0 * centres / A.constants.c, rtol=1e-12) + +# A narrow beam is uniformly filled by this horizontally homogeneous cloud. +# The antenna-weighted mean must then equal the pencil-beam return. +ws.mc_antenna.set_gaussian_fwhm(0.6, 0.6) +ws.MCRadar(**(kwargs | {"mc_max_iter": 16})) +gaussian = np.asarray(ws.radar_signal)[:, 0] +inside = (centres > 6e3) & (centres < 14e3) +np.testing.assert_allclose(gaussian[inside], mc_pencil[inside], rtol=1e-4) # Invalid higher scattering orders are rejected explicitly until the collision # sampler is added, rather than silently returning a single-scattering result. diff --git a/tests/core/scat/radar_single_scattering.py b/tests/core/scat/radar_single_scattering.py index f91ebe6494..e068beaffc 100644 --- a/tests/core/scat/radar_single_scattering.py +++ b/tests/core/scat/radar_single_scattering.py @@ -108,6 +108,18 @@ def manual_atmospheric_derivative(workspace, key, flat_index, step, kwargs): # polarization selector, matching old yRadar behavior. assert np.allclose(aux[3, interior], 0.0) # pext_scaling is zero +# The round-trip time of every gate of this vertical vacuum path. +centres = 0.5 * (range_bins[1:] + range_bins[:-1]) +np.testing.assert_allclose( + ws.radar_time, 2.0 * (20e3 - centres) / A.constants.c, rtol=1e-12 +) + +# The Magnus operator differs only for non-commuting propagation matrices. +ws.rte_option = "magop" +ws.measurement_vecFromRadarSingleScattering(**common, unit="1") +np.testing.assert_allclose(ws.measurement_vec, native, rtol=1e-12, atol=1e-30) +ws.rte_option = "linsrc" + # Ze and dBZe conversions, including the old clipping convention. k2 = 0.93 ws.measurement_vecFromRadarSingleScattering(**common, unit="Ze", k2=k2) From fa688ac0d20a0467f2af9f9e0fe6cf87567b2c13 Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Wed, 28 Oct 2026 08:33:43 +0100 Subject: [PATCH 2/3] reuse liebe code --- src/core/physics/physics_funcs.cc | 20 +++++++++++++ src/core/physics/physics_funcs.h | 16 ++++++++++ src/core/radar/radar.cc | 23 ++------------ src/python_interface/py_physics.cpp | 35 ++++++++++++++++++++++ tests/core/scat/radar_single_scattering.py | 15 +++------- 5 files changed, 78 insertions(+), 31 deletions(-) diff --git a/src/core/physics/physics_funcs.cc b/src/core/physics/physics_funcs.cc index 5026fc634a..0707ff5775 100644 --- a/src/core/physics/physics_funcs.cc +++ b/src/core/physics/physics_funcs.cc @@ -144,6 +144,26 @@ Numeric refractive_index_water_visible_nir_harvey98(const Numeric frequency, return std::sqrt((1.0 + 2.0 * lorentz_lorenz) / (1.0 - lorentz_lorenz)); } +Complex permittivity_water_liebe93(const Numeric frequency, const Numeric temperature, const bool check_validity) { + if (check_validity) { + if (not(temperature >= 233.15 and temperature <= 373.15)) + arts::user_error("Liebe93 water permittivity is valid for 233.15 K <= temperature <= 373.15 K. Got {} K.", + temperature); + if (not(frequency >= 10e9 and frequency <= 1000e9)) + arts::user_error("Liebe93 water permittivity is valid for 10 GHz <= frequency <= 1000 GHz. Got {} Hz.", + frequency); + } + + const Numeric theta = 1.0 - 300.0 / temperature; + const Numeric e0 = 77.66 - 103.3 * theta; + const Numeric e1 = 0.0671 * e0; + const Numeric f1 = 20.2 + 146.0 * theta + 316.0 * theta * theta; + const Numeric e2 = 3.52; + const Numeric f2 = 39.8 * f1; + const Complex ifGHz{0.0, frequency / 1e9}; + return e2 + (e1 - e2) / (1.0 - ifGHz / f2) + (e0 - e1) / (1.0 - ifGHz / f1); +} + /** fresnel * * Calculates complex AMPLITUDE reflection coeffcients for a specular diff --git a/src/core/physics/physics_funcs.h b/src/core/physics/physics_funcs.h index 9955ae83d8..1e997f49fa 100644 --- a/src/core/physics/physics_funcs.h +++ b/src/core/physics/physics_funcs.h @@ -97,4 +97,20 @@ Numeric refractive_index_water_visible_nir_harvey98(Numeric frequency, Numeric density, bool check_validity = true); +/** Complex permittivity of liquid water in the microwave. + * + * Implements the double-Debye model of Liebe, Hufford and Manabe (1991), "A + * model for the complex permittivity of water at frequencies below 1 THz", + * as ARTS2's complex_n_water_liebe93 and Atmlab's eps_water_liebe93.m. The + * relaxation-frequency coefficient is 146, not the 146.4 of epswater93.m. + * The imaginary part is positive. + * + * @param frequency Frequency [Hz]. + * @param temperature Temperature [K]. + * @param check_validity Enforce 10 GHz <= frequency <= 1000 GHz and + * 233.15 K <= temperature <= 373.15 K. + * @return Complex relative permittivity. + */ +Complex permittivity_water_liebe93(Numeric frequency, Numeric temperature, bool check_validity = true); + #endif // physics_h diff --git a/src/core/radar/radar.cc b/src/core/radar/radar.cc index 4210bd9314..b5653d96cc 100644 --- a/src/core/radar/radar.cc +++ b/src/core/radar/radar.cc @@ -2,34 +2,17 @@ #include #include +#include #include #include -#include #include #include namespace radar { Numeric liebe93_k2(const Numeric frequency, const Numeric temperature) { - if (temperature < Constant::temperature_at_0c - 40.0 or temperature > Constant::temperature_at_0c + 100.0) - arts::user_error("Automatic radar k2 requires {} K <= ze_tref <= {} K; got {} K", - Constant::temperature_at_0c - 40.0, - Constant::temperature_at_0c + 100.0, - temperature); - if (frequency < 10e9 or frequency > 1000e9) - arts::user_error("Automatic radar k2 uses Liebe-93 and requires 10 GHz <= frequency <= 1000 GHz; got {} Hz", - frequency); - - const Numeric theta = 1.0 - 300.0 / temperature; - const Numeric e0 = 77.66 - 103.3 * theta; - const Numeric e1 = 0.0671 * e0; - const Numeric f1 = 20.2 + 146.0 * theta + 316.0 * theta * theta; - const Numeric e2 = 3.52; - const Numeric f2 = 39.8 * f1; - const Complex ifGHz{0.0, frequency / 1e9}; - const Complex n = std::sqrt(e2 + (e1 - e2) / (1.0 - ifGHz / f2) + (e0 - e1) / (1.0 - ifGHz / f1)); - const Complex n2 = n * n; - return std::norm((n2 - 1.0) / (n2 + 2.0)); + const Complex permittivity = permittivity_water_liebe93(frequency, temperature); + return std::norm((permittivity - 1.0) / (permittivity + 2.0)); } Numeric ze_factor(const Numeric frequency, const Numeric ze_tref, const Numeric k2) { diff --git a/src/python_interface/py_physics.cpp b/src/python_interface/py_physics.cpp index 76aac593c5..91dcc88e89 100644 --- a/src/python_interface/py_physics.cpp +++ b/src/python_interface/py_physics.cpp @@ -1,4 +1,5 @@ #include +#include #include #include "hpy_numpy.h" @@ -51,6 +52,40 @@ n : Numeric or numpy.ndarray Real refractive index. )--"); + physics.def( + "permittivity_water_liebe93", + [](py::object frequency, py::object temperature, bool check_validity) { + return vectorize( + [check_validity](Numeric f, Numeric t) { return permittivity_water_liebe93(f, t, check_validity); }, + frequency, + temperature); + }, + "frequency"_a, + "temperature"_a, + "check_validity"_a = true, + R"--(Complex permittivity of liquid water in the microwave. + +Implements the double-Debye model of Liebe, Hufford and Manabe (1991), "A +model for the complex permittivity of water at frequencies below 1 THz", as +ARTS2 and Atmlab's eps_water_liebe93.m. The imaginary part is positive. The +function broadcasts scalar and array inputs using NumPy rules. + +Parameters +---------- +frequency : Numeric or numpy.ndarray + Frequency [Hz]. +temperature : Numeric or numpy.ndarray + Temperature [K]. +check_validity : bool, optional + Enforce the ranges 10 GHz <= frequency <= 1000 GHz and + 233.15 K <= temperature <= 373.15 K. + +Returns +------- +eps : complex or numpy.ndarray + Complex relative permittivity. +)--"); + physics.def( "number_density", [](py::object p, py::object t) { diff --git a/tests/core/scat/radar_single_scattering.py b/tests/core/scat/radar_single_scattering.py index e068beaffc..693137b303 100644 --- a/tests/core/scat/radar_single_scattering.py +++ b/tests/core/scat/radar_single_scattering.py @@ -219,18 +219,11 @@ def manual_atmospheric_derivative(workspace, key, flat_index, step, kwargs): grids=[fine_altitude, *old.atm_field["t"].data.grids[1:]], ) -# Same Liebe (1993) water model as Atmlab eps_water_liebe93.m. Its -# relaxation-frequency coefficient is 146, not epswater93.m's 146.4. -# Preserve the original frequency, temperatures, diameter and angular grid. +# The droplet uses the same Liebe (1993) water model as ARTS2 and as the +# automatic radar k2. Preserve the original frequency, temperatures, diameter +# and angular grid. temperatures = np.array([268.15, 273.15, 278.15]) -theta = 300.0 / temperatures - 1.0 -e0 = 77.66 + 103.3 * theta -e1 = 0.0671 * e0 -g1 = 20.2 - 146.0 * theta + 316.0 * theta**2 -permittivity = e0 - 94.0 * ( - (e0 - e1) / (94.0 + 1j * g1) - + (e1 - 3.52) / (94.0 + 1j * 39.8 * g1) -) +permittivity = A.physics.permittivity_water_liebe93(94e9, temperatures) habit = A.ParticleHabit.sphere( temperatures, [94e9], [50e-6], A.IrregularZenithAngleGrid(np.arange(181.0)), From 13e409791cda33652d6d0d06da0d7e051c724baa Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Fri, 30 Oct 2026 09:02:06 +0100 Subject: [PATCH 3/3] fix comments --- src/core/sensor/obsel.cpp | 14 +++++ src/core/sensor/obsel.h | 14 +++++ src/m_radar.cc | 22 +++++--- src/tests/test_ai_port_regressions.cc | 74 +++++++++++++++++++++++++++ src/workspace_methods.cpp | 6 ++- 5 files changed, 121 insertions(+), 9 deletions(-) diff --git a/src/core/sensor/obsel.cpp b/src/core/sensor/obsel.cpp index 0ca0534f08..d7d3b48a85 100644 --- a/src/core/sensor/obsel.cpp +++ b/src/core/sensor/obsel.cpp @@ -240,6 +240,20 @@ void Obsel::sumup(VectorView out, const StokvecMatrixView& j, Index ip) const { } } +void Obsel::geometry_weights(VectorView out, Index ip) const { + assert(out.size() == f->size()); + assert(ip < static_cast(poslos->size()) and ip >= 0); + + // w is a sparse sorted matrix of shape poslos->size() x f->size() + const SparseStokvec v0{.irow = static_cast(ip), .icol = static_cast(0)}; + const SparseStokvec vn{.irow = static_cast(ip), .icol = std::numeric_limits::max()}; + std::span span{stdr::lower_bound(w, v0), stdr::upper_bound(w, vn)}; + + out = 0.0; + for (const auto& ws : span) + out[ws.icol] = std::max(std::abs(ws.data.I()), std::hypot(ws.data.Q(), ws.data.U(), ws.data.V())); +} + Size Obsel::flat_size(const SensorKeyType& key) const { switch (key) { using enum SensorKeyType; diff --git a/src/core/sensor/obsel.h b/src/core/sensor/obsel.h index 20a488d427..de2ee3ef96 100644 --- a/src/core/sensor/obsel.h +++ b/src/core/sensor/obsel.h @@ -160,6 +160,20 @@ class Obsel { [[nodiscard]] Numeric sumup(const StokvecVectorView& i, Index ip) const; void sumup(VectorView out, const StokvecMatrixView& j, Index ip) const; + /** The geometric weight of each frequency at one poslos + * + * The weight of a polarization selector is the smallest total weight of + * physical receivers (|pol| <= I) whose signed sum gives it, that is + * max(|I|, sqrt(Q^2 + U^2 + V^2)). This equals I for physical selectors, + * but is also positive for the pure Q, U, and V difference selectors, whose + * I weight is zero. Use it to average scalar quantities over the sensor + * element independently of its polarization. + * + * @param out The weights [f_grid().size()] + * @param ip The poslos index + */ + void geometry_weights(VectorView out, Index ip) const; + [[nodiscard]] Size flat_size(const SensorKeyType& key) const; void flat(VectorView x, const SensorKeyType& key) const; [[nodiscard]] Vector flat(const SensorKeyType& key) const; diff --git a/src/m_radar.cc b/src/m_radar.cc index fc18ecdd20..df67b6923c 100644 --- a/src/m_radar.cc +++ b/src/m_radar.cc @@ -119,6 +119,7 @@ RadarResult radar_forward(const Workspace& ws, .time = Vector(ny, 0.0)}; std::vector covered(ny, false); Vector scalar_weight_sum(ny, 0.0); + Vector time_coverage(ny, 0.0); for (const auto& simulation : simulations) { if (simulation.samples.empty()) continue; @@ -188,6 +189,7 @@ RadarResult radar_forward(const Workspace& ws, particle_extinction[joker, iv] = radar::gate_sum(samples, optics.particle_extinction[iv], ngate); StokvecVector bin(nf); + Vector geometry(nf); for (Size ig = 0; ig < ngate; ++ig) { if (coverage[ig] == 0.0) continue; @@ -198,10 +200,15 @@ RadarResult radar_forward(const Workspace& ws, for (Size iv = 0; iv < nf; ++iv) bin[iv] = factor[iv] * signal[ig, iv]; out.measurement[iy] += obsel.sumup(bin, iposlos); - // Extinction, background and time are means over the covered part of the gate - const Numeric scalar_weight = obsel.sumup(Stokvec{1.0, 0.0, 0.0, 0.0}, iposlos); + // Extinction, background and time are means over the covered part of the + // gate. They are weighted by the geometry, not the polarization, of the + // sensor element. The time is also weighted by the coverage of every + // path, as in MCRadar. + obsel.geometry_weights(geometry, iposlos); + const Numeric scalar_weight = sum(geometry); scalar_weight_sum[iy] += scalar_weight; - out.time[iy] += scalar_weight * time[ig] / coverage[ig]; + out.time[iy] += scalar_weight * time[ig]; + time_coverage[iy] += scalar_weight * coverage[ig]; for (Size ia = 0; ia < aux_vars.size(); ++ia) { switch (aux_vars[ia]) { case RadarAuxiliary::RadiativeBackground: out.auxiliary[ia, iy] += scalar_weight * background; break; @@ -210,12 +217,10 @@ RadarResult radar_forward(const Workspace& ws, out.auxiliary[ia, iy] += obsel.sumup(bin, iposlos); break; case RadarAuxiliary::AbsSpeciesExtinction: - for (Size iv = 0; iv < nf; ++iv) bin[iv] = Stokvec{gas_extinction[ig, iv] / coverage[ig]}; - out.auxiliary[ia, iy] += obsel.sumup(bin, iposlos); + out.auxiliary[ia, iy] += dot(geometry, gas_extinction[ig]) / coverage[ig]; break; case RadarAuxiliary::ParticleExtinction: - for (Size iv = 0; iv < nf; ++iv) bin[iv] = Stokvec{particle_extinction[ig, iv] / coverage[ig]}; - out.auxiliary[ia, iy] += obsel.sumup(bin, iposlos); + out.auxiliary[ia, iy] += dot(geometry, particle_extinction[ig]) / coverage[ig]; break; } } @@ -246,8 +251,9 @@ RadarResult radar_forward(const Workspace& ws, continue; } + out.time[iy] = time_coverage[iy] == 0.0 ? nan : out.time[iy] / time_coverage[iy]; + if (scalar_weight_sum[iy] == 0.0) continue; - out.time[iy] /= scalar_weight_sum[iy]; for (Size ia = 0; ia < aux_vars.size(); ++ia) if (aux_vars[ia] != RadarAuxiliary::Backscattering) out.auxiliary[ia, iy] /= scalar_weight_sum[iy]; } diff --git a/src/tests/test_ai_port_regressions.cc b/src/tests/test_ai_port_regressions.cc index b4c554db9c..0a81ac0697 100644 --- a/src/tests/test_ai_port_regressions.cc +++ b/src/tests/test_ai_port_regressions.cc @@ -256,6 +256,80 @@ void radar_gates() { for (Numeric x : error[i]) close(x, 0.0, "Pencil beam is deterministic"); close(time[i] * Constant::speed_of_light / (2.0 * midpoint), 1.0, "Radar time must be the round trip to the gate"); } + + // One observation from two radar heights. The lower path ends halfway + // through the gate at the surface, so the gate time must be weighted by + // the covered part of each path, as in MCRadar. + const Agenda observer_path = [] { + Agenda agenda{"ray_path_observer_agenda"}; + CallbackOperator cb; + cb.inputs = {"obs_pos", "obs_los"}; + cb.outputs = {"ray_path"}; + cb.callback = CallbackOperator::func_t{[](Workspace& local) { + const Numeric start = local.get("obs_pos")[0]; + auto& path = local.get("ray_path"); + path.resize(11); + for (Size i = 0; i < path.size(); ++i) { + path[i].pos = Vector3{start * (1.0 - static_cast(i) / 10.0), 0.0, 0.0}; + path[i].los = Vector2{0.0, 0.0}; + path[i].pos_type = PathPositionType::atm; + path[i].los_type = PathPositionType::atm; + } + }}; + agenda.add(Method{"test observer path", Wsv{cb}}); + agenda.finalize(true); + return agenda; + }(); + Matrix limits(1, 2); + limits[0, 0] = 950.0; + limits[0, 1] = 1050.0; + Vector y, radar_time; + Matrix jac, aux; + const auto simulate = [&](const Stokvec& selector) { + sensor::SparseStokvecMatrix weights(2, 1); + weights[0, 0] = selector; + weights[1, 0] = selector; + const ArrayOfSensorObsel sensor{SensorObsel(AscendingGrid{94e9}, + SensorPosLosVector{{.pos = {1000.0, 0.0, 0.0}, .los = {180.0, 0.0}}, + {.pos = {1050.0, 0.0, 0.0}, .los = {180.0, 0.0}}}, + weights)}; + measurement_vecFromRadarSingleScattering(ws, + y, + jac, + aux, + radar_time, + sensor, + limits, + JacobianTargets{}, + atmosphere, + surface, + species, + observer_path, + gas, + TransmittanceOption::constant, + Stokvec{1.0, 1.0, 0.0, 0.0}, + RadarRangeMode::Distance, + RadarUnit::unit, + 273.15, + 0.93, + -99.0, + 1.0, + ArrayOfRadarAuxiliary{RadarAuxiliary::ParticleExtinction}); + }; + simulate(Stokvec{0.5, 0.5, 0.0, 0.0}); + // Half of the gate at a mean of 975 m, and all of it at a mean of 1000 m + const Numeric distance = (0.5 * 975.0 + 1.0 * 1000.0) / 1.5; + close(radar_time[0] * Constant::speed_of_light / (2.0 * distance), + 1.0, + "Radar time must be weighted by the coverage of every path"); + close(aux[0, 0] / extinction, 1.0, "Particle extinction must be the mean over the gate"); + + // A pure Q selector has no I weight but the same geometry. The time and + // the scalar auxiliary variables must not depend on the polarization. + const Numeric time_iv = radar_time[0]; + simulate(Stokvec{0.0, 1.0, 0.0, 0.0}); + close(radar_time[0] / time_iv, 1.0, "Radar time must not depend on the polarization selector"); + close(aux[0, 0] / extinction, 1.0, "Particle extinction must not depend on the polarization selector"); } } // namespace diff --git a/src/workspace_methods.cpp b/src/workspace_methods.cpp index 4150b2ccf3..21c54d0044 100644 --- a/src/workspace_methods.cpp +++ b/src/workspace_methods.cpp @@ -6986,7 +6986,11 @@ analytically through rtepack. There is no production finite-difference loop. The auxiliary quantities are listed by *RadarAuxiliary*. They are returned as rows of *radar_aux* in the requested order. Extinction and background are -means over the covered part of each gate. +means over the covered part of each gate. They and *radar_time* are weighted +by the geometry of the sensor elements, not by their polarization. The weight +of a polarization selector is the larger of the absolute I and the length of +(Q, U, V). This is I for a physical receiver, and positive also for pure Q, U, +and V selectors. )--", .author = {"Patrick Eriksson", "OpenAI Codex"}, .out = {"measurement_vec", "measurement_jac", "radar_aux", "radar_time"},