Add a method for rk4 3D CartesianGrid

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

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