From 686b26d348fc1c0add17a2cf9dbca8737e08374f Mon Sep 17 00:00:00 2001
From: Christian Haack <chr.hck@gmail.come>
Date: Thu, 13 Feb 2025 14:51:24 +0100
Subject: [PATCH] refactor scattering functions for improved type handling and
 clarity in documentation

---
 src/CherenkovMediumBase.jl |  4 ++--
 src/dispersion.jl          | 39 ++++++++------------------------------
 src/scattering.jl          |  8 +++++---
 3 files changed, 15 insertions(+), 36 deletions(-)

diff --git a/src/CherenkovMediumBase.jl b/src/CherenkovMediumBase.jl
index 51ed94b..e41caf6 100644
--- a/src/CherenkovMediumBase.jl
+++ b/src/CherenkovMediumBase.jl
@@ -65,7 +65,7 @@ radiation_length(medium::MediumProperties) = _not_implemented(medium)
 """
     sample_scattering_function(medium::MediumProperties)
 
-Return a scattering angle sampled from the scattering function of the medium.
+Return a cos(scattering angle) sampled from the scattering function of the medium.
 """
 function sample_scattering_function(medium::MediumProperties)
     model = get_scattering_model(medium)
@@ -101,7 +101,7 @@ Return the phase refractive index at `wavelength`.
 `wavelength` is expected to be in units nm.
 """
 function phase_refractive_index(medium::MediumProperties, wavelength)
-    model =get_dispersion_model(medium)
+    model = get_dispersion_model(medium)
     return phase_refractive_index(model, wavelength)
 end
 
diff --git a/src/dispersion.jl b/src/dispersion.jl
index 1847bec..59c9f12 100644
--- a/src/dispersion.jl
+++ b/src/dispersion.jl
@@ -25,9 +25,9 @@ Helper function to get the parameters for the Quan & Fry formula as function of
 salinity, temperature and pressure.
 """
 function calc_quan_fry_params(
-    salinity::Real,
-    temperature::Real,
-    pressure::Real)
+    salinity::T,
+    temperature::T,
+    pressure::T) where {T <: Real}
 
     n0 = 1.31405
     n1 = 1.45e-5
@@ -54,11 +54,12 @@ function calc_quan_fry_params(
     a3 = -n9
     a4 = n10
 
-    return a01, a2, a3, a4
+    return T(a01), T(a2), T(a3), T(a4)
 end
 
+
 """
-    refractive_index_fry(wavelength, salinity, temperature, pressure)
+    _refractive_index_fry(wavelength, quan_fry_params)
 
 The phase refractive index of sea water according to a model
 from Quan & Fry.
@@ -76,14 +77,6 @@ downloaded Jan 2011 from: http://www.physics.ox.ac.uk/Users/schuster/thesis0098m
 
 Adapted from clsim (https://github.com/claudiok/clsim)
 """
-function _refractive_index_fry(
-    wavelength::T;
-    salinity::Real,
-    temperature::Real,
-    pressure::Real) where {T<:Real}
-    refractive_index_fry(wavelength, T.(calc_quan_fry_params(salinity, temperature, pressure)))
-end
-
 function _refractive_index_fry(
     wavelength::Real,
     quan_fry_params::Tuple{U,U,U,U}
@@ -94,26 +87,10 @@ function _refractive_index_fry(
 end
 
 """
-    dispersion_fry(
-        wavelength::T;
-        salinity::Real,
-        temperature::Real,
-        pressure::Real) where {T <: Real}
+    dispersion_fry(wavelength, a2, a3, a4)
 
-Calculate the dispersion (dn/dλ) for the Quan & Fry model.
-Wavelength is given in nm, salinity in permille, temperature in °C and pressure in atm
+Calculate the dispersion for the Quan & Fry dispersion model.
 """
-function dispersion_fry(
-    wavelength::T;
-    salinity::Real,
-    temperature::Real,
-    pressure::Real) where {T<:Real}
-
-    quan_fry_params = T.(calc_quan_fry_params(salinity, temperature, pressure))
-
-    dispersion_fry(wavelength, quan_fry_params[2], quan_fry_params[3], quan_fry_params[4])
-end
-
 function dispersion_fry(wavelength::T, a2, a3, a4) where {T<:Number}
     x = one(T) / wavelength
 
diff --git a/src/scattering.jl b/src/scattering.jl
index 455ff85..1d00533 100644
--- a/src/scattering.jl
+++ b/src/scattering.jl
@@ -2,8 +2,8 @@ using Distributions
 using Random: AbstractRNG
 using Polynomials: fit, Polynomial, ImmutablePolynomial
 
+# Scattering functions
 export AbstractScatteringFunction
-export AbstractScatteringModel
 export HenyeyGreenStein
 export PolynomialScatteringFunction
 export EinsteinSmoluchowsky
@@ -12,6 +12,8 @@ export TwoComponentScatteringFunction
 export MixedHGES
 export MixedHGSL
 
+# Scattering models
+export AbstractScatteringModel
 export KopelevichScatteringModel
 export WavelengthIndependentScatteringModel
 
@@ -69,7 +71,7 @@ Implementation from: https://user-web.icecube.wisc.edu/~dima/work/WISC/ppc/spice
 function sl_scattering_func(rng::AbstractRNG, g::T) where {T <: Real}
     eta = Base.rand(rng, T)
     beta = (1-g) / (1+g)
-    costheta::T = 2 * eta^beta - 1
+    costheta::T = T(2) * eta^beta - T(1)
     return clamp(costheta, T(-1), T(1))
 end
 
@@ -138,7 +140,7 @@ end
 """
     TwoComponentScatteringFunction{F1<:AbstractScatteringFunction, F2<:AbstractScatteringFunction}
 
-Struct for two-component scattering model.
+Struct for two-component mixture scattering model.
 """
 struct TwoComponentScatteringFunction{F1<:AbstractScatteringFunction, F2<:AbstractScatteringFunction,T} <: AbstractScatteringFunction
     f1::F1