From 0ce4b8a8699c2bb0f5cc81b9189bf1285e9b2034 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Mon, 30 Oct 2023 14:17:10 -0400
Subject: [PATCH] Add a method for rk4 3D CartesianGrid

---
 src/structured3d.jl | 21 +++++++++++++++++++++
 test/runtests.jl    |  6 ++++++
 2 files changed, 27 insertions(+)

diff --git a/src/structured3d.jl b/src/structured3d.jl
index 43fe5c6..5fba956 100644
--- a/src/structured3d.jl
+++ b/src/structured3d.jl
@@ -334,4 +334,25 @@ function trace3d_euler(fieldx::F, fieldy::F, fieldz::F, startx::T, starty::T, st
 
    trace3d_euler(fieldx, fieldy, fieldz, startx, starty, startz, gridx, gridy, gridz;
       kwargs...)
+end
+
+"""
+    trace3d_rk4(fieldx, fieldy, fieldz, startx, starty, startz, grid::CartesianGrid;
+		 maxstep=20000, ds=0.01, gridtype="ndgrid", direction="both")
+
+See also [`trace3d_euler`](@ref).
+"""
+function trace3d_rk4(fieldx::F, fieldy::F, fieldz::F, startx::T, starty::T, startz::T,
+   grid::CartesianGrid; kwargs...) where {F, T}
+
+   gridmin = coordinates(minimum(grid))
+   gridmax = coordinates(maximum(grid))
+   Δx = spacing(grid)
+
+   gridx = range(gridmin[1], gridmax[1], step=Δx[1])
+   gridy = range(gridmin[2], gridmax[2], step=Δx[2])
+   gridz = range(gridmin[3], gridmax[3], step=Δx[3])
+
+   trace3d_rk4(fieldx, fieldy, fieldz, startx, starty, startz, gridx, gridy, gridz;
+      kwargs...)
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index ab3d684..2a49540 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -96,6 +96,12 @@ using Test
          alg=Euler(), ds=0.2, maxstep=200)
 
       @test length(x1) == 170 && x1[end] ≈ y1[end] ≈ z1[end]
+
+      # RK4 by default
+      x1, y1, z1 = trace(bx, bz, bz, xs, ys, zs, x, y, z;
+         ds=0.2, maxstep=200, direction="forward")
+   
+      @test length(x1) == 152 && x1[end] ≈ y1[end] ≈ z1[end]
    end
 
    @testset "2D unstructured mesh" begin