internal eval and integral added.

Project.toml

@@ -1,10 +1,18 @@
 name = "JPEC"
 uuid = "462872dd-e066-4d2e-b993-6468b5239634"
-authors = ["Matthew Pharr <matthewcpharr@gmail.com>"]
+authors = ["Matthew Pharr <matthewcpharr@gmail.com>", "JayMyoung Lee <jmlmir@snu.ac.kr>"]
 version = "0.1.0"
 
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
 [compat]
+BenchmarkTools = "1.6.0"
+DifferentialEquations = "7.16.1"
+LinearAlgebra = "1.11.0"
 Plots = "1.40.15"
+Printf = "1.11.0"

src/JPEC.jl

@@ -1,7 +1,10 @@
+# JPEC.jl
 module JPEC
 
 include("Splines/Splines.jl")
 using .SplinesMod
 
+#include("Equilibrium/Equilibrium.jl")
+#using .Equilibrium
 
-end
+end # module JPEC

src/Splines/BicubicSpline.jl

@@ -3,31 +3,60 @@ module BicubicSpline
 const libdir = joinpath(@__DIR__, "..", "..", "deps")
 const libspline = joinpath(libdir, "libspline")
 
+using ..Helper: parse_bctype, ReadOnlyArray, @expose_fields
+
 export bicube_setup, bicube_eval
 
 mutable struct BicubicSplineType
 	handle::Ptr{Cvoid}
-	xs::Vector{Float64}
-	ys::Vector{Float64}
-	fs::Array{Float64, 3} # 3D array for bicubic spline values
+	_xs::Vector{Float64}
+	_ys::Vector{Float64}
+	_fs::Array{Float64, 3} # 3D array for bicubic spline values
 	mx::Int64
 	my::Int64
 	nqty::Int64
-	ix::Int32  # Index of x position in the spline
-	iy::Int32  # Index of y position in the spline
 	bctypex::Int32  # Boundary condition type for x
 	bctypey::Int32  # Boundary condition type for y
+
+
+	_fsx::Array{Float64, 3}
+	_fsy::Array{Float64, 3}
+	_fsxy::Array{Float64, 3}
+
+
+
 end
 
-function _MakeBicubicSpline(mx::Int64, my::Int64, nqty::Int64)
+@expose_fields BicubicSplineType xs ys fs fsx fsy fsxy
+
+
+function _destroy_bicubic_spline(bicube::BicubicSplineType)
+    if bicube.handle != C_NULL
+        ccall((:bicube_c_destroy, libspline), Cvoid, (Ptr{Cvoid},), bicube.handle)
+		Core.setfield!(bicube, :handle, C_NULL)
+    end
+end
+
+function _MakeBicubicSpline(mx::Int64, my::Int64, nqty::Int64, bctypex::Int32, bctypey::Int32)
 	h = Ref{Ptr{Cvoid}}()
 	ccall((:bicube_c_create, libspline), Cvoid,
 		(Int64, Int64, Int64, Ref{Ptr{Cvoid}}), mx, my, nqty, h)
-	return BicubicSplineType(h[], Vector{Float64}(undef, mx), Vector{Float64}(undef, my), 
-		Array{Float64,3}(undef, mx, my, nqty), mx, my, nqty, 0, 0, 0, 0)
+
+	fsx = Array{Float64, 3}(undef, 0,0,0)
+	fsy = Array{Float64, 3}(undef, 0,0,0)
+	fsxy = Array{Float64, 3}(undef, 0,0,0)
+
+
+
+
+	return BicubicSplineType(h[], Vector{Float64}(undef, 0), Vector{Float64}(undef, 0), 
+		Array{Float64,3}(undef, 0, 0, 0), mx, my, nqty, bctypex, bctypey,
+		fsx,fsy, fsxy)
 end
 
 
+
+
 function _bicube_setup(xs::Vector{Float64}, ys::Vector{Float64}, fs::Array{Float64, 3}, bctypex::Int32, bctypey::Int32)
 	# xs -> Float64 (mx)
 	# ys -> Float64 (my)
