rustmatrix.spectra#

Doppler and polarimetric spectra.

Evaluate Doppler-velocity-resolved scattering matrices S_spec(v) and Z_spec(v) for a hydrometeor PSD (single species) or a HydroMix, then derive the spectral radar observables (sZ_h, sZ_dr, sK_dp, sρ_hv, sδ_hv, sLDR) by applying the bulk rustmatrix.radar formulas bin-by-bin.

Physics#

Line-of-sight velocity for a vertical-pointing radar is v = v_t(D) + w, where v_t(D) is terminal fall speed and w is the mean vertical air motion.

Sign convention: positive velocity = toward a downward-pointing radar = fall direction. For an up-looking radar flip the sign of w and v_bins (it is a pure sign flip — the kernels are even).

Finite beamwidth combined with horizontal wind broadens the spectrum by a deterministic Gaussian of width σ_beam = |u_h| · θ_b / (2 √(2 ln 2)) (Doviak & Zrnić 1993, small-θ_b vertical-pointing limit). This σ_beam adds in quadrature to the diameter-dependent turbulence width σ_t(D) when the spectral Gaussian kernel is built.

S_spec and Z_spec are linear in N(D), so a HydroMix is handled by summing per-component spectra on a shared velocity grid.

System noise#

Radar receivers have a thermal noise floor that appears in every velocity bin. The noise is uncorrelated between H and V channels, so it adds incoherently to sZ_h and sZ_v but does not feed the underlying S_spec / Z_spec matrices — those remain the signal-only scattering matrices and continue to round-trip to the bulk radar.* helpers exactly. Noise biases sZ_dr, sρ_hv, sLDR through the per-bin SNR.

Pass noise= to SpectralIntegrator:

  • None (default) — no noise, the spectrum is signal-only.

  • "realistic" or True — use realistic_noise_floor() for the scatterer’s wavelength (a -20 dBZ total reflectivity by default).

  • a float — total noise reflectivity [mm⁶ m⁻³] distributed uniformly across the velocity grid and applied equally to H and V.

  • a 2-tuple (noise_h, noise_v) [mm⁶ m⁻³] — separate H and V noise.

References

Doviak, R. J. & Zrnić, D. S. (1993). Doppler Radar and Weather Observations, 2nd ed., Academic Press. (§5.3 beam-broadening.)

Zeng, Y., Janjić, T., de Lozar, A., Blahak, U., Reich, H., Keil, C., Seifert, A., Hagen, M., Tridon, F., Kneifel, S. (2023). “Interpreting Doppler spectra measured by cloud radars using a particle inertia model.” Atmos. Meas. Tech., 16, 3727–3757. doi:10.5194/amt-16-3727-2023

class rustmatrix.spectra.fall_speed[source]#

Bases: object

Namespace of fall-speed presets. Call any as fall_speed.preset(D).

static atlas_srivastava_sekhon_1973(rho_ratio=1.0)#

Atlas–Srivastava–Sekhon (1973) rain terminal velocity [m/s].

v_t(D) = (9.65 − 10.3 · exp(−0.6 · D)) · rho_ratio

Parameters:

  • D (array_like) – Equivalent drop diameter [mm].

  • rho_ratio (float) – (ρ_0 / ρ_air)^0.4 air-density correction (1.0 at sea level).

static brandes_et_al_2002(D)#

Brandes et al. (2002) 4th-order polynomial rain fall speed [m/s].

v_t(D) = −0.1021 + 4.932 D − 0.9551 D² + 0.07934 D³ − 0.002362 D⁴ (valid 0.1 ≤ D ≤ 8 mm; clamped to zero below 0.1 mm).

static beard_1976(T=293.15, P=101325.0)#

Beard-style rain terminal velocity with explicit T, P dependence [m/s].

Uses the Brandes et al. (2002) sea-level fit as the reference fall speed and applies Beard’s (1977) density correction (ρ₀ / ρ)^0.5 to shift to arbitrary (T, P). This is the practical form used in most cloud microphysics codes when “Beard 1976” is requested with a non-standard ambient state — the explicit three-regime Beard (1976) formulation is notoriously sensitive to coefficient transcription and gains little accuracy over the much-simpler Brandes fit below 7 mm.

Parameters:

  • D (array_like) – Equivalent drop diameter [mm].

  • T (float) – Air temperature [K]. Default 293.15 K (20 °C).

  • P (float) – Air pressure [Pa]. Default 101325 Pa.

References

Beard, K. V. (1977). Terminal velocity adjustment for cloud and precipitation drops aloft. J. Atmos. Sci., 34, 1293–1298. Brandes, E. A., Zhang, G., & Vivekanandan, J. (2002). Experiments in rainfall estimation with a polarimetric radar in a subtropical environment. J. Appl. Meteor., 41, 674–685.

static locatelli_hobbs_1974_aggregates(D)#

