diff --git a/src/TransitionMatrices.jl b/src/TransitionMatrices.jl index 11fc0c7..11f3eb3 100644 --- a/src/TransitionMatrices.jl +++ b/src/TransitionMatrices.jl @@ -74,7 +74,8 @@ export calc_T, calc_T_iitm, calc_S, calc_Z, calc_F, calc_ω, # Shape related exports export AbstractShape, AbstractAxisymmetricShape, AbstractNFoldShape, volume, volume_equivalent_radius, has_symmetric_plane, refractive_index, - rmin, rmax, Spheroid, Cylinder, Chebyshev, Prism + rmin, rmax, Spheroid, Cylinder, Chebyshev, Prism, + SuperSpheroid, SuperEllipsoid, SuperSpheroidRevolved # Re-exports export RotZYZ, Double64, Float128, ComplexF128, Arb, Acb diff --git a/src/shapes/index.jl b/src/shapes/index.jl index 7e93eed..ceb9a51 100644 --- a/src/shapes/index.jl +++ b/src/shapes/index.jl @@ -15,6 +15,12 @@ include("spheroid.jl") include("cylinder.jl") include("chebyshev.jl") +# Axisymmetric super-spheroid (surface of revolution) + +include("superspheroidrevolved.jl") + # NFold shapes include("prism.jl") +include("superspheroid.jl") +include("superellipsoid.jl") diff --git a/src/shapes/superellipsoid.jl b/src/shapes/superellipsoid.jl new file mode 100644 index 0000000..2efdc1a --- /dev/null +++ b/src/shapes/superellipsoid.jl @@ -0,0 +1,221 @@ +@doc raw""" +A super-ellipsoid scatterer (Bi et al. 2018 D2h symmetry). + +Defined by the implicit surface + +```math +\left[\left|\frac{x}{a}\right|^{2/e} + \left|\frac{y}{b}\right|^{2/e}\right]^{e/n} ++ \left|\frac{z}{c}\right|^{2/n} = 1 +``` + +For ``a = b`` this reduces to a 4-fold (D4h) shape; for ``a \neq b`` it is D2h (2-fold). +At ``e = n = 1`` it reduces to an ordinary ellipsoid. + +Conventions (Bi et al. 2018 Opt. Express, doi:10.1364/OE.26.001726): +- ``e``: east–west (equatorial cross-section) roundness exponent. +- ``n``: north–south (meridional profile) roundness exponent. +- ``e = n = 1``: ordinary ellipsoid. +- ``e = n = 2``: octahedron-like shape. + +Attributes: + +- `a`: the ``x`` semi-axis. +- `b`: the ``y`` semi-axis. +- `c`: the ``z`` (polar) semi-axis. +- `e`: the equatorial roundness exponent. +- `n`: the meridional roundness exponent. +- `m`: the relative complex refractive index. +""" +struct SuperEllipsoid{T, CT} <: AbstractNFoldShape{2, T, CT} + a::T + b::T + c::T + e::T + n::T + m::CT +end + +SuperEllipsoid(a, b, c, e, n, m) = SuperEllipsoid{typeof(a), typeof(m)}(a, b, c, e, n, m) + +# Helper shared across super-spheroid shapes (defined in superspheroid.jl): +# _beta_lgamma(x, y) — uses Arblib lgamma! + +@doc raw""" +Volume of a super-ellipsoid: + +```math +V = a\,b\,c\, e\,n \, B\!\left(\tfrac{e}{2},\, \tfrac{e}{2}\right) + B\!\left(\tfrac{n}{2},\, n+1\right) +``` + +where ``B(x,y) = \Gamma(x)\Gamma(y)/\Gamma(x+y)`` is the beta function. +At ``e = n = 1`` this reduces to ``\tfrac{4}{3}\pi abc`` (ellipsoid). + +Derivation: the cross-section at height ``z`` is the super-ellipse +``\{|x/a|^{2/e}+|y/b|^{2/e} \le (1-|z/c|^{2/n})^{n/e}\}`` +whose area is ``ab\,(1-|z/c|^{2/n})^n \cdot e\,B(e/2,e/2)``; +integrating gives the formula above. +""" +function volume(s::SuperEllipsoid{T}) where {T} + e = s.e + n = s.n + b1 = _beta_lgamma(e / 2, e / 2) # B(e/2, e/2) = Γ(e/2)²/Γ(e) + b2 = _beta_lgamma(n / 2, n + 1) # B(n/2, n+1) + return s.a * s.b * s.c * e * n * b1 * b2 +end + +volume_equivalent_radius(s::SuperEllipsoid) = ∛(3 * volume(s) / (4 * oftype(s.a, π))) + +has_symmetric_plane(::SuperEllipsoid) = true + +@testitem "SuperEllipsoid utility functions" begin + using TransitionMatrices: SuperEllipsoid, volume, volume_equivalent_radius, + has_symmetric_plane + + # e=n=1 reduces to ellipsoid: V = 4π/3 * a * b * c + @testset "e=n=1 reduction to ellipsoid" begin + a, b, c = 1.0, 1.5, 2.0 + s = SuperEllipsoid(a, b, c, 1.0, 1.0, complex(1.5)) + @test volume(s) ≈ 4π / 3 * a * b * c rtol=1e-8 + @test volume_equivalent_radius(s) ≈ ∛(a * b * c) rtol=1e-8 + end + + # a=b, e=n=1 reduces to spheroid + @testset "a=b, e=n=1 reduces to spheroid" begin + a, c = 1.0, 2.0 + s = SuperEllipsoid(a, a, c, 1.0, 1.0, complex(1.5)) + @test volume(s) ≈ 4π / 3 * a^2 * c rtol=1e-8 + end + + # Volume vs independent numerical integration via cross-sectional slices + # V = ∫_{-c}^{c} A(z) dz, A(z) = ab·R(z)^e·A_unit + # R(z) = (1-(|z|/c)^(2/n))^(n/e), A_unit = area({|u|^(2/e)+|v|^(2/e)≤1}) = 4·∫_0^1(1-u^(2/e))^(e/2)du + @testset "Volume vs cross-section integration for e=$e, n=$n" for (e, n) in [(1.5, 1.2), (2.0, 2.0)] + a, b, c = 1.0, 1.2, 0.8 + s = SuperEllipsoid(a, b, c, e, n, complex(1.5)) + V_analytic = volume(s) + + Npts = 5000 + du = 1.0 / Npts + A_unit = 4 * sum((1 - ((i - 0.5) * du)^(2 / e))^(e / 2) for i in 1:Npts) * du + + dz = 2c / Npts + V_num = 0.0 + for i in 1:Npts + z = -c + (i - 0.5) * dz + t = abs(z) / c + t < 1 || continue + R = (1 - t^(2 / n))^(n / e) + V_num += a * b * R^e * A_unit * dz + end + + @test abs(V_analytic - V_num) / V_analytic < 0.005 + end + + # Fully independent 3D Monte-Carlo via `∈` (shares no derivation with the closed-form + # beta-function volume — guards the corrected SuperEllipsoid formula). Deterministic + # LCG seed so the test is reproducible, not flaky. + @testset "Volume vs 3D Monte-Carlo for e=$e, n=$n" for (e, n) in [(0.8, 1.6), (1.5, 1.2), (2.0, 2.5)] + a, b, c = 1.0, 1.2, 0.8 + s = SuperEllipsoid(a, b, c, e, n, complex(1.5)) + V_analytic = volume(s) + N = 400_000 + hits = 0 + let seed = UInt64(20240607) + x = zeros(3) + for _ in 1:N + seed = seed * 6364136223846793005 + 1442695040888963407 + x[1] = (Float64(seed >> 11) / 2^53 * 2 - 1) * a + seed = seed * 6364136223846793005 + 1442695040888963407 + x[2] = (Float64(seed >> 11) / 2^53 * 2 - 1) * b + seed = seed * 6364136223846793005 + 1442695040888963407 + x[3] = (Float64(seed >> 11) / 2^53 * 2 - 1) * c + hits += x ∈ s ? 1 : 0 + end + end + V_mc = 2a * 2b * 2c * hits / N + @test abs(V_analytic - V_mc) / V_analytic < 0.01 + end + + @test has_symmetric_plane(SuperEllipsoid(1.0, 1.2, 0.8, 1.5, 1.2, complex(1.5))) +end + +function rmax(s::SuperEllipsoid) + # For e,n ≥ 1: the surface never extends beyond the cardinal semi-axes, + # so rmax = max(a, b, c). + # For e or n < 1 (cube-like bumpy shapes): corners can extend beyond axes; + # hypot gives a safe conservative upper bound in that case. + if s.e >= 1 && s.n >= 1 + return max(s.a, s.b, s.c) + else + return hypot(s.a, s.b, s.c) + end +end + +function rmin(s::SuperEllipsoid) + # Inscribed-sphere radius = min over directions û of the surface radius + # r(û) = ((|ux/a|^(2/e)+|uy/b|^(2/e))^(e/n) + |uz/c|^(2/n))^(-n/2). + # There is NO simple closed form: for concave superquadrics the closest surface point + # sits at a transcendental, off-symmetry direction (the symmetry-direction minimum + # over-estimates by a few % for e≠n). So minimise r(û) over a dense direction grid. + # An over-large rₘᵢₙ — e.g. min(a,b,c), which ignores the body-diagonal/face dimples — + # makes the IITM Mie-init sphere poke outside the particle and gives wrong cross + # sections, so a small shrink (below) keeps it safely inscribed. + e = Float64(s.e) + n = Float64(s.n) + a = Float64(s.a) + b = Float64(s.b) + c = Float64(s.c) + M = 90 + rmn = Inf + for i in 0:M, j in 0:M + ϑ = i * (π / 2) / M + φ = j * (π / 2) / M + ux = sin(ϑ) * cos(φ) + uy = sin(ϑ) * sin(φ) + uz = cos(ϑ) + Σ = (ux / a)^(2 / e) + (uy / b)^(2 / e) + g = Σ^(e / n) + (uz / c)^(2 / n) + rmn = min(rmn, g^(-n / 2)) + end + # The grid minimum can only OVER-estimate the true minimum (it may miss the exact + # closest direction between nodes). rₘᵢₙ must never exceed the true inscribed radius — + # an over-large value pokes the Mie-init sphere outside the particle — so shrink by a + # margin comfortably larger than the grid error. + return oftype(s.a, 0.99 * rmn) +end + +function Base.:∈(x, s::SuperEllipsoid) + inner = (abs(x[1]) / s.a)^(2 / s.e) + (abs(x[2]) / s.b)^(2 / s.e) + return inner^(s.e / s.n) + (abs(x[3]) / s.c)^(2 / s.n) <= 1 +end + +refractive_index(s::SuperEllipsoid, x) = x ∈ s ? s.m : one(s.m) + +@testitem "SuperEllipsoid e=n=1 (a=b) reduces to a spheroid (IITM)" begin + using TransitionMatrices: SuperEllipsoid, Spheroid, transition_matrix, IITM, + calc_Csca, calc_Cext + + m = complex(1.5) + # a=b, e=n=1 ⇒ the super-ellipsoid IS a spheroid. Compare both through the SAME IITM + # solver with a matched explicit rₘᵢₙ: tests the geometric reduction directly (to + # machine precision), robust to the rₘᵢₙ default and cross-method discretisation. + slv = IITM(5, 60, 200, 100; rₘᵢₙ = 0.6) + 𝐓sph = transition_matrix(Spheroid(1.0, 2.0, m), 2π, slv) + 𝐓se = transition_matrix(SuperEllipsoid(1.0, 1.0, 2.0, 1.0, 1.0, m), 2π, slv) + + @test abs(calc_Csca(𝐓sph) - calc_Csca(𝐓se)) < 1e-10 + @test abs(calc_Cext(𝐓sph) - calc_Cext(𝐓se)) < 1e-10 +end + +@testitem "SuperEllipsoid geometry predicates (incl. e,n<1 rmax branch)" begin + using TransitionMatrices: SuperEllipsoid + for (e, n) in ((0.7, 1.5), (1.0, 1.0), (2.0, 2.5)) # e<1 hits the hypot rmax branch + s = SuperEllipsoid(1.0, 1.2, 0.8, e, n, complex(1.5)) + @test TransitionMatrices.rmax(s) ≥ TransitionMatrices.rmin(s) > 0 + @test [0.0, 0.0, 0.0] ∈ s + @test [5.0, 0.0, 0.0] ∉ s + @test TransitionMatrices.refractive_index(s, [0.0, 0.0, 0.0]) == complex(1.5) + @test TransitionMatrices.refractive_index(s, [5.0, 0.0, 0.0]) == one(complex(1.5)) + end +end diff --git a/src/shapes/superspheroid.jl b/src/shapes/superspheroid.jl new file mode 100644 index 0000000..badd02b --- /dev/null +++ b/src/shapes/superspheroid.jl @@ -0,0 +1,223 @@ +@doc raw""" +A super-spheroid scatterer (Bi/Lin dust model, D4h symmetry). + +Defined by the implicit surface + +```math +\left|\frac{x}{a}\right|^{2/n} + \left|\frac{y}{a}\right|^{2/n} + \left|\frac{z}{c}\right|^{2/n} = 1 +``` + +where the cross-sections at constant ``z`` are super-circles (rounded squares for +``n \neq 1``), giving **4-fold** (D4h) rotational symmetry — NOT axial symmetry. +Use `SuperSpheroidRevolved` for the true surface of revolution. + +Conventions (Lin & Bi 2018 JGR, doi:10.1029/2018JD029464): +- ``n = 1``: ordinary spheroid ``(a = b)``. +- ``n = 2``: regular octahedron. +- ``1 < n < 2``: convex rounded shape. +- ``n > 2``: concave (dust best-fit region). + +Attributes: + +- `a`: the equatorial semi-axis (``x`` and ``y`` directions). +- `c`: the polar semi-axis (``z`` direction). +- `n`: the roundness exponent. +- `m`: the relative complex refractive index. +""" +struct SuperSpheroid{T, CT} <: AbstractNFoldShape{4, T, CT} + a::T + c::T + n::T + m::CT +end + +SuperSpheroid(a, c, n, m) = SuperSpheroid{typeof(a), typeof(m)}(a, c, n, m) + +# Beta function B(x, y) = Γ(x)Γ(y)/Γ(x+y) for the superquadric volume formulas. +# Computed in Arb (lgamma) at the inputs' own precision and converted straight back to +# `T` — the same `T(Arblib...(Arb(...); prec))` idiom the package uses elsewhere (see +# `special_functions/factorial.jl`). No BigFloat intermediate: that would both be a +# needless detour and, because `BigFloat(::Arb)` uses the *global* precision, silently +# truncate high-precision `T` (Arb, large BigFloat). `precision(x)` reads the instance's +# precision, so Float64/Double64/Float128/BigFloat/Arb shapes each compute at their own. +# Working precision (bits) from the value's own precision, so high-precision shapes +# (Double64, BigFloat, Arb instances) are not truncated. `precision(::Float128)` (the +# instance method) is unimplemented in the current Quadmath/Julia (it routes to a +# missing `_precision_with_base_2`), so Float128 is pinned to its 113-bit mantissa via +# the working type method `precision(Float128)`. +_beta_prec(::Float128) = precision(Float128) +_beta_prec(x::Real) = precision(x) + +function _beta_lgamma(x::T, y::T) where {T <: Real} + prec = max(_beta_prec(x), _beta_prec(y)) + lΓx = Arblib.lgamma!(Arb(; prec), Arb(x; prec); prec) + lΓy = Arblib.lgamma!(Arb(; prec), Arb(y; prec); prec) + lΓxy = Arblib.lgamma!(Arb(; prec), Arb(x + y; prec); prec) + return T(exp(lΓx + lΓy - lΓxy)) +end + +@doc raw""" +Volume of a super-spheroid: + +```math +V = 4\,a^2 c\, n^2 \, B\!\left(\tfrac{n}{2}+1,\, n\right) B\!\left(\tfrac{n}{2},\, \tfrac{n+2}{2}\right) +``` + +where ``B`` is the beta function. At ``n=1`` this reduces to ``\tfrac{4}{3}\pi a^2 c`` +(spheroid). +""" +function volume(s::SuperSpheroid{T}) where {T} + n = s.n + b1 = _beta_lgamma(n / 2 + 1, n) + b2 = _beta_lgamma(n / 2, (n + 2) / 2) + return 4 * s.a^2 * s.c * n^2 * b1 * b2 +end + +volume_equivalent_radius(s::SuperSpheroid) = ∛(3 * volume(s) / (4 * oftype(s.a, π))) + +has_symmetric_plane(::SuperSpheroid) = true + +@testitem "SuperSpheroid utility functions" begin + using TransitionMatrices: SuperSpheroid, volume, volume_equivalent_radius, + has_symmetric_plane, Double64, Float128, ComplexF128 + + # n=1 reduces to spheroid: V = 4π/3 * a^2 * c + @testset "n=1 reduction to spheroid" begin + a, c = 1.0, 2.0 + s = SuperSpheroid(a, c, 1.0, complex(1.5)) + @test volume(s) ≈ 4π / 3 * a^2 * c rtol=1e-8 + @test volume_equivalent_radius(s) ≈ ∛(a^2 * c) rtol=1e-8 + end + + # `volume` keeps the shape's real type and its full precision for every supported + # numeric type (regression guard for the Arb→T beta computation: no BigFloat detour, + # instance precision honoured, and the Quadmath `precision(::Float128)` workaround). + @testset "volume type & precision for $T" for (T, CT, mk) in ( + (Float64, ComplexF64, 1e-14), + (Double64, Complex{Double64}, 1e-28), + (Float128, ComplexF128, 1e-30), + (BigFloat, Complex{BigFloat}, 1e-70)) + a, c = T(12) / 10, T(8) / 10 + s = SuperSpheroid(a, c, one(T), complex(CT(3) / 2)) + v = volume(s) + @test v isa T + @test abs(v - 4 * T(π) / 3 * a^2 * c) < mk * abs(v) + end + + # n=2 reduces to octahedron: V = 4/3 * a^2 * c + @testset "n=2 reduction to octahedron" begin + a, c = 1.0, 1.0 + s = SuperSpheroid(a, c, 2.0, complex(1.5)) + @test volume(s) ≈ 4 / 3 * a^2 * c rtol=1e-8 + end + + # Volume vs Monte-Carlo ∈ test + @testset "Volume vs Monte-Carlo for n=$n" for n in [1.2, 1.7, 2.5] + a, c = 1.2, 0.8 + s = SuperSpheroid(a, c, n, complex(1.5)) + V_analytic = volume(s) + + rng_state = 42 + N = 200_000 + # Sample within bounding box [-a,a]^2 x [-c,c] + hits = 0 + let seed = UInt64(rng_state) + # Simple LCG for reproducibility + x = zeros(3) + for _ in 1:N + seed = seed * 6364136223846793005 + 1442695040888963407 + x[1] = (Float64(seed >> 11) / 2^53 * 2 - 1) * a + seed = seed * 6364136223846793005 + 1442695040888963407 + x[2] = (Float64(seed >> 11) / 2^53 * 2 - 1) * a + seed = seed * 6364136223846793005 + 1442695040888963407 + x[3] = (Float64(seed >> 11) / 2^53 * 2 - 1) * c + hits += x ∈ s ? 1 : 0 + end + end + V_box = 2a * 2a * 2c + V_mc = V_box * hits / N + @test abs(V_analytic - V_mc) / V_analytic < 0.01 + end + + @test has_symmetric_plane(SuperSpheroid(1.0, 1.0, 1.5, complex(1.5))) +end + +function rmax(s::SuperSpheroid) + # Equatorial max is a (at cardinal phi=0,pi/2), polar max is c. + # For n>=1 the equatorial maximum is a (super-circle corners are at axes). + # For n<1 the diagonal is larger, but IITM convergence typically uses n>=1. + n = Float64(s.n) + if n >= 1 + return max(s.a, s.c) + else + # equatorial diagonal: a * 2^((1-n)/2) + return max(s.a * oftype(s.a, 2)^((1 - n) / 2), s.c) + end +end + +function rmin(s::SuperSpheroid) + # Inscribed-sphere radius = min over directions û of the surface radius r(û). There is + # NO simple closed form: for concave superquadrics (n>2) with a≠c the closest surface + # point sits at a transcendental, off-symmetry direction (the symmetry-direction minimum + # over-estimates by a few %). So minimise r(û) over a dense direction grid. An over-large + # rₘᵢₙ — e.g. min(a,c), which ignores the body-diagonal/face dimples — makes the IITM + # Mie-init sphere poke outside the particle and gives wrong cross sections, so a small + # shrink (below) keeps it safely inscribed. + n = Float64(s.n) + a = Float64(s.a) + c = Float64(s.c) + q = 2 / n + M = 90 + rmn = Inf + for i in 0:M, j in 0:M + ϑ = i * (π / 2) / M + φ = j * (π / 2) / M + ux = sin(ϑ) * cos(φ) + uy = sin(ϑ) * sin(φ) + uz = cos(ϑ) + g = (ux / a)^q + (uy / a)^q + (uz / c)^q + rmn = min(rmn, g^(-n / 2)) + end + # The grid minimum can only OVER-estimate the true minimum (it may miss the exact + # closest direction between nodes, by ≲1e-3 here). rₘᵢₙ must never exceed the true + # inscribed radius — an over-large value pokes the Mie-init sphere outside the + # particle — so shrink by a margin comfortably larger than the grid error. + return oftype(s.a, 0.99 * rmn) +end + +function Base.:∈(x, s::SuperSpheroid) + return (abs(x[1]) / s.a)^(2 / s.n) + + (abs(x[2]) / s.a)^(2 / s.n) + + (abs(x[3]) / s.c)^(2 / s.n) <= 1 +end + +refractive_index(s::SuperSpheroid, x) = x ∈ s ? s.m : one(s.m) + +@testitem "SuperSpheroid n=1 reduces to a spheroid (IITM)" begin + using TransitionMatrices: SuperSpheroid, Spheroid, transition_matrix, IITM, + calc_Csca, calc_Cext + + m = complex(1.5) + # n=1 ⇒ the super-spheroid IS a spheroid. Compare both through the SAME IITM solver + # with a matched explicit rₘᵢₙ: this checks the geometric reduction directly (to + # machine precision) and is robust to the rₘᵢₙ default and to cross-method + # (EBCM-vs-IITM) discretisation differences. + slv = IITM(5, 60, 200, 100; rₘᵢₙ = 0.6) + 𝐓sph = transition_matrix(Spheroid(1.0, 2.0, m), 2π, slv) + 𝐓ss = transition_matrix(SuperSpheroid(1.0, 2.0, 1.0, m), 2π, slv) + + @test abs(calc_Csca(𝐓sph) - calc_Csca(𝐓ss)) < 1e-10 + @test abs(calc_Cext(𝐓sph) - calc_Cext(𝐓ss)) < 1e-10 +end + +@testitem "SuperSpheroid geometry predicates (incl. n<1 rmax branch)" begin + using TransitionMatrices: SuperSpheroid + for n in (0.7, 1.0, 2.5) # n<1 hits the diagonal-extended rmax branch + s = SuperSpheroid(1.0, 2.0, n, complex(1.5)) + @test TransitionMatrices.rmax(s) ≥ TransitionMatrices.rmin(s) > 0 + @test [0.0, 0.0, 0.0] ∈ s + @test [5.0, 0.0, 0.0] ∉ s + @test TransitionMatrices.refractive_index(s, [0.0, 0.0, 0.0]) == complex(1.5) + @test TransitionMatrices.refractive_index(s, [5.0, 0.0, 0.0]) == one(complex(1.5)) + end +end diff --git a/src/shapes/superspheroidrevolved.jl b/src/shapes/superspheroidrevolved.jl new file mode 100644 index 0000000..27a3443 --- /dev/null +++ b/src/shapes/superspheroidrevolved.jl @@ -0,0 +1,269 @@ +@doc raw""" +An axisymmetric super-spheroid scatterer (surface of revolution of a superellipse). + +This is the surface of revolution + +```math +\left(\frac{\varrho}{a}\right)^{2/n} + \left(\frac{z}{c}\right)^{2/n} = 1, +\qquad \varrho = \sqrt{x^2 + y^2} +``` + +Unlike `SuperSpheroid` (which has 4-fold D4h symmetry from super-circular cross-sections), +this shape is **truly axisymmetric** and is compatible with the EBCM / Sh-matrix solvers. + +At ``n = 1`` it reduces exactly to `Spheroid(a, c)`. + +Conventions: +- ``n = 1``: ordinary spheroid. +- ``n = 2``: bicone (diamond/spindle shape). +- ``1 < n < 2``: convex rounded. +- ``n > 2``: concave (waisted), analogous to the Lin/Bi dust regime. + +Attributes: + +- `a`: the equatorial semi-axis. +- `c`: the polar semi-axis. +- `n`: the roundness exponent. +- `m`: the relative complex refractive index. +""" +struct SuperSpheroidRevolved{T, CT} <: AbstractAxisymmetricShape{T, CT} + a::T + c::T + n::T + m::CT +end + +SuperSpheroidRevolved(a, c, n, m) = SuperSpheroidRevolved{typeof(a), typeof(m)}(a, c, n, m) + +@doc raw""" +Volume of the axisymmetric super-spheroid (surface of revolution): + +```math +V = \pi a^2 c \, n \, B\!\left(\tfrac{n}{2},\, n+1\right) +``` + +At ``n = 1``: ``B(1/2, 2) = 4/3``, so ``V = \tfrac{4}{3}\pi a^2 c`` (spheroid). +At ``n = 2``: ``B(1, 3) = 1/3``, so ``V = \tfrac{2}{3}\pi a^2 c`` (bicone). +""" +function volume(s::SuperSpheroidRevolved{T}) where {T} + n = s.n + b = _beta_lgamma(n / 2, n + 1) + return oftype(s.a, π) * s.a^2 * s.c * n * b +end + +volume_equivalent_radius(s::SuperSpheroidRevolved) = ∛(3 * volume(s) / (4 * oftype(s.a, π))) + +has_symmetric_plane(::SuperSpheroidRevolved) = true + +@testitem "SuperSpheroidRevolved utility functions" begin + using TransitionMatrices: SuperSpheroidRevolved, Spheroid, volume, + volume_equivalent_radius, has_symmetric_plane + + # n=1 must exactly match Spheroid volume + @testset "n=1 reduction to spheroid" begin + a, c = 1.0, 2.0 + s = SuperSpheroidRevolved(a, c, 1.0, complex(1.5)) + sph = Spheroid(a, c, complex(1.5)) + @test volume(s) ≈ volume(sph) rtol=1e-8 + @test volume_equivalent_radius(s) ≈ volume_equivalent_radius(sph) rtol=1e-8 + end + + # n=2 → bicone: V = 2π/3 * a^2 * c + @testset "n=2 reduction to bicone" begin + a, c = 1.0, 1.5 + s = SuperSpheroidRevolved(a, c, 2.0, complex(1.5)) + @test volume(s) ≈ 2π / 3 * a^2 * c rtol=1e-8 + end + + # Volume vs numeric integration (Simpson on the revolved surface) + @testset "Volume vs quadrature for n=$n" for n in [1.2, 1.7, 2.5] + a, c = 1.2, 0.8 + s = SuperSpheroidRevolved(a, c, n, complex(1.5)) + V_analytic = volume(s) + + # V = π ∫_{-c}^{c} ϱ(z)^2 dz where (ϱ/a)^(2/n) + (z/c)^(2/n) = 1 + # ϱ(z) = a * (1 - (|z|/c)^(2/n))^(n/2) + Npts = 10_000 + dz = 2c / Npts + V_num = 0.0 + for i in 0:(Npts - 1) + z = -c + (i + 0.5) * dz + t = abs(z) / c + t < 1 || continue + ϱ = a * (1 - t^(2 / n))^(n / 2) + V_num += π * ϱ^2 * dz + end + @test abs(V_analytic - V_num) / V_analytic < 1e-3 + end + + @test has_symmetric_plane(SuperSpheroidRevolved(1.0, 1.0, 1.5, complex(1.5))) +end + +@doc raw""" +``` +gaussquad(s::SuperSpheroidRevolved{T}, ngauss) where {T} +``` + +Evaluate the quadrature points, weights and the corresponding radius ``r(\vartheta)`` +and its derivative ``r'(\vartheta) = \mathrm{d}r/\mathrm{d}\vartheta`` for an +axisymmetric super-spheroid. + +The surface in spherical coordinates is defined by ``r(\vartheta)`` satisfying + +```math +\left(\frac{r\sin\vartheta}{a}\right)^{2/n} + \left(\frac{r\cos\vartheta}{c}\right)^{2/n} = 1 +``` + +Let ``p = 2/n`` and ``g(\vartheta) = (\sin\vartheta/a)^p + (\cos\vartheta/c)^p``. Then + +```math +r(\vartheta) = g(\vartheta)^{-1/p}, \qquad +r'(\vartheta) = -g^{-1/p-1}\left[\frac{\sin^{p-1}\vartheta\cos\vartheta}{a^p} + - \frac{\cos^{p-1}\vartheta\sin\vartheta}{c^p}\right] +``` + +Returns ``(x, w, r, r')`` where ``x = \cos\vartheta`` are Gauss–Legendre nodes on ``[-1,1]``. +""" +function gaussquad(s::SuperSpheroidRevolved{T}, ngauss) where {T} + x, w = gausslegendre(T, ngauss) + r = similar(x) + r′ = similar(x) + + p = 2 / s.n # exponent + ap = s.a^p + cp = s.c^p + + @simd for i in 1:ngauss + cosϑ = x[i] + sinϑ = √(1 - cosϑ^2) + + # g(ϑ) = (sinϑ/a)^p + (cosϑ/c)^p (using |cosϑ|: pole symmetry handled below) + # Use absolute values for generality; z-symmetry from has_symmetric_plane + sinp = abs(sinϑ)^p / ap + cosp = abs(cosϑ)^p / cp + g = sinp + cosp + + # r = g^(-1/p) = g^(-n/2) + r[i] = g^(-s.n / 2) + + # r′ = -g^(-1/p - 1) * [sinϑ^(p-1)*cosϑ/a^p - cosϑ^(p-1)*sinϑ/c^p] + # = -g^(-n/2 - 1) * [sign(cosϑ)*|sinϑ|^(p-1)*|cosϑ|^1/a^p + # - sign(sinϑ)*|cosϑ|^(p-1)*|sinϑ|^1/c^p] + # For ϑ ∈ (0,π): sinϑ ≥ 0, sign(sinϑ)=+1 + # cosϑ can be negative (ϑ > π/2) + if sinϑ == 0 + # Defensive pole guard: gausslegendre nodes are interior (|cosϑ|<1), so + # sinϑ>0 always and this branch is unreachable via `gaussquad` — excluded + # from coverage rather than removed (keeps r′ finite if ever called at a pole). + r′[i] = zero(T) # COV_EXCL_LINE + else + # sinϑ > 0 for i not at exactly the pole + dg_dϑ = (p * abs(sinϑ)^(p - 1) * cosϑ / ap - + p * abs(cosϑ)^(p - 1) * sinϑ * sign(cosϑ) / cp) + # r′ = dr/dϑ = -1/p * g^(-1/p-1) * dg/dϑ = -(n/2) * g^(-n/2-1) * dg/dϑ + r′[i] = -(s.n / 2) * g^(-s.n / 2 - 1) * dg_dϑ + end + end + + return x, w, r, r′ +end + +@testitem "SuperSpheroidRevolved gaussquad matches Spheroid at n=1" begin + using TransitionMatrices: SuperSpheroidRevolved, Spheroid, gaussquad + + a, c = 1.0, 2.0 + s = SuperSpheroidRevolved(a, c, 1.0, complex(1.5)) + sph = Spheroid(a, c, complex(1.5)) + ngauss = 60 + + _, _, r_s, rp_s = gaussquad(s, ngauss) + _, _, r_sph, rp_sph = gaussquad(sph, ngauss) + + @test all(r_s .≈ r_sph) + @test all(rp_s .≈ rp_sph) +end + +@testitem "SuperSpheroidRevolved gaussquad r' sign" begin + using TransitionMatrices: SuperSpheroidRevolved, gaussquad + + # Prolate (c > a): at theta=pi (costheta=-1, south pole), dr/dtheta > 0 + # (r decreases from c to a as theta goes from pi to pi/2) + @testset "prolate a=$a c=$c n=$n" for (a, c, n) in [(1.0, 2.0, 1.0), (1.0, 2.0, 1.5)] + s = SuperSpheroidRevolved(a, c, n, complex(1.5)) + _, _, _, r′ = gaussquad(s, 100) + @test r′[1] > 0 # near theta=pi (costheta ≈ -1) + @test r′[end] < 0 # near theta=0 (costheta ≈ +1) + end +end + +function rmax(s::SuperSpheroidRevolved) + # r(ϑ) = g(ϑ)^(-n/2); extremes occur at poles (ϑ=0,π) and equator (ϑ=π/2) + # r at pole: g = (1/c)^p → r = c + # r at equator: g = (1/a)^p → r = a + # For n>2 (concave), rmax may occur at an interior ϑ; compute numerically. + return _superspheroidrevolved_rmax_rmin(s)[1] +end + +function rmin(s::SuperSpheroidRevolved) + return _superspheroidrevolved_rmax_rmin(s)[2] +end + +function _superspheroidrevolved_rmax_rmin(s::SuperSpheroidRevolved) + p = 2 / s.n + ap = s.a^p + cp = s.c^p + # Sample r over ϑ ∈ [0, π/2] (symmetric about equator and pole) + N = 1000 + rvals = [begin + ϑ = i * π / 2 / N + sinϑ = sin(ϑ) + cosϑ = cos(ϑ) + g = sinϑ^p / ap + cosϑ^p / cp + g^(-s.n / 2) + end + for i in 0:N] + return maximum(rvals), minimum(rvals) +end + +function Base.:∈(x, s::SuperSpheroidRevolved) + ϱ = √(x[1]^2 + x[2]^2) + return (ϱ / s.a)^(2 / s.n) + (abs(x[3]) / s.c)^(2 / s.n) <= 1 +end + +refractive_index(s::SuperSpheroidRevolved, x) = x ∈ s ? s.m : one(s.m) + +@testitem "SuperSpheroidRevolved n=1 EBCM T-matrix matches Spheroid" begin + using TransitionMatrices: SuperSpheroidRevolved, Spheroid, EBCM, + transition_matrix, calc_Csca, calc_Cext + + a, c = 1.0, 2.0 + m = complex(1.5) + + sph = Spheroid(a, c, m) + ssr = SuperSpheroidRevolved(a, c, 1.0, m) + + nₘₐₓ = 5 + Ng = 80 + + 𝐓sph = transition_matrix(sph, 2π, EBCM(nₘₐₓ, Ng)) + 𝐓ssr = transition_matrix(ssr, 2π, EBCM(nₘₐₓ, Ng)) + + @test abs(calc_Csca(𝐓sph) - calc_Csca(𝐓ssr)) < 1e-10 + @test abs(calc_Cext(𝐓sph) - calc_Cext(𝐓ssr)) < 1e-10 +end + +@testitem "SuperSpheroidRevolved geometry predicates (∈, rmax, rmin, refractive_index)" begin + using TransitionMatrices: SuperSpheroidRevolved + # span n<1, ellipse, convex, concave so rmax/rmin sampling + ∈ are exercised + for n in (0.8, 1.0, 1.5, 2.5) + s = SuperSpheroidRevolved(1.0, 2.0, n, complex(1.5)) + rmx = TransitionMatrices.rmax(s) + rmn = TransitionMatrices.rmin(s) + @test rmx ≥ rmn > 0 + @test [0.0, 0.0, 0.0] ∈ s # centre inside + @test [5.0, 0.0, 0.0] ∉ s # far outside + @test [0.0, 0.0, 0.99 * s.c] ∈ s # near a pole, inside + @test TransitionMatrices.refractive_index(s, [0.0, 0.0, 0.0]) == complex(1.5) + @test TransitionMatrices.refractive_index(s, [5.0, 0.0, 0.0]) == one(complex(1.5)) + end +end