src/splines.f90: remove spline_simple
This commit is contained in:
Michele Guerini Rocco
parent
6e82aeeba7
commit
97906348ba
GPG Key ID:
BFBAF4C975F76450
136
src/beams.f90
136
src/beams.f90
@ -95,7 +95,7 @@ contains
|
||||
! ellipses in the transverse plane at the launch point (deg)
|
||||
|
||||
use gray_params, only : antenna_parameters
|
||||
use splines, only : spline_simple
|
||||
use splines, only : spline_1d, linear_1d
|
||||
use utils, only : locate
|
||||
use logger, only : log_error
|
||||
|
||||
@ -105,13 +105,11 @@ contains
|
||||
|
||||
! local variables
|
||||
integer :: u, nisteer, i, k, ii
|
||||
real(wp_) :: steer, dal
|
||||
real(wp_), dimension(:), allocatable :: &
|
||||
alphastv, betastv, x00v, y00v, &
|
||||
character(256) :: msg
|
||||
real(wp_) :: steer
|
||||
real(wp_), dimension(:), allocatable :: &
|
||||
alphastv, betastv, x00v, y00v, &
|
||||
z00v, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v
|
||||
type(spline_simple) :: beta, waist1, waist2, &
|
||||
rci1, rci2, phi1, phi2, &
|
||||
x0, y0, z0
|
||||
|
||||
open(newunit=u, file=params%filenm, status='old', action='read', iostat=err)
|
||||
if (err /= 0) then
|
||||
@ -145,47 +143,95 @@ contains
|
||||
rci1v = 10 * rci1v
|
||||
rci2v = 10 * rci2v
|
||||
|
||||
call beta%init(alphastv, betastv, nisteer)
|
||||
call waist1%init(alphastv, waist1v, nisteer)
|
||||
call waist2%init(alphastv, waist2v, nisteer)
|
||||
call rci1%init(alphastv, rci1v, nisteer)
|
||||
call rci2%init(alphastv, rci2v, nisteer)
|
||||
call phi1%init(alphastv, phi1v, nisteer)
|
||||
call phi2%init(alphastv, phi2v, nisteer)
|
||||
call x0%init(alphastv, x00v, nisteer)
|
||||
call y0%init(alphastv, y00v, nisteer)
|
||||
call z0%init(alphastv, z00v, nisteer)
|
||||
call locate(alphastv, params%alpha, k)
|
||||
|
||||
if((params%alpha > alphastv(1)) .and. (params%alpha < alphastv(nisteer))) then
|
||||
call locate(alphastv, params%alpha, k)
|
||||
dal = params%alpha - alphastv(k)
|
||||
params%beta = beta%raw_eval(k, dal)
|
||||
params%pos(1) = x0%raw_eval(k, dal)
|
||||
params%pos(2) = y0%raw_eval(k, dal)
|
||||
params%pos(3) = z0%raw_eval(k, dal)
|
||||
params%w(1) = waist1%raw_eval(k, dal)
|
||||
params%w(2) = waist2%raw_eval(k, dal)
|
||||
params%ri(1) = rci1%raw_eval(k, dal)
|
||||
params%ri(2) = rci2%raw_eval(k, dal)
|
||||
params%phi(1) = phi1%raw_eval(k, dal)
|
||||
params%phi(2) = phi2%raw_eval(k, dal)
|
||||
else
|
||||
! params%alpha outside table range
|
||||
if(params%alpha >= alphastv(nisteer)) ii=nisteer
|
||||
if(params%alpha <= alphastv(1)) ii=1
|
||||
params%alpha = alphastv(ii)
|
||||
params%beta = betastv(ii)
|
||||
params%pos(1) = x00v(ii)
|
||||
params%pos(2) = y00v(ii)
|
||||
params%pos(3) = z00v(ii)
|
||||
params%w(1) = waist1v(ii)
|
||||
params%w(2) = waist2v(ii)
|
||||
params%ri(1) = rci1v(ii)
|
||||
params%ri(2) = rci2v(ii)
|
||||
params%phi(1) = phi1v(ii)
|
||||
params%phi(2) = phi2v(ii)
|
||||
! No interpolation: 0D table or out of range
|
||||
if (nisteer == 1 .or. k == 0 .or. k == nisteer) then
|
||||
if (k == 0) k = k + 1
|
||||
params%alpha = alphastv(k)
|
||||
params%beta = betastv(k)
|
||||
params%pos(1) = x00v(k)
|
||||
params%pos(2) = y00v(k)
|
||||
params%pos(3) = z00v(k)
|
||||
params%w(1) = waist1v(k)
|
||||
params%w(2) = waist2v(k)
|
||||
params%ri(1) = rci1v(k)
|
||||
params%ri(2) = rci2v(k)
|
||||
params%phi(1) = phi1v(k)
|
||||
params%phi(2) = phi2v(k)
|
||||
return
|
||||
end if
|
||||
|
||||
if (nisteer < 4) then
|
||||
! Linear interpolation for small tables
|
||||
block
|
||||
type(linear_1d) :: beta, waist1, waist2, &
|
||||
rci1, rci2, phi1, phi2, &
|
||||
x0, y0, z0
|
||||
|
||||
call beta%init(alphastv, betastv, nisteer)
|
||||
call waist1%init(alphastv, waist1v, nisteer)
|
||||
call waist2%init(alphastv, waist2v, nisteer)
|
||||
call rci1%init(alphastv, rci1v, nisteer)
|
||||
call rci2%init(alphastv, rci2v, nisteer)
|
||||
call phi1%init(alphastv, phi1v, nisteer)
|
||||
call phi2%init(alphastv, phi2v, nisteer)
|
||||
call x0%init(alphastv, x00v, nisteer)
|
||||
call y0%init(alphastv, y00v, nisteer)
|
||||
call z0%init(alphastv, z00v, nisteer)
|
||||
|
||||
params%beta = beta%raw_eval(k, params%alpha)
|
||||
params%pos(1) = x0%raw_eval(k, params%alpha)
|
||||
params%pos(2) = y0%raw_eval(k, params%alpha)
|
||||
params%pos(3) = z0%raw_eval(k, params%alpha)
|
||||
params%w(1) = waist1%raw_eval(k, params%alpha)
|
||||
params%w(2) = waist2%raw_eval(k, params%alpha)
|
||||
params%ri(1) = rci1%raw_eval(k, params%alpha)
|
||||
params%ri(2) = rci2%raw_eval(k, params%alpha)
|
||||
params%phi(1) = phi1%raw_eval(k, params%alpha)
|
||||
params%phi(2) = phi2%raw_eval(k, params%alpha)
|
||||
return
|
||||
end block
|
||||
end if
|
||||
|
||||
! Cubic interpolation otherwise
|
||||
block
|
||||
type(spline_1d) :: beta, waist1, waist2, &
|
||||
rci1, rci2, phi1, phi2, &
|
||||
x0, y0, z0
|
||||
|
||||
call beta%init(alphastv, betastv, nisteer, err=err)
|
||||
call waist1%init(alphastv, waist1v, nisteer, err=err)
|
||||
call waist2%init(alphastv, waist2v, nisteer, err=err)
|
||||
call rci1%init(alphastv, rci1v, nisteer, err=err)
|
||||
call rci2%init(alphastv, rci2v, nisteer, err=err)
|
||||
call phi1%init(alphastv, phi1v, nisteer, err=err)
|
||||
call phi2%init(alphastv, phi2v, nisteer, err=err)
|
||||
call x0%init(alphastv, x00v, nisteer, err=err)
|
||||
call y0%init(alphastv, y00v, nisteer, err=err)
|
||||
call z0%init(alphastv, z00v, nisteer, err=err)
|
||||
|
||||
if (err > 0) then
|
||||
write (msg, '(a, a, g0)') 'interpolation failed', &
|
||||
'`curfit` returned ier=', err
|
||||
call log_error(msg, mod='beams', proc='read_beam1')
|
||||
return
|
||||
end if
|
||||
err = 0
|
||||
|
||||
call locate(beta%knots, params%alpha, k)
|
||||
params%beta = beta%raw_eval(k, params%alpha)
|
||||
params%pos(1) = x0%raw_eval(k, params%alpha)
|
||||
params%pos(2) = y0%raw_eval(k, params%alpha)
|
||||
params%pos(3) = z0%raw_eval(k, params%alpha)
|
||||
params%w(1) = waist1%raw_eval(k, params%alpha)
|
||||
params%w(2) = waist2%raw_eval(k, params%alpha)
|
||||
params%ri(1) = rci1%raw_eval(k, params%alpha)
|
||||
params%ri(2) = rci2%raw_eval(k, params%alpha)
|
||||
params%phi(1) = phi1%raw_eval(k, params%alpha)
|
||||
params%phi(2) = phi2%raw_eval(k, params%alpha)
|
||||
end block
|
||||
|
||||
end subroutine read_beam1
|
||||
|
||||
|
||||
|
||||
@ -4,7 +4,7 @@
|
||||
!
|
||||
module gray_equil
|
||||
use const_and_precisions, only : wp_, comp_huge
|
||||
use splines, only : spline_simple, spline_1d, spline_2d, linear_1d
|
||||
use splines, only : spline_1d, spline_2d, linear_1d
|
||||
use types, only : contour
|
||||
|
||||
implicit none
|
||||
@ -26,12 +26,12 @@ module gray_equil
|
||||
real(wp_) :: z_boundary(2) ! z range of the plasma boundary
|
||||
|
||||
! Flux average splines (see `flux_average`)
|
||||
type(spline_simple) :: spline_area, spline_volume
|
||||
type(spline_simple) :: spline_dAdrho_t, spline_dVdrho_t
|
||||
type(spline_simple) :: spline_R_J, spline_B_avg, spline_B_min, spline_B_max
|
||||
type(spline_simple) :: spline_f_c, spline_q, spline_I_p, spline_J_phi_avg
|
||||
type(spline_simple) :: spline_astra_tor, spline_jintrac_tor, spline_paral_tor
|
||||
type(spline_2d) :: spline_cd_eff
|
||||
type(spline_1d) :: spline_area, spline_volume
|
||||
type(spline_1d) :: spline_dAdrho_t, spline_dVdrho_t
|
||||
type(spline_1d) :: spline_R_J, spline_B_avg, spline_B_min, spline_B_max
|
||||
type(spline_1d) :: spline_f_c, spline_q, spline_I_p, spline_J_phi_avg
|
||||
type(spline_1d) :: spline_astra_tor, spline_jintrac_tor, spline_paral_tor
|
||||
type(spline_2d) :: spline_cd_eff
|
||||
|
||||
contains
|
||||
procedure(pol_flux_sub), deferred :: pol_flux
|
||||
@ -147,11 +147,10 @@ module gray_equil
|
||||
private
|
||||
real(wp_) :: fpol_a ! Poloidal current at the edge (r=a), F_a
|
||||
! Splines
|
||||
type(spline_2d) :: psi_spline
|
||||
type(contour) :: psi_domain
|
||||
type(spline_1d) :: fpol_spline
|
||||
type(spline_simple) :: q_spline
|
||||
type(linear_1d) :: rhop_spline, rhot_spline
|
||||
type(spline_2d) :: psi_spline
|
||||
type(contour) :: psi_domain
|
||||
type(spline_1d) :: fpol_spline, q_spline
|
||||
type(linear_1d) :: rhop_spline, rhot_spline
|
||||
contains
|
||||
procedure :: pol_flux => numeric_pol_flux
|
||||
procedure :: pol_curr => numeric_pol_curr
|
||||
@ -773,6 +772,7 @@ contains
|
||||
use gray_params, only : EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
|
||||
use gray_params, only : X_AT_TOP, X_AT_BOTTOM, X_IS_MISSING
|
||||
use utils, only : vmaxmin
|
||||
use numint, only : simpson
|
||||
use logger, only : log_debug, log_info
|
||||
|
||||
! subroutine arguments
|
||||
@ -1017,13 +1017,11 @@ contains
|
||||
|
||||
call self%q_spline%init(data%psi, data%q, npsi)
|
||||
|
||||
! Toroidal flux as Φ(ψ) = 2π ∫ q(ψ)dψ
|
||||
! Specifically, Φ(ψ_n) = 2π|ψ_a| ∫₀ψ_n q(ψ_n)dψ_n, with ψ_n ∈ [0,1]
|
||||
! Toroidal flux as Φ(ψ) = 2π ∫₀^ψ q(ψ) dψ
|
||||
phi(1) = 0
|
||||
do k = 1, npsi-1
|
||||
dx = self%q_spline%data(k+1) - self%q_spline%data(k)
|
||||
phi(k+1) = phi(k) + dx*(self%q_spline%coeffs(k,1) + dx*(self%q_spline%coeffs(k,2)/2 + &
|
||||
dx*(self%q_spline%coeffs(k,3)/3 + dx* self%q_spline%coeffs(k,4)/4) ) )
|
||||
do k = 2, npsi
|
||||
call simpson(2, h=data%psi(2)-data%psi(1), fi=data%q(k-1:k), s=phi(k))
|
||||
phi(k) = phi(k) + phi(k-1)
|
||||
end do
|
||||
|
||||
self%phi_a = phi(npsi)
|
||||
|
||||
@ -8,7 +8,7 @@
|
||||
!
|
||||
module gray_plasma
|
||||
use const_and_precisions, only : wp_, zero
|
||||
use splines, only : spline_simple, spline_1d
|
||||
use splines, only : spline_1d
|
||||
|
||||
implicit none
|
||||
|
||||
@ -79,10 +79,10 @@ module gray_plasma
|
||||
! Numerical plasma description
|
||||
type, extends(abstract_plasma) :: numeric_plasma
|
||||
private
|
||||
type(spline_1d) :: dens_spline
|
||||
type(spline_simple) :: temp_spline
|
||||
type(spline_simple) :: zeff_spline
|
||||
type(density_tail) :: tail
|
||||
type(spline_1d) :: dens_spline
|
||||
type(spline_1d) :: temp_spline
|
||||
type(spline_1d) :: zeff_spline
|
||||
type(density_tail) :: tail
|
||||
contains
|
||||
procedure :: init => init_splines
|
||||
procedure :: density => numeric_density
|
||||
|
||||
228
src/splines.f90
228
src/splines.f90
@ -1,7 +1,6 @@
|
||||
! This module provides a high-level interface for creating and evaluating
|
||||
! several kind of splines:
|
||||
!
|
||||
! `spline_simple` is a simple interpolating cubic spline,
|
||||
! `spline_1d` and `spline_2d` are wrapper around the DIERCKX cubic splines,
|
||||
! `linear_1d` is a linear interpolation.
|
||||
module splines
|
||||
@ -9,27 +8,16 @@ module splines
|
||||
|
||||
implicit none
|
||||
|
||||
! A 1D interpolating cubic spline
|
||||
type spline_simple
|
||||
integer :: ndata ! Number of data points
|
||||
real(wp_), allocatable :: data(:) ! Data points (ndata)
|
||||
real(wp_), allocatable :: coeffs(:,:) ! Spline coefficients (ndata, 4)
|
||||
contains
|
||||
procedure :: init => spline_simple_init
|
||||
procedure :: eval => spline_simple_eval
|
||||
procedure :: raw_eval => spline_simple_raw_eval
|
||||
procedure :: deriv => spline_simple_deriv
|
||||
end type
|
||||
|
||||
! A 1D smoothing/interpolating cubic spline
|
||||
type spline_1d
|
||||
integer :: nknots ! Number of spline knots
|
||||
real(wp_), allocatable :: knots(:) ! Knots positions
|
||||
real(wp_), allocatable :: coeffs(:) ! B-spline coefficients
|
||||
contains
|
||||
procedure :: init => spline_1d_init
|
||||
procedure :: eval => spline_1d_eval
|
||||
procedure :: deriv => spline_1d_deriv
|
||||
procedure :: init => spline_1d_init
|
||||
procedure :: eval => spline_1d_eval
|
||||
procedure :: raw_eval => spline_1d_raw_eval
|
||||
procedure :: deriv => spline_1d_deriv
|
||||
end type
|
||||
|
||||
! Wrapper to store pointers in an array
|
||||
@ -70,178 +58,10 @@ module splines
|
||||
end type
|
||||
|
||||
private
|
||||
public spline_simple, spline_1d, spline_2d, linear_1d
|
||||
public spline_1d, spline_2d, linear_1d
|
||||
|
||||
contains
|
||||
|
||||
subroutine spline_simple_init(self, x, y, n)
|
||||
! Initialises the spline
|
||||
|
||||
! subroutine arguments
|
||||
class(spline_simple), intent(inout) :: self
|
||||
integer, intent(in) :: n
|
||||
real(wp_), dimension(n), intent(in) :: x, y
|
||||
|
||||
self%ndata = n
|
||||
allocate(self%data(n))
|
||||
allocate(self%coeffs(n, 4))
|
||||
|
||||
self%data = x
|
||||
call spline_simple_coeffs(x, y, n, self%coeffs)
|
||||
end subroutine spline_simple_init
|
||||
|
||||
|
||||
pure function spline_simple_eval(self, x) result(y)
|
||||
! Evaluates the spline at x
|
||||
use utils, only : locate
|
||||
|
||||
! subroutine arguments
|
||||
class(spline_simple), intent(in) :: self
|
||||
real(wp_), intent(in) :: x
|
||||
real(wp_) :: y
|
||||
|
||||
! local variables
|
||||
integer :: i
|
||||
real(wp_) :: dx
|
||||
|
||||
call locate(self%data, x, i)
|
||||
i = min(max(1, i), self%ndata - 1)
|
||||
dx = x - self%data(i)
|
||||
y = self%raw_eval(i, dx)
|
||||
end function spline_simple_eval
|
||||
|
||||
|
||||
pure function spline_simple_raw_eval(self, i, dx) result(y)
|
||||
! Evaluates the i-th polynomial of the spline at dx
|
||||
|
||||
! subroutine arguments
|
||||
class(spline_simple), intent(in) :: self
|
||||
integer, intent(in) :: i
|
||||
real(wp_), intent(in) :: dx
|
||||
real(wp_) :: y
|
||||
|
||||
y = self%coeffs(i,1) + dx*(self%coeffs(i,2) &
|
||||
+ dx*(self%coeffs(i,3) + dx*self%coeffs(i,4)))
|
||||
end function spline_simple_raw_eval
|
||||
|
||||
|
||||
pure function spline_simple_deriv(self, x) result(y)
|
||||
! Computes the derivative of the spline at x
|
||||
use utils, only : locate
|
||||
|
||||
|
||||
! subroutine arguments
|
||||
class(spline_simple), intent(in) :: self
|
||||
real(wp_), intent(in) :: x
|
||||
real(wp_) :: y
|
||||
|
||||
! local variables
|
||||
integer :: i
|
||||
real(wp_) :: dx
|
||||
|
||||
call locate(self%data, x, i)
|
||||
i = min(max(1, i), self%ndata - 1)
|
||||
dx = x - self%data(i)
|
||||
y = self%coeffs(i,2) + dx*(2*self%coeffs(i,3) + 3*dx*self%coeffs(i,4))
|
||||
end function spline_simple_deriv
|
||||
|
||||
|
||||
pure subroutine spline_simple_coeffs(x, y, n, c)
|
||||
! Computes the cubic coefficients of all n polynomials
|
||||
|
||||
! subroutine arguments
|
||||
integer, intent(in) :: n
|
||||
real(wp_), intent(in) :: x(n), y(n)
|
||||
real(wp_), intent(inout) :: c(n*4)
|
||||
|
||||
! local variables
|
||||
integer :: jmp, i, ii, j, j1, j2, j3, k
|
||||
real(wp_) :: xb, xc, ya, yb, h, a, r, dya, dyb, dy2
|
||||
jmp = 1
|
||||
|
||||
if (n <= 1) return
|
||||
|
||||
! initialisation
|
||||
xc = x(1)
|
||||
yb = y(1)
|
||||
h = 0
|
||||
a = 0
|
||||
r = 0
|
||||
dyb = 0
|
||||
dy2 = c(2)
|
||||
|
||||
! form the system of linear equations
|
||||
! and eliminate subsequentially
|
||||
j = 1
|
||||
do i = 1, n
|
||||
j2 = n + i
|
||||
j3 = j2 + n
|
||||
a = h*(2 - a)
|
||||
dya = dyb + h*r
|
||||
if (i >= n) then
|
||||
! set derivative dy2 at last point
|
||||
dyb = dy2
|
||||
h = 0
|
||||
dyb = dya
|
||||
goto 13
|
||||
else
|
||||
j = j+jmp
|
||||
xb = xc
|
||||
xc = x(j)
|
||||
h = xc-xb
|
||||
|
||||
! ii=0 - increasing abscissae
|
||||
! ii=1 - decreasing abscissae
|
||||
ii = 0
|
||||
if (h == 0) return
|
||||
if (h < 0) ii = 1
|
||||
ya = yb
|
||||
yb = y(j)
|
||||
dyb = (yb - ya)/h
|
||||
|
||||
if (i <= 1) then
|
||||
j1 = ii
|
||||
goto 13
|
||||
end if
|
||||
end if
|
||||
|
||||
if (j1-ii /= 0) return
|
||||
a = 1 / (2*h + a)
|
||||
13 continue
|
||||
r = a*(dyb - dya)
|
||||
c(j3) = r
|
||||
a = h*a
|
||||
c(j2) = a
|
||||
c(i) = dyb
|
||||
end do
|
||||
|
||||
! back substitution of the system of linear equations
|
||||
! and computation of the other coefficients
|
||||
a = 1
|
||||
j1 = j3+n+ii-ii*n
|
||||
i = n
|
||||
do k = 1, n
|
||||
xb = x(j)
|
||||
h = xc - xb
|
||||
xc = xb
|
||||
a = a+h
|
||||
yb = r
|
||||
r = c(j3)-r*c(j2)
|
||||
ya = 2*r
|
||||
c(j3) = ya + r
|
||||
c(j2) = c(i) - h*(ya+yb)
|
||||
c(j1) = (yb - r)/a
|
||||
c(i) = y(j)
|
||||
a = 0
|
||||
j = j-jmp
|
||||
i = i-1
|
||||
j2 = j2-1
|
||||
j3 = j3-1
|
||||
j1 = j3+n+ii
|
||||
end do
|
||||
end subroutine spline_simple_coeffs
|
||||
|
||||
|
||||
subroutine spline_1d_init(self, x, y, n, range, weights, tension, err)
|
||||
! Initialises a spline_1d.
|
||||
! Takes:
|
||||
@ -260,7 +80,7 @@ contains
|
||||
real(wp_), intent(in) :: x(n)
|
||||
real(wp_), intent(in) :: y(n)
|
||||
integer, intent(in) :: n
|
||||
real(wp_), intent(in) :: range(2)
|
||||
real(wp_), intent(in), optional :: range(2)
|
||||
real(wp_), intent(in), optional :: weights(n)
|
||||
real(wp_), intent(in), optional :: tension
|
||||
integer, intent(out), optional :: err
|
||||
@ -273,10 +93,12 @@ contains
|
||||
|
||||
! default values
|
||||
integer :: err_def
|
||||
real(wp_) :: weights_def(n), tension_def
|
||||
real(wp_) :: range_def(2), weights_def(n), tension_def
|
||||
|
||||
range_def = [x(1), x(n)]
|
||||
weights_def = 1
|
||||
tension_def = 0
|
||||
if (present(range)) range_def = range
|
||||
if (present(weights)) weights_def = weights
|
||||
if (present(tension)) tension_def = tension
|
||||
|
||||
@ -286,9 +108,10 @@ contains
|
||||
allocate(self%knots(nknots_est), self%coeffs(nknots_est))
|
||||
end if
|
||||
|
||||
call curfit(0, n, x, y, weights_def, range(1), range(2), 3, tension_def, &
|
||||
nknots_est, self%nknots, self%knots, self%coeffs, residuals, &
|
||||
work_real, size(work_real), work_int, err_def)
|
||||
call curfit(0, n, x, y, weights_def, range_def(1), range_def(2), &
|
||||
3, tension_def, nknots_est, self%nknots, self%knots, &
|
||||
self%coeffs, residuals, work_real, size(work_real), &
|
||||
work_int, err_def)
|
||||
|
||||
if (present(err)) err = err_def
|
||||
end subroutine spline_1d_init
|
||||
@ -312,6 +135,31 @@ contains
|
||||
end function spline_1d_eval
|
||||
|
||||
|
||||
pure function spline_1d_raw_eval(self, i, x) result(y)
|
||||
! Evaluates the spline at: t(i) ≤ x < t(i+1), where t are the knots
|
||||
use dierckx, only : fpbspl
|
||||
|
||||
! subroutine arguments
|
||||
class(spline_1d), intent(in) :: self
|
||||
integer, intent(in) :: i
|
||||
real(wp_), intent(in) :: x
|
||||
real(wp_) :: y
|
||||
|
||||
! local variables
|
||||
integer :: err, j
|
||||
real(wp_) :: h(6)
|
||||
|
||||
! Evaluate the basis splines
|
||||
call fpbspl(self%knots, self%nknots, 3, x, i, h)
|
||||
|
||||
! Take the linear combination with the coefficients
|
||||
y = 0
|
||||
do j = 1, 4
|
||||
y = y + self%coeffs(i-4+j)*h(j)
|
||||
end do
|
||||
end function spline_1d_raw_eval
|
||||
|
||||
|
||||
pure function spline_1d_deriv(self, x, order) result(y)
|
||||
! Evaluates the spline n-th order derivative at x
|
||||
use dierckx, only : splder
|
||||
|
||||
Loading…
Reference in New Issue
Block a user