Locatelli & Hobbs (1974) unrimed aggregates (mixed habits) [m/s].

v_t(D) = 0.69 · D^0.41 (D in mm, valid 1 ≤ D ≤ 10 mm).

static locatelli_hobbs_1974_graupel_hex(D)#

Locatelli & Hobbs (1974) hexagonal graupel [m/s].

v_t(D) = 1.1 · D^0.57 (D in mm, valid 0.5 ≤ D ≤ 2 mm).

static power_law(b, D_ref=1.0, c=0.0)#

Return a user-parametrised power-law fall-speed callable.

v_t(D) = a · (D / D_ref)^b + c [m/s, D in mm]

Parameters:

  • a (float) – Prefactor and exponent of the power law.

  • b (float) – Prefactor and exponent of the power law.

  • D_ref (float) – Reference diameter used to normalise D (mm). Default 1.0.

  • c (float) – Additive offset [m/s]. Default 0.0.

Return type:

Callable

class rustmatrix.spectra.turbulence[source]#

Bases: object

Namespace for turbulence models. turbulence.from_params(...) picks the right concrete model from simple scalar inputs.

class NoTurbulence#

Bases: _TurbulenceModel

σ_t(D) = 0 everywhere (delta-function binning).

class GaussianTurbulence(sigma_t)#

Bases: _TurbulenceModel

Constant Gaussian broadening: σ_t(D) = σ_t.

The canonical spectral-polarimetry assumption. Matches classical papers that parameterise broadening with a single bulk σ_t.

Parameters:

sigma_t (float)

class InertialZeng2023(sigma_air, epsilon, L_o=100.0, v_t_ref=None)#

Bases: _TurbulenceModel

Particle-inertia-corrected turbulence response (Zeng et al. 2023).

Heavy drops cannot follow small-scale eddies, so their effective turbulence broadening is a low-pass-filtered version of the ambient air turbulence. We use a first-order response:

σ_t(D) = σ_air / √(1 + St(D)²)

with Stokes number St(D) = τ_p(D) / τ_eddy, particle relaxation time τ_p(D) = v_t_ref(D) / g (g = 9.80665 m/s²), and eddy turnover time τ_eddy = (L_o² / ε)^(1/3).

Limits: * Small drops (St → 0) → GaussianTurbulence(σ_air). * Large drops (St → ∞) → NoTurbulence().

Parameters:

  • sigma_air (float) – Ambient turbulence intensity [m/s].

  • epsilon (float) – Turbulent kinetic energy dissipation rate [m²/s³].

  • L_o (float, optional) – Outer (integral) length scale [m]. Default 100 m.

  • v_t_ref (callable, optional) – D -> v_t(D) used to estimate τ_p. Defaults to atlas_srivastava_sekhon_1973(). Pass the same fall-speed model used in the integrator for consistency.

References

Zeng et al. (2023), Atmos. Meas. Tech., 16, 3727.

static from_params(sigma=None, epsilon=None, L_o=100.0, v_t_ref=None)[source]#

Smart constructor.

Return type:

_TurbulenceModel

Parameters:
class rustmatrix.spectra.SpectralIntegrator(source, fall_speed=None, turbulence=None, *, component_kinematics=None, v_bins=None, v_min=None, v_max=None, n_bins=None, w=0.0, u_h=0.0, beamwidth=0.0, geometry_backscatter=(0.0, 180.0, 0.0, 0.0, 0.0, 0.0), geometry_forward=None, noise=None)[source]#

Bases: object

Build spectral S and Z matrices on a Doppler-velocity grid.

Parameters:

  • source (Scatterer or HydroMix) – The scattering target. A Scatterer must have an initialised psd_integrator and a psd attached. A HydroMix carries its own components (each with scatterer.psd).

  • fall_speed (callable, optional) – D_mm -> v_t_m_per_s for the single-species case. Ignored when source is a HydroMix.

  • turbulence (_TurbulenceModel or callable or None, optional) – Turbulence model for the single-species case. Ignored when source is a HydroMix.

  • component_kinematics (mapping, optional) – For HydroMix inputs only. Maps component label (or index) to (fall_speed, turbulence). Every component in the mixture must be represented.

  • v_bins (array_like, optional) – Explicit 1-D velocity grid [m/s]. Mutually exclusive with (v_min, v_max, n_bins).

  • v_min (float, float, int, optional) – Convenience triple that builds v_bins = linspace(v_min, v_max, n_bins).

  • v_max (float, float, int, optional) – Convenience triple that builds v_bins = linspace(v_min, v_max, n_bins).

  • n_bins (float, float, int, optional) – Convenience triple that builds v_bins = linspace(v_min, v_max, n_bins).

  • w (float) – Mean vertical air motion [m/s], positive-downward. Default 0.

  • u_h (float) – Horizontal wind speed magnitude [m/s]. Drives beam broadening. Default 0.

  • beamwidth (float) – One-way half-power beamwidth θ_b [rad]. Default 0 (pencil beam).

  • geometry_backscatter (6-tuple) – Backscatter geometry. Default geom_vert_back.

  • geometry_forward (6-tuple, optional) – Forward-scatter geometry. Required if you want sK_dp.

  • noise (None / "realistic" / True / float / (float, float), optional) – System noise floor. None (default) disables noise and the spectrum is signal-only — SpectralResult.collapse_to_bulk() then round-trips to the bulk radar observables exactly. Pass "realistic" or True to use realistic_noise_floor() for the scatterer’s wavelength; a scalar for equal H/V noise in mm⁶ m⁻³; or a 2-tuple (noise_h, noise_v). Noise is distributed uniformly across v_bins, added to sZ_h / sZ_v, and biases sZ_dr / sρ_hv / sLDR through per-bin SNR. The underlying S_spec and Z_spec stay signal-only.

