From fb79b8d5f3baf6db37c21d06519d2d5f8615c3ce Mon Sep 17 00:00:00 2001
From: Aaron Kaw <aaronjkaw@gmail.com>
Date: Tue, 6 Aug 2024 19:59:31 +1000
Subject: [PATCH] Ray tracing works

---
 Project.toml                               |  4 ++
 docs/Project.toml                          |  1 +
 src/02_acoustics/03_propagation/tracing.jl | 56 ++++++++++++++++++++++
 src/OceanSonar.jl                          | 11 ++++-
 4 files changed, 71 insertions(+), 1 deletion(-)
 create mode 100644 src/02_acoustics/03_propagation/tracing.jl

diff --git a/Project.toml b/Project.toml
index f12cb91..f5627af 100644
--- a/Project.toml
+++ b/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"
diff --git a/docs/Project.toml b/docs/Project.toml
index 66d4c6b..e877161 100644
--- a/docs/Project.toml
+++ b/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"
diff --git a/src/02_acoustics/03_propagation/tracing.jl b/src/02_acoustics/03_propagation/tracing.jl
new file mode 100644
index 0000000..71ba7f6
--- /dev/null
+++ b/src/02_acoustics/03_propagation/tracing.jl
@@ -0,0 +1,56 @@
+function TracingODESystem(
+    c_ocn::Function,
+    r_max::Real;
+    name
+)
+    @parameters begin
+        s, [description = "Ray arc length"]
+    end
+
+    pars = @parameters begin
+        r₀ = 0.0
+        z₀
+        θ₀
+    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) = cos(θ₀) / c_ocn(r₀, z₀)
+		ζ(s) = sin(θ₀) / c_ocn(r₀, z₀)
+    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
diff --git a/src/OceanSonar.jl b/src/OceanSonar.jl
index d9e1ff5..d9dc947 100644
--- a/src/OceanSonar.jl
+++ b/src/OceanSonar.jl
@@ -4,6 +4,15 @@ export subtypes
 
 using InteractiveUtils: subtypes
 
+using ModelingToolkit:
+    ModelingToolkit,
+    Differential,
+    ODESystem,
+    terminate!,
+    @mtkbuild,
+    @parameters,
+    @variables
+
 """
 ```
 include_subroots(path::AbstractString)
@@ -37,6 +46,6 @@ end
 
 include_subroots(@__FILE__)
 
-const TRIGGER_RECOMPILATION_BY_CHANGING_THIS_NUMBER = 0
+const TRIGGER_RECOMPILATION_BY_CHANGING_THIS_NUMBER = 1
 
 end
\ No newline at end of file