Started interpolate module + tests

CMakeSources.in

@@ -12,6 +12,7 @@ set(SOURCE_FILES
     src/solve_systems.f90
     src/system_utility.f90
     src/util.f90
+    src/interpolate.f90
 
     # VMEC
     src/canonical_coordinates_mod.f90

src/interpolate.f90

@@ -0,0 +1,73 @@
+module interpolate
+    implicit none
+
+contains
+
+    subroutine construct_splines_3d(spl_data, spl_order, n_r, n_z, n_phi, h_r, h_z, h_phi)
+        integer, intent(in) :: spl_order
+        integer, intent(in) :: n_r, n_z, n_phi
+        real(8), intent(in) :: h_r, h_z, h_phi
+        real(8), intent(inout) :: spl_data(:,:,:,:,:,:,:)
+
+        real(8), dimension(:,:), allocatable  :: splcoe
+
+        integer :: n_data
+        integer :: i_r, i_z, i_phi            ! Loop counters for grid
+        integer :: i_r_z, i_r_phi, i_data, k  ! Loop counters for splines
+
+        n_data = size(spl_data, 1)
+
+        ! Spline over $\varphi$
+        allocate(splcoe(0:spl_order,n_phi))
+        do i_r=1,n_r
+        do i_z=1,n_z
+            do i_data=1,n_data
+                splcoe(0,:)=spl_data(i_data,1,1,1,i_r,i_z,:)
+                call spl_per(spl_order,n_phi,h_phi,splcoe)
+                do k=1,spl_order
+                    spl_data(i_data,1,1,k+1,i_r,i_z,:)=splcoe(k,:)
+                enddo
+            enddo
+        enddo
+        enddo
+        deallocate(splcoe)
+
+        ! Spline over $Z$
+        allocate(splcoe(0:spl_order,n_z))
+        do i_r=1,n_r
+        do i_phi=1,n_phi
+            do i_r_phi=1,spl_order+1
+                do i_data=1,n_data
+                    splcoe(0,:)=spl_data(i_data,1,1,i_r_phi,i_r,:,i_phi)
+                    call spl_reg(spl_order,n_z,h_z,splcoe)
+                    do k=1,spl_order
+                        spl_data(i_data,1,k+1,i_r_phi,i_r,:,i_phi)=splcoe(k,:)
+                    enddo
+                enddo
+            enddo
+        enddo
+        enddo
+        deallocate(splcoe)
+
+        ! Spline over $R$
+        allocate(splcoe(0:spl_order,n_r))
+        do i_z=1,n_z
+        do i_phi=1,n_phi
+            do i_r_z=1,spl_order+1
+            do i_r_phi=1,spl_order+1
+                do i_data=1,n_data
+                    splcoe(0,:)=spl_data(i_data,1,i_r_z,i_r_phi,:,i_z,i_phi)
+                    call spl_reg(spl_order,n_r,h_r,splcoe)
+                    do k=1,spl_order
+                        spl_data(i_data,k+1,i_r_z,i_r_phi,:,i_z,i_phi)=splcoe(k,:)
+                    enddo
+                enddo
+            enddo
+            enddo
+        enddo
+        enddo
+        deallocate(splcoe)
+
+    end subroutine construct_splines_3d
+
+end module interpolate

tests/CMakeLists.txt

@@ -76,6 +76,9 @@ target_link_libraries(test_system_utility.x
 add_executable(test_util.x source/test_util.f90)
 target_link_libraries(test_util.x ${COMMON_LIBS})
 
+add_executable(test_interpolate.x source/test_interpolate.f90)
+target_link_libraries(test_interpolate.x ${COMMON_LIBS})
+
 
 ## Standard tests according to standard libneo Fortran test convention.
 
@@ -90,6 +93,7 @@ add_fortran_test(test_efit_class)
 add_fortran_test(test_geqdsk_tools)
 add_fortran_test(test_hdf5_tools)
 add_fortran_test(test_util)
+add_fortran_test(test_interpolate)
 
 ## Custom tests
 

tests/source/test_interpolate.f90

@@ -1,32 +1,65 @@
 program test_interpolate
-    use interpolate
+    use math_constants
+    use util
 
     implicit none
 
     integer, parameter :: N_POINTS = 100
-    real(8), parameter :: X_MIN = 0.0d0, X_MAX = 2.0d0 * pi
+    real(8), parameter :: X_MIN = 0.0d0, X_MAX = TWOPI
 
-    real(8), dimension(N_POINTS) :: x, y, dy, d2y
-
-    call linspace(x, 0.0d0, 2.0d0 * pi, 100)
-
-    call test_interpolate_1d
+    call test_spl_reg(spline_order=3)
 
 contains
 
     subroutine generate_test_data_1d(x, y, dy, d2y)
-        real, dimension(:), intent(in)  :: x
-        real, dimension(:), intent(out) :: y, dy, d2y
-        integer :: i
+        real(8), dimension(:), intent(in)  :: x
+        real(8), dimension(:), intent(out) :: y, dy, d2y
 
         y = cos(x)
         dy = -sin(x)
         d2y = -cos(x)
 
     end subroutine generate_test_data_1d
 
-    subroutine test_interpolate_1d
+    subroutine test_spl_reg(spline_order)
+
+        integer, intent(in) :: spline_order
+
+        real(8), dimension(N_POINTS) :: x, y, dy, d2y
+
+        real(8) :: expected, actual
+
+        real(kind=real_kind), dimension(0:spline_order, N_POINTS) :: spline_coeff
+
+        real(8) :: h_step, x_norm
+        integer :: i, k, point_index
+
+        call linspace(0.0d0, 2.0d0 * pi, 100, x)
+        call generate_test_data_1d(x, y, dy, d2y)
+
+        h_step = x(2) - x(1)
+
+        spline_coeff(0,:) = y
+        call spl_reg(spline_order, N_POINTS, h_step, spline_coeff)
+
+        i = 30
+
+        expected = cos(x(i))
+
+        x_norm = x(i)/h_step
+        point_index = max(0, min(N_POINTS-1, int(x_norm)))
+        x_norm = (x_norm - dble(point_index))*h_step
+        point_index = point_index + 1
+
+        actual = spline_coeff(spline_order+1, point_index)
+
+        do k = spline_order, 0, -1
+            actual = spline_coeff(k, point_index) + x_norm*actual
+        enddo
+
+        print *, 'expected = ', expected
+        print *, 'actual = ', actual
 
-    end subroutine test_interpolate_1d
+    end subroutine test_spl_reg
 
 end program test_interpolate