run()[source]#

Evaluate the spectrum and derived observables.

Return type:

SpectralResult

class rustmatrix.spectra.SpectralResult(v, S_spec, Z_spec, sZ_h, sZ_v, sZ_dr, sKdp, srho_hv, sdelta_hv, sLDR, wavelength, Kw_sqr, noise_h=0.0, noise_v=0.0, _S_spec_forward=None)[source]#

Bases: object

Doppler-resolved scattering matrices and derived observables.

Parameters:
v#

Doppler-velocity bin centres [m/s].

Type:

ndarray (M,)

S_spec#

Spectral amplitude matrix.

Type:

ndarray (2, 2, M) complex

Z_spec#

Spectral phase matrix.

Type:

ndarray (4, 4, M) float

sZ_h, sZ_v

Spectral linear reflectivity [mm⁶ m⁻³ / (m/s)], horizontal / vertical polarisation.

Type:

ndarray (M,)

sZ_dr#

Spectral differential reflectivity (linear ratio).

Type:

ndarray (M,)

sKdp#

Spectral specific differential phase [° / km]. None when no forward-scatter geometry was supplied.

Type:

ndarray (M,) or None

srho_hv#

Spectral copolar correlation coefficient (0–1).

Type:

ndarray (M,)

sdelta_hv#

Spectral backscatter differential phase [rad].

Type:

ndarray (M,)

sLDR#

Spectral linear depolarisation ratio (linear).

Type:

ndarray (M,)

wavelength#

Wavelength [mm] (copied from the scatterer for collapse_to_bulk).

Type:

float

Kw_sqr#

|K_w|² (copied for collapse_to_bulk).

Type:

float

collapse_to_bulk()[source]#

Integrate S_spec / Z_spec over v and return a Scatterer-shaped shim so bulk radar helpers can be called on the integrated matrices.

The returned object has .get_S(), .get_Z(), .wavelength, .Kw_sqr, and a backscatter geometry — enough for radar.refl, radar.Zdr, radar.rho_hv, radar.delta_hv to work. If the spectrum was built with a forward geometry, .get_S_forward() returns the forward-summed S so you can hand-construct radar.Kdp-style quantities.

Return type:

SimpleNamespace

class rustmatrix.spectra.NoTurbulence[source]#

Bases: _TurbulenceModel

σ_t(D) = 0 everywhere (delta-function binning).

class rustmatrix.spectra.GaussianTurbulence(sigma_t)[source]#

Bases: _TurbulenceModel

Constant Gaussian broadening: σ_t(D) = σ_t.

The canonical spectral-polarimetry assumption. Matches classical papers that parameterise broadening with a single bulk σ_t.

Parameters:

sigma_t (float)

class rustmatrix.spectra.InertialZeng2023(sigma_air, epsilon, L_o=100.0, v_t_ref=None)[source]#

Bases: _TurbulenceModel

Particle-inertia-corrected turbulence response (Zeng et al. 2023).

Heavy drops cannot follow small-scale eddies, so their effective turbulence broadening is a low-pass-filtered version of the ambient air turbulence. We use a first-order response:

σ_t(D) = σ_air / √(1 + St(D)²)

with Stokes number St(D) = τ_p(D) / τ_eddy, particle relaxation time τ_p(D) = v_t_ref(D) / g (g = 9.80665 m/s²), and eddy turnover time τ_eddy = (L_o² / ε)^(1/3).

Limits: * Small drops (St → 0) → GaussianTurbulence(σ_air). * Large drops (St → ∞) → NoTurbulence().

Parameters:

  • sigma_air (float) – Ambient turbulence intensity [m/s].

  • epsilon (float) – Turbulent kinetic energy dissipation rate [m²/s³].

  • L_o (float, optional) – Outer (integral) length scale [m]. Default 100 m.

  • v_t_ref (callable, optional) – D -> v_t(D) used to estimate τ_p. Defaults to atlas_srivastava_sekhon_1973(). Pass the same fall-speed model used in the integrator for consistency.

