From daf3d500af644d3b2a8432c18e993672d80f380e Mon Sep 17 00:00:00 2001
From: Michele Guerini Rocco <m.guerinirocco2@campus.unimib.it>
Date: Wed, 19 Apr 2023 14:32:32 +0200
Subject: [PATCH] use linear interpolation for monotonic data
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The ρ_p/ρ_t mapping is 1:1, so the interpolation must always preserve
monotonicity, of which cubic splines generally make no guarantee.

Note: Linear interpolation does not provide even C¹ continuity, but
these data is not directly used in the numerical integration, so it
should be fine. Ideally this should be replaced with cubic splines
computed with the Fritsch–Carlson algorithm.
---
 src/equilibrium.f90 | 10 ++---
 src/splines.f90     | 90 ++++++++++++++++++++++++++++++++++++++++++++-
 2 files changed, 92 insertions(+), 8 deletions(-)

diff --git a/src/equilibrium.f90 b/src/equilibrium.f90
index 7feb2b9..c53e186 100644
--- a/src/equilibrium.f90
+++ b/src/equilibrium.f90
@@ -7,7 +7,7 @@
 ! the data is interpolated using splines.
 module equilibrium
   use const_and_precisions, only : wp_
-  use splines,              only : spline_simple, spline_1d, spline_2d
+  use splines,              only : spline_simple, spline_1d, spline_2d, linear_1d
 
   implicit none
 
@@ -27,7 +27,7 @@ module equilibrium
   type(spline_1d),     save :: fpol_spline              ! Poloidal current function F(ψ)
   type(spline_2d),     save :: psi_spline               ! Poloidal flux ψ(R,z)
   type(spline_simple), save :: q_spline                 ! Safey factor q(ψ)
-  type(spline_simple), save :: rhop_spline, rhot_spline ! Normalised radii ρ_p(ρ_t), ρ_t(ρ_p)
+  type(linear_1d),     save :: rhop_spline, rhot_spline ! Normalised radii ρ_p(ρ_t), ρ_t(ρ_p)
 
   ! Analytic model
   type(analytic_model), save :: model
@@ -951,14 +951,12 @@ contains
 
     ! local variables
     integer :: i, i0, j
-    real(wp_) :: dr
 
     i0 = 1
     do j = 1, size(rhot)
-      call locate(rhop_spline%data(i0:), rhop_spline%ndata-i0+1, rhot(j), i)
+      call locate(rhop_spline%xdata(i0:), rhop_spline%ndata-i0+1, rhot(j), i)
       i = min(max(1,i), rhop_spline%ndata - i0) + i0 - 1
-      dr = rhot(j) - rhop_spline%data(i)
-      frhopolv(j) = rhop_spline%raw_eval(i, dr)
+      frhopolv(j) = rhop_spline%raw_eval(i, rhot(j))
       i0 = i
     end do
   end function frhopolv
diff --git a/src/splines.f90 b/src/splines.f90
index 7174323..097200c 100644
--- a/src/splines.f90
+++ b/src/splines.f90
@@ -2,7 +2,8 @@
 ! several kind of splines:
 !
 !   `spline_simple` is a simple interpolating cubic spline,
-!   `spline_1d` and `spline_2d` are wrapper around the DIERCKX cubic splines.
+!   `spline_1d` and `spline_2d` are wrapper around the DIERCKX cubic splines,
+!   `linear_1d` is a linear interpolation.
 module splines
   use const_and_precisions, only : wp_
 
@@ -57,8 +58,20 @@ module splines
     real(wp_), pointer :: ptr(:) => null()
   end type
 
+  ! A simple 1D linear interpolation
+  type linear_1d
+    integer :: ndata                   ! Number of data points
+    real(wp_), allocatable :: xdata(:) ! X data (ndata)
+    real(wp_), allocatable :: ydata(:) ! Y data (ndata)
+    contains
+      procedure :: init     => linear_1d_init
+      procedure :: deinit   => linear_1d_deinit
+      procedure :: eval     => linear_1d_eval
+      procedure :: raw_eval => linear_1d_raw_eval
+  end type
+
   private
-  public spline_simple, spline_1d, spline_2d
+  public spline_simple, spline_1d, spline_2d, linear_1d
 
 contains
 
@@ -543,4 +556,77 @@ contains
     z = zv(1)
   end function spline_2d_deriv
 
+
+  subroutine linear_1d_init(self, x, y, n)
+    ! Initialises the spline
+
+    implicit none
+
+    ! subroutine arguments
+    class(linear_1d),        intent(inout) :: self
+    integer,                 intent(in)    :: n
+    real(wp_), dimension(n), intent(in)    :: x, y
+
+    call self%deinit
+    self%ndata = n
+    allocate(self%xdata(n))
+    allocate(self%ydata(n))
+
+    self%xdata = x
+    self%ydata = y
+  end subroutine linear_1d_init
+
+
+  subroutine linear_1d_deinit(self)
+    ! Deinitialises a linear_1d
+    implicit none
+
+    ! subroutine arguments
+    class(linear_1d), intent(inout) :: self
+
+    if (allocated(self%xdata)) deallocate(self%xdata)
+    if (allocated(self%ydata)) deallocate(self%ydata)
+    self%ndata = 0
+  end subroutine linear_1d_deinit
+
+
+  function linear_1d_eval(self, x) result(y)
+    ! Evaluates the linear interpolated data at x
+    use utils, only : locate
+
+    implicit none
+
+    ! subroutine arguments
+    class(linear_1d), intent(in) :: self
+    real(wp_),        intent(in) :: x
+    real(wp_)                    :: y
+
+    ! local variables
+    integer :: i
+
+    call locate(self%xdata, self%ndata, x, i)
+    i = min(max(1, i), self%ndata - 1)
+    y = self%raw_eval(i, x)
+  end function linear_1d_eval
+
+
+  function linear_1d_raw_eval(self, i, x) result(y)
+    ! Performs the linear interpolation in the i-th interval
+
+    implicit none
+
+    ! subroutine arguments
+    class(linear_1d), intent(in) :: self
+    integer,          intent(in) :: i
+    real(wp_),        intent(in) :: x
+    real(wp_)                    :: y
+
+    ! local variables
+    real(wp_) :: dx, dy
+
+    dx = self%xdata(i+1) - self%xdata(i)
+    dy = self%ydata(i+1) - self%ydata(i)
+    y = self%ydata(i) + dy/dx * (x - self%xdata(i))
+  end function linear_1d_raw_eval
+
 end module splines