Skip to content

Commit 3d68a95

Browse files
committed
switch from using Dierckx to SmoothingSplines for cubic interpolation
1 parent 838056a commit 3d68a95

1 file changed

Lines changed: 15 additions & 11 deletions

File tree

src/Domain.jl

Lines changed: 15 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
11
using Parameters
22
using LinearAlgebra
33
using StaticArrays
4-
using Dierckx
54
using Calculus
5+
using SmoothingSplines
66

77
abstract type AbstractDomain end
88
export AbstractDomain
@@ -263,16 +263,14 @@ function ParametrizedTPDomain(;phase::Z,initialconds::Dict{X,Any},constantspecie
263263
end
264264
@assert V != 0.0 || (T != 0.0 && P != 0.0)
265265
if isa(T,AbstractArray)
266-
q = Spline1D(ts,T;k=3,s=1e-11)
267-
Tfcn(x::Float64) = q(x)
266+
Tfcn = getspline(ts,T)
268267
elseif isa(T,Function)
269268
Tfcn = T
270269
else
271270
throw(error("ParametrizedTPDomain must take \"T\" as a function or if an array of times for \"ts\" is supplied as an array of volumes"))
272271
end
273272
if isa(P,AbstractArray)
274-
v = Spline1D(ts,P;k=3,s=1e-11)
275-
Pfcn(x::Float64) = v(x)
273+
Pfcn = getspline(ts,P)
276274
elseif isa(P,Function)
277275
Pfcn = P
278276
else
@@ -346,8 +344,7 @@ function ParametrizedVDomain(;phase::Z,initialconds::Dict{X,Any},constantspecies
346344
end
347345
@assert isa(V,Function) || isa(V,AbstractArray)
348346
if isa(V,AbstractArray)
349-
q = Spline1D(ts,V;k=3,s=1e-11)
350-
Vfcn = f(x::Float64) = q(x)
347+
Vfcn = getspline(ts,V)
351348
elseif isa(V,Function)
352349
Vfcn = V
353350
else
@@ -422,8 +419,7 @@ function ParametrizedPDomain(;phase::Z,initialconds::Dict{X,Any},constantspecies
422419
end
423420
@assert isa(P,Function) || isa(P,AbstractArray)
424421
if isa(P,AbstractArray)
425-
q = Spline1D(ts,P;k=3,s=1e-11)
426-
Pfcn = f(x::Float64) = q(x)
422+
Pfcn = getspline(ts,P)
427423
elseif isa(P,Function)
428424
Pfcn = P
429425
else
@@ -577,8 +573,7 @@ function ParametrizedTConstantVDomain(;phase::IdealDiluteSolution,initialconds::
577573
end
578574
end
579575
if isa(T,AbstractArray)
580-
q = Spline1D(ts,T;k=3,s=1e-11)
581-
Tfcn = f(x::Float64) = q(x)
576+
Tfcn = getspline(ts,T)
582577
elseif isa(T,Function)
583578
Tfcn = T
584579
else
@@ -943,3 +938,12 @@ function getreactionindices(ig::Q) where {Q<:AbstractPhase}
943938
return arr
944939
end
945940
export getreactionindices
941+
942+
"""
943+
fit a cubic spline to data and return a function evaluating that spline
944+
"""
945+
function getspline(xs,vals;s=1e-10)
946+
smspl = fit(SmoothingSpline,xs,vals,s)
947+
F(x::T) where {T} = predict(smspl,x)
948+
return F
949+
end

0 commit comments

Comments
 (0)