@@ -38,40 +67,57 @@ function _bicube_setup(xs::Vector{Float64}, ys::Vector{Float64}, fs::Array{Float
 	mx = length(xs)-1
 	my = length(ys)-1
 	nqty = size(fs, 3)
-	bicube = _MakeBicubicSpline(mx, my, nqty)
-	bicube.xs = xs
-	bicube.ys = ys
-	bicube.fs = fs
-	bicube.bctypex = Int32(bctypex)
-	bicube.bctypey = Int32(bctypey)
+	bicube = _MakeBicubicSpline(mx, my, nqty, Int32(bctypex), Int32(bctypey))
+	bicube._xs = xs
+	bicube._ys = ys
+	bicube._fs = fs
+
+	bicube._fsx = Array{Float64, 3}(undef, mx+1, my+1, nqty)
+	bicube._fsy = Array{Float64, 3}(undef, mx+1, my+1, nqty)
+	bicube._fsxy = Array{Float64, 3}(undef, mx+1, my+1, nqty)
+
 
 	ccall((:bicube_c_setup, libspline), Cvoid,
 		(Ptr{Cvoid}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
 		bicube.handle, xs, ys, fs)
 
 	ccall((:bicube_c_fit, libspline), Cvoid,
-		(Ptr{Cvoid}, Int32, Int32), bicube.handle, bctypex, bctypey)
+		(Ptr{Cvoid}, Int32, Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}), 
+		bicube.handle, bicube.bctypex, bicube.bctypey, bicube._fsx, bicube._fsy, bicube._fsxy)
+
+
+
+
 	return bicube
 end
 
-function bicube_setup(xs, ys, fs, bctypex::Int=1, bctypey::Int=1)
+function bicube_setup(xs, ys, fs; bctypex::Union{String, Int}="not-a-knot",
+	 bctypey::Union{String, Int}="not-a-knot")
 	"""
 	# bicube_setup(xs, ys, fs, bctypex=0, bctypey=0)
 	## Arguments:
 	- `xs`: A vector of Float64 values representing the x-coordinates.
 	- `ys`: A vector of Float64 values representing the y-coordinates.
 	- `fs`: A 3D array of Float64 values representing the function values at the (x,y) coordinates.
-	- `bctypex`: An integer specifying the boundary condition type for x (default is 0).
-	- `bctypey`: An integer specifying the boundary condition type for y (default is 0).
+	## Keyword Arguments:
+	- `bctypex`: An integer specifying the boundary condition type for x (Default is 4, not a knot)
+	- `bctypey`: An integer specifying the boundary condition type for y  (Default is 4, not a knot) 
 	## Returns:
 	- A `BicubicSpline` object containing the spline handle, x-coordinates, y-coordinates,
 	function values, number of x-coordinates, number of y-coordinates, number of quantities,
 	and boundary condition types.
 	"""
+
+	local bctype_code_x::Int = parse_bctype(bctypex)
+	local bctype_code_y::Int = parse_bctype(bctypey)
+
 	if !isa(xs, Vector{Float64}) || !isa(ys, Vector{Float64}) || !isa(fs, Array{Float64, 3})
 		error("xs must be a vector of Float64, ys must be a vector of Float64, and fs must be a 3D array of Float64")
 	end
-	bicube = _bicube_setup(xs, ys, fs, Int32(bctypex), Int32(bctypey))
+	bicube = _bicube_setup(xs, ys, fs, Int32(bctype_code_x), Int32(bctype_code_y))
+
+	finalizer(_destroy_bicubic_spline, bicube)
+
 	return bicube
 end
 
@@ -87,16 +133,16 @@ function _bicube_eval(bicube::BicubicSplineType, x::Float64, y::Float64, derivs:
 	fyy = Vector{Float64}(undef, bicube.nqty)
 	if derivs == 0
 		ccall((:bicube_c_eval, libspline), Cvoid,
-			(Ptr{Cvoid}, Float64, Float64, Ptr{Float64}, Int32, Int32),
-			bicube.handle, x, y, f, bicube.ix, bicube.iy)
+			(Ptr{Cvoid}, Float64, Float64, Ptr{Float64}),
+			bicube.handle, x, y, f)
 	elseif derivs == 1
 		ccall((:bicube_c_eval_deriv, libspline), Cvoid,
-			(Ptr{Cvoid}, Float64, Float64, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int32, Int32),
-			bicube.handle, x, y, f, fx, fy, bicube.ix, bicube.iy)
+			(Ptr{Cvoid}, Float64, Float64, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
+			bicube.handle, x, y, f, fx, fy)
 	elseif derivs == 2
 		ccall((:bicube_c_eval_deriv2, libspline), Cvoid,
-			(Ptr{Cvoid}, Float64, Float64, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int32, Int32),
-			bicube.handle, x, y, f, fx, fy, fxx, fxy, fyy, bicube.ix, bicube.iy)
+			(Ptr{Cvoid}, Float64, Float64, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
+			bicube.handle, x, y, f, fx, fy, fxx, fxy, fyy)
 	else
 		error("Invalid number of derivatives requested: $derivs. Must be 0, 1, or 2.")
 	end
@@ -111,6 +157,7 @@ function _bicube_eval(bicube::BicubicSplineType, x::Float64, y::Float64, derivs:
 	end
 end
 
+
 function _bicube_eval(bicube::BicubicSplineType, xs::Vector{Float64}, ys::Vector{Float64}, derivs::Int=0)
 	# xs -> Float64 (any length)
 	# ys -> Float64 (any length)
@@ -180,6 +227,9 @@ function _bicube_eval(bicube::BicubicSplineType, xs::Vector{Float64}, ys::Vector
 	end
 end
 
+
+
+
 function bicube_eval(bicube::BicubicSplineType, x, y, derivs::Int=0)
 	"""
 	# bicube_eval(bicube, x, y)
@@ -195,4 +245,6 @@ function bicube_eval(bicube::BicubicSplineType, x, y, derivs::Int=0)
 	return _bicube_eval(bicube, x, y, derivs)
 end
 
+
+
 end # module BicubicSpline