Ensemble ray tracing

src/01_oceanography/oceanography.jl

@@ -0,0 +1,2 @@
+"Circumference [m] of earth's equator."
+const EQUATORIAL_EARTH_CIRCUMFERENCE = 40_075e3

src/02_acoustics/03_propagation/tracing/03_beam.jl

@@ -0,0 +1,91 @@
+export Beam
+
+const DEFAULT_RAY_ARC_LENGTH_SPAN = (0, EQUATORIAL_EARTH_CIRCUMFERENCE)
+
+function AcousticTracingODESystem(
+    c_ocn::Function,
+    r_max::Real;
+    name
+)
+    @parameters begin
+        s, [description = "Ray arc length"]
+    end
+
+    pars = @parameters begin
+        r₀ = 0.0
+        z₀
+        θ₀
+    end
+
+    deps = @variables begin
+		r(s) = r₀
+		z(s) = z₀
+		ξ(s) = cos(θ₀) / c_ocn(r₀, z₀)
+		ζ(s) = sin(θ₀) / c_ocn(r₀, z₀)
+        θ(s), [description = "Ray angle, positive downwards"]
+        c(s), [description = "Ray celerity"]
+    end
+
+    Ds = Differential(s)
+
+    c²(r′, z′) = c_ocn(r′, z′)^2
+    ∂c_∂r(r′, z′) = ModelingToolkit.derivative(c_ocn(r, z′), r′)
+    ∂c_∂z(r′, z′) = ModelingToolkit.derivative(c_ocn(r′, z), z′)
+
+    eqns = [
+        Ds(r) ~ c_ocn(r, z) * ξ
+        Ds(z) ~ c_ocn(r, z) * ζ
+		Ds(ξ) ~ -∂c_∂r(r, z) / c_ocn(r, z)^2
+		Ds(ζ) ~ -∂c_∂z(r, z) / c_ocn(r, z)^2
+        θ ~ atan(ζ, ξ)
+        c ~ c_ocn(r, z)
+    ]
+
+    maximum_range = [r ~ r_max] => (
+        (ntg, _, _, _) -> terminate!(ntg),
+        [], [], []
+    )
+
+    events = [
+        maximum_range
+    ]
+
+    return ODESystem(eqns, s, deps, pars;
+        name = name,
+        continuous_events = events
+    )
+end
+
+struct Beam{MN <: ModelName}
+    model::MN
+
+    s_max::Float64
+
+    c::Function
+    θ::Function
+    r::Function
+    z::Function
+    ξ::Function
+    ζ::Function
+
+    sol::ODESolution
+
+    function Beam(model::ModelName, sys::ODESystem, sol::ODESolution)
+        r(s::Real) = sol(s, idxs = sys.r)
+        r(s::AbstractVector{<:Real}) = sol(s, idxs = sys.r) |> collect
+        z(s::Real) = sol(s, idxs = sys.z)
+        z(s::AbstractVector{<:Real}) = sol(s, idxs = sys.z) |> collect
+        ξ(s::Real) = sol(s, idxs = sys.ξ)
+        ξ(s::AbstractVector{<:Real}) = sol(s, idxs = sys.ξ) |> collect
+        ζ(s::Real) = sol(s, idxs = sys.ζ)
+        ζ(s::AbstractVector{<:Real}) = sol(s, idxs = sys.ζ) |> collect
+        c(s::Real) = sol(s, idxs = sys.c)
+        c(s::AbstractVector{<:Real}) = sol(s, idxs = sys.c) |> collect
+        θ(s::Real) = sol(s, idxs = sys.θ)
+        θ(s::AbstractVector{<:Real}) = sol(s, idxs = sys.θ) |> collect
+
+        new{model |> typeof}(model, sol.t[end], c, θ, r, z, ξ, ζ, sol)
+    end
+end
+
+show(io::IO, beam::Beam) = print(io, "Beam($(beam.θ(0) |> rad2deg)°)")

src/02_acoustics/03_propagation/tracing/04_fan.jl

@@ -0,0 +1,45 @@
+export Fan
+
+struct Fan{MN <: ModelName}
+    beams::Vector{Beam{MN}}
+    sys::ODESystem
+    prob::EnsembleProblem
+
+    function Fan(
+        model::ModelName,
+        θ₀s::AbstractVector{<:Real},
+        z₀::Real,
+        r_max::Real,
+        c_ocn::Function
+    )
+    @mtkbuild sys = AcousticTracingODESystem(c_ocn, r_max)
+    prob = ODEProblem(sys, [], DEFAULT_RAY_ARC_LENGTH_SPAN, [sys.θ₀ => 0.0, sys.z₀ => z₀])
+
+    # Attempt 1: Errors on `init`.
+    # tiltray! = setp(sys, sys.θ₀)
+    # prob_func(internal_prob, n, _) = tiltray!(internal_prob, θ₀s[n])
+
+    # Attempt 2: Doesn't change launch angle parameter
+    # function prob_func(internal_prob, n, _)
+    #     new_prob = remake(internal_prob, p = [sys.θ₀ => θ₀s[n], sys.z₀ => z₀])
+    #     @show prob.p[1]
+    #     return new_prob
+    # end
+
+    # Attempt 3: Make an entirely new ODEProblem
+    prob_func(_, n, _) = ODEProblem(sys, [], DEFAULT_RAY_ARC_LENGTH_SPAN, [sys.θ₀ => θ₀s[n], sys.z₀ => z₀])
+
+    ens_prob = EnsembleProblem(prob, prob_func = prob_func)
+    ens_sol = solve(ens_prob, Tsit5(), EnsembleThreads(),
+        reltol = 1e-50,
+        trajectories = length(θ₀s)
+    )
+        beams = [Beam(model, sys, sol) for sol in ens_sol]
+        new{model |> typeof}(beams, sys, ens_prob)
+    end
+end
+
+function show(io::IO, fan::Fan)
+    xtr = [beam.θ(0) for beam in fan.beams] |> extrema .|> rad2deg
+    print(io, "Fan($(fan.beams |> length): $(xtr |> minimum)° .. $(xtr |> maximum)°)")
+end

src/OceanSonar.jl

@@ -2,17 +2,30 @@ module OceanSonar
 
 export subtypes
 
+import Base: getproperty, show
+
 using InteractiveUtils: subtypes
 
 using ModelingToolkit:
     ModelingToolkit,
     Differential,
+    EnsembleProblem,
+    EnsembleSolution,
+    EnsembleThreads,
+    ODEProblem,
+    ODESolution,
     ODESystem,
+    remake,
+    # setp,
     terminate!,
     @mtkbuild,
     @parameters,
     @variables
 
+using OrdinaryDiffEq:
+    solve,
+    Tsit5
+
 """
 ```
 include_subroots(path::AbstractString)