Cardinal B-Spline

c **********************************************************************
c     spline.F
c **********************************************************************
      function spline(u,order)
c **********************************************************************
      implicit none
      real u, spline, M2,M3,M4,M5,M6
      integer order
      select case (order)
      case (6)
         spline = M6(u)
      case (5)
         spline = M5(u)
      case (4)
         spline = M4(u)
      case (3)
         spline = M3(u)
      case (2)
         spline = M2(u)
      case default
         print *, "spline not defined for order ",order
      end select

      return
      end

c **********************************************************************
      function M2(u)
c **********************************************************************
      implicit none
      integer, parameter :: one = 1.0d0
      real u, M2, M3
      M2 = 0
      if (u>0) then
         if (u<2) then
            M2 = one - abs(u-one)
         end if
      end if
      return
      end            
c **********************************************************************
      function M3(u)
      implicit none
      integer, parameter :: one = 1.0d0 
      real u, M2, M3
     
      M3 = u/2.0d0 * M2(u) + (3.0d0-u)/2.0d0 * M2(u-1.0d0)
      return
      end
c **********************************************************************
      function M4(u)
      implicit none
      integer, parameter :: one = 1.0d0 
      real u, M3, M4

      M4 = u/3.0d0 * M3(u) + (4.0d0-u)/3.0d0 * M3(u-1.0d0)
      return
      end
c **********************************************************************
      function M5(u)
      implicit none
      integer, parameter :: one = 1.0d0
      real u, M4, M5

      M5 = u/4.0d0 * M4(u) + (5.0d0-u)/4.0d0 * M4(u-one)
      return
      end
c **********************************************************************
      function M6(u)
      implicit none
      integer, parameter :: one = 1.0d0 
      real u, M5, M6

      M6 = u/5.0d0 * M5(u) + (6.0d0-u)/5.0d0 * M5(u-one)
      return
      end