Merge branch 'main' into he/refactor-2M-tests

Author: haakon-e

docs/src/API.md

@@ -71,6 +71,16 @@ Microphysics2M.rain_evaporation
 Microphysics2M.conv_q_liq_to_q_rai
 ```
 
+## Distribution tools for 2-moment microphysics
+
+```@docs
+DistributionTools
+DistributionTools.generalized_gamma_quantile
+DistributionTools.generalized_gamma_cdf
+DistributionTools.exponential_cdf
+DistributionTools.exponential_quantile
+```
+
 # P3 scheme
 
 ```@docs

src/CloudMicrophysics.jl

@@ -11,6 +11,7 @@ function conv_q_liq_to_q_rai end
 function accretion end
 
 include("Common.jl")
+include("DistributionTools.jl")
 include("Microphysics0M.jl")
 include("Microphysics1M.jl")
 include("Microphysics2M.jl")

src/DistributionTools.jl

@@ -0,0 +1,116 @@
+"""
+    DistributionTools
+
+A module containing tools for working with size distributions.
+
+Currently, it contains tools for working with the generalized gamma distribution
+and the exponential distribution.
+
+Mostly used for the 2-moment microphysics.
+"""
+module DistributionTools
+
+import SpecialFunctions as SF
+
+"""
+    generalized_gamma_quantile(ν, μ, B, Y)
+
+Calculate the quantile (inverse cumulative distribution function) for a 
+    generalized gamma distribution parameterized in the form:
+
+    g(x) = A ⋅ x^ν ⋅ exp(-B ⋅ x^μ)
+
+# Arguments
+- `ν, μ, B`: The PDF parameters
+- `Y`: The probability level (0 ≤ Y ≤ 1) for which to compute the quantile
+
+# Returns
+- `x`: The value x such that P(X ≤ x) = Y
+"""
+function generalized_gamma_quantile(ν, μ, B, Y)
+    # Compute the inverse of the regularized incomplete gamma function
+    z = SF.gamma_inc_inv((ν + 1) / μ, Y, 1 - Y)
+    return (z / B)^(1 / μ)
+end
+
+"""
+    generalized_gamma_cdf(ν, μ, B, x)
+
+Calculate the cumulative distribution function (CDF) for a generalized gamma distribution
+parameterized in the form:
+
+    g(x) = A * x^ν * exp(-B * x^μ)
+
+The CDF gives the probability P(X ≤ x).
+
+# Arguments
+ - `ν, μ, B`: The PDF parameters
+ - `x`: The value at which to evaluate the CDF
+
+# Returns
+ - `p`: The probability P(X ≤ x)
+"""
+function generalized_gamma_cdf(ν, μ, B, x)
+    # Check input validity
+    μ > 0 || throw(DomainError(μ, "Parameter μ must be positive"))
+    B > 0 || throw(DomainError(B, "Parameter B must be positive"))
+    # Handle edge cases
+    x ≤ 0 && return zero(x)
+
+    # Compute the regularized incomplete gamma function
+    p, _ = SF.gamma_inc((ν + 1) / μ, (B * x)^μ)
+    return p
+end
+
+"""
+   exponential_cdf(D_mean, D)
+
+Calculate the cumulative distribution function (CDF) for an exponential distribution
+parameterized in the form:
+
+   n(D) = N₀ * exp(-D / D_mean)
+
+where N₀ is a normalizing constant such that the total probability is 1.
+
+# Arguments
+- `D_mean`: The mean value of the distribution
+- `D`: The point at which to evaluate the CDF (must be ≥ 0)
+
+# Returns
+- `p`: The probability P(X ≤ D)
+"""
+function exponential_cdf(D_mean, D)
+    # Check input validity
+    D_mean > 0 || throw(DomainError(D_mean, "Mean parameter must be positive"))
+    # Handle edge cases
+    D < 0 && return zero(D)
+    # Calculate CDF: P(X ≤ D) = 1 - exp(-D/D_mean)
+    return 1 - exp(-D / D_mean)
+end
+
+"""
+    exponential_quantile(D_mean, Y)
+
+Calculate the quantile (inverse cumulative distribution function) for an exponential distribution
+parameterized in the form:
+
+    n(D) = N₀ * exp(-D / D_mean)
+
+where N₀ is a normalizing constant such that the total probability is 1.
+
+# Arguments
+- `Y`: The probability level (0 ≤ Y ≤ 1) for which to compute the quantile
+- `D_mean`: The mean value of the distribution
+
+# Returns
+- `D`: The value D such that P(X ≤ D) = Y
+"""
+function exponential_quantile(D_mean, Y)
+    # Check input validity
+    (0 ≤ Y ≤ 1) || throw(DomainError(Y, "Probability Y must be in [0,1]"))
+    D_mean > 0 || throw(DomainError(D_mean, "Mean parameter must be positive"))
+    # Calculate quantile: x = -D_mean * ln(1-Y)
+    return -D_mean * log(1 - Y)
+end
+
+end # module DistributionTools

test/DistributionTools_tests.jl

@@ -0,0 +1,49 @@
+using Test
+import CloudMicrophysics.DistributionTools as DT
+
+@testset "DistributionTools" begin
+    @testset "Generalized Gamma Distribution" begin
+        # Test parameters
+        ν = 2.0
+        μ = 3.0
+        B = 1.0
+
+        # Test CDF and quantile correspondence
+        for Y in [0.1, 0.25, 0.5, 0.75, 0.9]
+            x = DT.generalized_gamma_quantile(ν, μ, B, Y)
+            p = DT.generalized_gamma_cdf(ν, μ, B, x)
+            @test p ≈ Y rtol = 1e-10
+        end
+
+        # Test edge cases
+        @test DT.generalized_gamma_cdf(ν, μ, B, 0) == 0
+        @test DT.generalized_gamma_cdf(ν, μ, B, -1) == 0
+
+        # Test error handling
+        @test_throws DomainError DT.generalized_gamma_cdf(ν, -1, B, 1)  # μ < 0
+        @test_throws DomainError DT.generalized_gamma_cdf(ν, μ, -1, 1)  # B < 0
+    end
+
+    @testset "Exponential Distribution" begin
+        # Test parameters
+        D_mean = 2.0
+
+        # Test CDF and quantile correspondence
+        for Y in [0.1, 0.25, 0.5, 0.75, 0.9]
+            D = DT.exponential_quantile(D_mean, Y)
+            p = DT.exponential_cdf(D_mean, D)
+            @test p ≈ Y rtol = 1e-10
+        end
+
+        # Test edge cases
+        @test DT.exponential_cdf(D_mean, 0) == 0
+        @test DT.exponential_cdf(D_mean, -1) == 0
+        @test DT.exponential_cdf(D_mean, Inf) ≈ 1 rtol = 1e-10
+
+        # Test error handling
+        @test_throws DomainError DT.exponential_cdf(-1, 1)  # D_mean < 0
+        @test_throws DomainError DT.exponential_quantile(D_mean, -0.1)  # Y < 0
+        @test_throws DomainError DT.exponential_quantile(D_mean, 1.1)   # Y > 1
+        @test_throws DomainError DT.exponential_quantile(-1, 0.5)     # D_mean < 0
+    end
+end

test/runtests.jl

@@ -18,3 +18,4 @@ include("performance_tests.jl")
 include("aerosol_activation_calibration.jl")
 include("ice_nucleation_calibration.jl")
 include("ventilation_tests.jl")
+include("DistributionTools_tests.jl")