Finished cubic spline and tweaked errors.

Author: nshaffer

error.f90

@@ -1,22 +1,22 @@
 module error
   implicit none
 contains
-  subroutine throw(msg, u)
-    character(len=*), intent(in) :: msg
-    integer, optional, intent(in) :: u
-    if (present(u)) then 
-       write(u, *) "Error: " // trim(msg)
+  subroutine throw(msg, unit)
+    character(*), intent(in) :: msg
+    integer, optional, intent(in) :: unit
+    if (present(unit)) then
+       write(unit, *) "Error: " // trim(msg)
     else
        print *, "Error: " // trim(msg)
     end if
     stop
   end subroutine throw
 
-  subroutine warn(msg, u)
-    character(len=*), intent(in) :: msg
-    integer, optional, intent(in) :: u
-    if (present(u)) then 
-       write(u, *) "Warning: " // trim(msg)
+  subroutine warn(msg, unit)
+    character(*), intent(in) :: msg
+    integer, optional, intent(in) :: unit
+    if (present(unit)) then
+       write(unit, *) "Warning: " // trim(msg)
     else
        print *, "Warning: " // trim(msg)
     end if

error.mod

@@ -20,26 +20,6 @@ subroutine linear_coeff(xx, yy, bb, n)
     end do
   end subroutine linear_coeff
 
-  real(dp) function linear_eval(x, xx, yy, bb, n) result(s)
-    ! Evaluates the piecewise linear interpolant
-    !   S_i(x) = y_i + b_i*(x-x_i)
-    ! at a point x.
-    integer, intent(in) :: n
-    real(dp), intent(in) :: x
-    real(dp), dimension(n), intent(in) :: xx, yy
-    real(dp), dimension(n-1), intent(in) :: bb
-
-    integer :: i
-
-    if ( (x.lt.xx(1)).or.(x.gt.xx(n)) ) then
-       call warn("Requested abcissa lies outside the given range. Defaulting to 0.d0.")
-       s = 0.d0
-    else
-       i = isearch(x, xx, n)
-       s = yy(i) + bb(i)*(x - xx(i))
-    end if
-  end function linear_eval
-    
   subroutine cubic_coeff(xx, yy, bb, cc, dd, n)
     ! Computes the N coefficients of the natural cubic spline fitting
     ! the set of N knots (xx, y). 
@@ -110,6 +90,26 @@ subroutine cubic_coeff(xx, yy, bb, cc, dd, n)
     cc(n) = 3.d0*cc(n)
     dd(n) = dd(nm1)
   end subroutine cubic_coeff
+
+  real(dp) function linear_eval(x, xx, yy, bb, n) result(s)
+    ! Evaluates the piecewise linear interpolant
+    !   S_i(x) = y_i + b_i*(x-x_i)
+    ! at a point x.
+    integer, intent(in) :: n
+    real(dp), intent(in) :: x
+    real(dp), dimension(n), intent(in) :: xx, yy
+    real(dp), dimension(n-1), intent(in) :: bb
+
+    integer :: i
+
+    if ( (x.lt.xx(1)).or.(x.gt.xx(n)) ) then
+       call warn("Requested abcissa lies outside the given range. Defaulting to 0.d0.")
+       s = 0.d0
+    else
+       i = isearch(x, xx, n)
+       s = yy(i) + bb(i)*(x - xx(i))
+    end if
+  end function linear_eval
 
   real(dp) function cubic_eval(x, xx, yy, bb, cc, dd, n) result(s)
     ! Evaluates the cubic spline
@@ -121,8 +121,8 @@ real(dp) function cubic_eval(x, xx, yy, bb, cc, dd, n) result(s)
     real(dp), dimension(n), intent(in) :: bb, cc, dd
     integer :: i
     if ( (x.lt.xx(1)).or.(x.gt.xx(n)) ) then
+       call warn("Requested abcissa lies outside the given range. Defaulting to 0.")
        s = 0.d0
-       call warn("Requested abcissa lies outside the given range. Defaulting to 0.d0.")
     else
        i = isearch(x, xx, n)
        s = yy(i) + bb(i)*(x-xx(i)) + cc(i)*(x-xx(i))**2 + dd(i)*(x-xx(i))**3

error.o

@@ -0,0 +1,30 @@
+program test_cubic
+  use types, only: dp
+  use interp, only: cubic_coeff, cubic_eval
+  implicit none
+  
+  integer, parameter :: n = 5
+  real(dp), dimension(n) :: xx, yy
+  real(dp), dimension(n) :: bb, cc, dd
+  integer :: i
+  real(dp) :: x, y
+  real(dp), dimension(500) :: rr, ss
+
+  xx(1) = 0.d0
+  do i = 2, n
+     xx(i) = xx(i-1) + 0.5d0
+  end do
+  yy = sin(xx)
+
+  call cubic_coeff(xx, yy, bb, cc, dd, n)
+  
+  x = 1.2d0
+  y = cubic_eval(x, xx, yy, bb, cc, dd, n) ! Should be 0.93148217019426938
+  print *, y
+
+  x = 4.d0
+  y = cubic_eval(x, xx, yy, bb, cc, dd, n) ! Should warn and be 0.d0
+  print *, y
+
+end program test_cubic
+