src/splines.f90: add mollifier_1d

src/splines.f90

@@ -3,6 +3,8 @@
 !
 !   `spline_1d` and `spline_2d` are wrapper around the DIERCKX cubic splines,
 !   `linear_1d` is a linear interpolation.
+!   `mollifier_1d` is a smooth interpolation scheme based on a C² mollifier
+
 module splines
   use const_and_precisions, only : wp_
 
@@ -57,8 +59,19 @@ module splines
       procedure :: raw_eval => linear_1d_raw_eval
   end type
 
+  ! A 1D smooth interpolation based on a C² mollifier
+  type mollifier_1d
+    integer                :: ndata    ! Number of data points
+    real(wp_), allocatable :: xdata(:) ! X data
+    real(wp_), allocatable :: ydata(:) ! Y data
+    real(wp_)              :: w        ! kernel width
+    contains
+      procedure :: init => mollifier_1d_init
+      procedure :: eval => mollifier_1d_eval
+  end type
+
   private
-  public spline_1d, spline_2d, linear_1d
+  public spline_1d, spline_2d, linear_1d, mollifier_1d
 
 contains
 
@@ -495,4 +508,137 @@ contains
     y = self%ydata(i) + dy/dx * (x - self%xdata(i))
   end function linear_1d_raw_eval
 
+
+  subroutine mollifier_1d_init(self, x, y, width)
+    ! Initialises the mollifier
+
+    ! subroutine arguments
+    class(mollifier_1d),     intent(inout) :: self
+    real(wp_), dimension(:), intent(in)    :: x, y
+    real(wp_),               intent(in)    :: width
+
+    self%ndata = size(x)
+    self%xdata = x
+    self%ydata = y
+    self%w = width
+  end subroutine mollifier_1d_init
+
+
+  pure function mollifier_1d_eval(self, x, order) result(y)
+    ! Evaluates the interpolated data, or the derivative
+    ! of given `order`, at x
+    use utils, only : locate
+
+    ! subroutine arguments
+    class(mollifier_1d), intent(in) :: self
+    real(wp_),           intent(in) :: x
+    integer, optional,   intent(in) :: order
+    real(wp_)                       :: y
+
+    ! local variables
+    real(wp_) :: m, m1, q, a, b
+    real(wp_) :: J, K
+    integer :: i, left, right
+    integer :: order_
+
+    order_ = 0
+    if (present(order)) order_ = order
+
+    ! We compute the convolution ∫ f(x)⋅g((x₀-x)/w)/w dx
+    ! where f is a linear interpolation of the data (x, y),
+    !       g(x) = 70⋅(1/2 + x/2)³(1/2 - x/2)³ χ_[‐1,+1](x)
+    !
+    ! The integral can be rewritten as:
+    !   I(x₀) = ∫₋₁⁺¹ f(x₀-tw)⋅g(t) dt
+    !         = Σ_i ∫_aᵢ^bᵢ (mᵢ(x₀-tw) + qᵢ)⋅g(t) dt
+    ! where the sum is over the samples within [x₀-w, x₀+w] with
+    !   aᵢ = max((x₀-xᵢ)/w, -1)
+    !   bᵢ = min((x₀-xᵢ₋₁])/w, +1)
+    !
+    ! The convolution reduces to only 2 integrals:
+    !   I(x₀) = Σ_i { (mᵢx₀ + qᵢ)⋅Jᵢ - mᵢw Kᵢ}
+    ! where:
+    !   Jᵢ = ∫_aᵢ^bᵢ g(t) dt,
+    !   Kᵢ = ∫_aᵢ^bᵢ t⋅g(t) dt.
+
+    ! Find data within [x₀-w, x₀+w]
+    call locate(self%xdata, x - self%w, left)
+    call locate(self%xdata, x + self%w, right)
+
+    ! if x(left) = x₀-w increase left by one to get
+    ! the line segment actually within the window.
+    if (x - self%w == self%xdata(left)) left = left + 1
+
+    y = 0
+    do i = left+1, right+1
+      ! Compute f
+      if (i <= 1) then
+        ! f=y(1) for x≤x(1)
+        m = 0
+        q = self%ydata(1)
+        a = max((x-self%xdata(1))/self%w, -1.0)
+        b = 1
+      else if (i == self%ndata + 1) then
+        ! f=y(n) for x≥x(n)
+        if (order_ == 2) exit
+        m = 0
+        q = self%ydata(self%ndata)
+        a = -1
+        b = min((x-self%xdata(self%ndata))/self%w, +1.0)
+      else
+        ! linear interpolation: f(x) = mx + q
+        m = (self%ydata(i) - self%ydata(i-1)) / (self%xdata(i) - self%xdata(i-1))
+        q = self%ydata(i-1) - m * self%xdata(i-1)
+        a = max((x-self%xdata(i))/self%w, -1.0)
+        b = min((x-self%xdata(i-1))/self%w, +1.0)
+      end if
+
+      ! Convolve with mollifier
+      select case (order_)
+        case (0) ! Compute f*g
+          J = g_integral(b) - g_integral(a)
+          K = g_x_integral(b) - g_x_integral(a)
+          y = y + (m*x + q)*J - (m*self%w)*K
+
+        case (1) ! Compute f'*g
+          J = g_integral(b) - g_integral(a)
+          y = y + m*J
+
+        case (2) ! Compute f"*g
+          m1 = 0
+          if (i < size(self%xdata)) &
+            m1 = (self%ydata(i+1) - self%ydata(i)) / (self%xdata(i+1) - self%xdata(i))
+          y = y + g((x-self%xdata(i))/self%w)/self%w * (m1 - m)
+      end select
+    end do
+
+  contains
+
+    pure function g(x)
+      ! A C²_c([-1,1]) polynomial mollifier
+      real(wp_), intent(in) :: x
+      real(wp_)             :: g
+      if (abs(x) >= 1) then
+        g = 0
+      else
+        g = 35 * (1 + x)**3 * (1 - x)**3 / 32
+      end if
+    end function
+
+    pure function g_integral(x)
+      ! Returns ∫₀^x g(t)dt
+      real(wp_), intent(in) :: x
+      real(wp_)             :: g_integral
+      g_integral = x*(x**2*(x**2*(21 - 5*x**2) - 35) + 35)/32
+    end function
+
+    pure function g_x_integral(x)
+      ! Returns ∫₀^x t⋅g(t)dt
+      real(wp_), intent(in) :: x
+      real(wp_)             :: g_x_integral
+      g_x_integral = -35*x**2*(x**2 - 2)*(x**2*(x**2 - 2) + 2)/256
+    end function
+
+  end function mollifier_1d_eval
+
 end module splines