References

Zeng et al. (2023), Atmos. Meas. Tech., 16, 3727.

rustmatrix.spectra.atlas_srivastava_sekhon_1973(D, rho_ratio=1.0)[source]#

Atlas–Srivastava–Sekhon (1973) rain terminal velocity [m/s].

v_t(D) = (9.65 − 10.3 · exp(−0.6 · D)) · rho_ratio

Parameters:

  • D (array_like) – Equivalent drop diameter [mm].

  • rho_ratio (float) – (ρ_0 / ρ_air)^0.4 air-density correction (1.0 at sea level).

rustmatrix.spectra.brandes_et_al_2002(D)[source]#

Brandes et al. (2002) 4th-order polynomial rain fall speed [m/s].

v_t(D) = −0.1021 + 4.932 D − 0.9551 D² + 0.07934 D³ − 0.002362 D⁴ (valid 0.1 ≤ D ≤ 8 mm; clamped to zero below 0.1 mm).

rustmatrix.spectra.beard_1976(D, T=293.15, P=101325.0)[source]#

Beard-style rain terminal velocity with explicit T, P dependence [m/s].

Uses the Brandes et al. (2002) sea-level fit as the reference fall speed and applies Beard’s (1977) density correction (ρ₀ / ρ)^0.5 to shift to arbitrary (T, P). This is the practical form used in most cloud microphysics codes when “Beard 1976” is requested with a non-standard ambient state — the explicit three-regime Beard (1976) formulation is notoriously sensitive to coefficient transcription and gains little accuracy over the much-simpler Brandes fit below 7 mm.

Parameters:

  • D (array_like) – Equivalent drop diameter [mm].

  • T (float) – Air temperature [K]. Default 293.15 K (20 °C).

  • P (float) – Air pressure [Pa]. Default 101325 Pa.

References

Beard, K. V. (1977). Terminal velocity adjustment for cloud and precipitation drops aloft. J. Atmos. Sci., 34, 1293–1298. Brandes, E. A., Zhang, G., & Vivekanandan, J. (2002). Experiments in rainfall estimation with a polarimetric radar in a subtropical environment. J. Appl. Meteor., 41, 674–685.

rustmatrix.spectra.locatelli_hobbs_1974_aggregates(D)[source]#

Locatelli & Hobbs (1974) unrimed aggregates (mixed habits) [m/s].

v_t(D) = 0.69 · D^0.41 (D in mm, valid 1 ≤ D ≤ 10 mm).

rustmatrix.spectra.locatelli_hobbs_1974_graupel_hex(D)[source]#

Locatelli & Hobbs (1974) hexagonal graupel [m/s].

v_t(D) = 1.1 · D^0.57 (D in mm, valid 0.5 ≤ D ≤ 2 mm).

rustmatrix.spectra.power_law(a, b, D_ref=1.0, c=0.0)[source]#

Return a user-parametrised power-law fall-speed callable.

v_t(D) = a · (D / D_ref)^b + c [m/s, D in mm]

Parameters:

  • a (float) – Prefactor and exponent of the power law.

  • b (float) – Prefactor and exponent of the power law.

  • D_ref (float) – Reference diameter used to normalise D (mm). Default 1.0.

  • c (float) – Additive offset [m/s]. Default 0.0.

Return type:

Callable

rustmatrix.spectra.realistic_noise_floor(wavelength_mm=33.3, Z_dBZ=-20.0)[source]#

Return a realistic total noise reflectivity in mm⁶ m⁻³.

Converts a noise-equivalent reflectivity in dBZ to linear units. The default (-20 dBZ) is a realistic per-range-gate noise floor for a research radar; operational S-band systems are typically quieter (-30 to -40 dBZ), while cloud radars at Ka/W with short integration times can be noisier.

The wavelength argument is accepted for signature-compatibility so callers can write realistic_noise_floor(scatterer.wavelength) even though the default is currently wavelength-independent.

Parameters:

  • wavelength_mm (float) – Radar wavelength [mm]. Currently unused by the default model but reserved for future band-dependent presets.

  • Z_dBZ (float) – Total noise reflectivity [dBZ] across the Nyquist range.

Returns:

Total noise reflectivity in linear mm⁶ m⁻³.

Return type:

float

rustmatrix.spectra.REALISTIC_NOISE_DBZ = -20.0#

Default “realistic” total noise reflectivity [dBZ] when the user asks for a realistic default. Picked to be a sensible mid-range value that is representative of research radars without being wavelength-specific (operational S-band is quieter, cloud radars noisier in absolute dBZ but equivalent per-range-gate).

See also rustmatrix.spectra.beam for the explicit beam-pattern × scene integrator.