! $###########################################
! Translate from Mr.Michael's code
! http://mysite.verizon.net/~vze2vrva/  
! ############################################
 
module BicubicSplineMod
  implicit none
  integer in_width, in_height, out_width, out_height
  real,dimension(:,:),allocatable :: f,g
  real a
contains

subroutine splineInit()
  integer u,v, n, d, i,j, k, m, index, il
  real,allocatable,dimension(:) :: L,h
  real,allocatable,dimension(:,:) :: c
  real x

  call random_seed () 
  in_width=4; in_height=in_width
  out_width=40; out_height=out_width
  a=-0.5; d=in_width; n=out_width

  allocate(f(0:in_width-1,0:in_height-1))
  allocate(g(0:out_width-1,0:out_width-1))
  allocate(L(0:out_width-1),c(0:3,0:out_width-1),h(0:out_width-1))

  do u=0,in_width-1
     do v=0,in_width-1
        ipoint=ipoint + 1
        points(1,ipoint) = u; points(2,ipoint) = v
        call random_number(x)
        f(u,v)=points(3,ipoint)        
     end do
  end do
  do k=0,out_width-1
     L(k) = k * d/ n
  end do
  do k=0, n-1
     i=mod(k*d,n)
     x = real(i) / real(n)
     c(0,k) = C0(x); c(1,k) = C1(x)
     c(2,k) = C2(x); c(3,k) = C3(x)
  end do
  do k=n,out_width-1
     do il=0,3
        i=mod(k,n)
        c(il,k) = c(il,i)
     end do
  end do
  do k=0,out_height-1
     do j=0,in_width-1
        h(j) = 0.0;
        do il=0,3
           index = L(k) + il - 1;
           if ((index >= 0) .and. (index < in_height) ) then
              h(j) = h(j) + f(index,j) * c(3 - il,k)
           end if
        end do
     end do
     do m=0,out_width-1
        x = 0.5
        do il=0,3
           index = L(m) + il - 1;
           if ((index >= 0) .and. (index < in_width) ) then
              x = x + h(index) * c(3 - il,m)
           end if
        end do       
        g(k,m) = x
     end do
  end do

! output g vector
! x,y,z is just the coordinate of interpoints points
  do u=0,out_width-1
     x=real(u)*real(in_width-1)/real(out_width-1)
     do v=0,out_width-1
        y=real(v)*real(in_width-1)/real(out_width-1)
        z=g(u,v)
     end do
  end do

end subroutine

function C0(t)
  real t,C0
  C0 = -a * t * t * t + a * t * t
end function

function C1(t)
  real t,C1
  C1 = -(a + 2.0) * t * t * t + (2.0 * a + 3.0) * t * t - a * t
end function

function C2(t)
  real t,C2
  C2 = (a + 2.0) * t * t * t - (a + 3.0) * t * t + 1.0
end function

function C3(t)
  real t,C3
  C3 = a * t * t * t - 2.0 * a * t * t + a * t
end function

end module BicubicSplineMod