Ray tracing works

Project.toml

@@ -5,6 +5,8 @@ version = "1.0.0-DEV"
 
 [deps]
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 
 [weakdeps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -15,4 +17,6 @@ OceanSonarAbstractTreesExt = "AbstractTrees"
 [compat]
 AbstractTrees = "0.4"
 InteractiveUtils = "1.11.0"
+ModelingToolkit = "9.30.0"
+OrdinaryDiffEq = "6.87.0"
 julia = "1.10"

docs/Project.toml

@@ -1,4 +1,5 @@
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 OceanSonar = "03019ade-4524-4ecd-af79-46d4f04a1b56"

src/02_acoustics/03_propagation/tracing.jl

@@ -0,0 +1,58 @@
+function TracingODESystem(
+    c_ocn::Function,
+    r_max::Real;
+    name
+)
+    @parameters begin
+        s, [description = "Ray arc length"]
+    end
+
+    pars = @parameters begin
+        θ₀
+        r₀ = 0.0
+        z₀
+        c₀ = c_ocn(r₀, z₀)
+        ξ₀ = cos(θ₀) / c₀
+        ζ₀ = sin(θ₀) / c₀
+    end
+
+    @variables begin
+        θ(s), [description = "Ray angle, positive downwards"]
+        c(s), [description = "Ray celerity"]
+    end
+
+    deps = @variables begin
+		r(s) = r₀
+		z(s) = z₀
+		ξ(s) = ξ₀
+		ζ(s) = ζ₀
+    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
+    )
+end

src/OceanSonar.jl

@@ -4,6 +4,14 @@ export subtypes
 
 using InteractiveUtils: subtypes
 
+using ModelingToolkit:
+    ModelingToolkit,
+    Differential,
+    ODESystem,
+    @mtkbuild,
+    @parameters,
+    @variables
+
 """
 ```
 include_subroots(path::AbstractString)
@@ -37,6 +45,6 @@ end
 
 include_subroots(@__FILE__)
 
-const TRIGGER_RECOMPILATION_BY_CHANGING_THIS_NUMBER = 0
+const TRIGGER_RECOMPILATION_BY_CHANGING_THIS_NUMBER = 1
 
 end