!!.. This module generates the cubic B-spline basis functions
module base_spline_mod
  use OW_type
  implicit none

  type(axis), private :: axd, ayd, azd
  type(axis), private :: axk, ayk, azk

contains

  integer function base_spln1d_init(axd_in,axk_in) result(stat)
    type(axis) :: axd_in
    type(axis) :: axk_in
    axd  = axd_in
    axk  = axk_in
    stat = 0
  end function

  integer function base_spln2d_init(axd_in,ayd_in,axk_in,ayk_in) result(stat)
    type(axis) :: axd_in,ayd_in
    type(axis) :: axk_in,ayk_in
    axd = axd_in; ayd = ayd_in
    axk = axk_in; ayk = ayk_in
    stat = 0
  end function

  integer function base_spln3d_init(axd_in,ayd_in,azd_in,axk_in,ayk_in,azk_in) result(stat)
    type(axis) :: axd_in,ayd_in,azd_in
    type(axis) :: axk_in,ayk_in,azk_in
    axd = axd_in; ayd = ayd_in; azd = azd_in
    axk = axk_in; ayk = ayk_in; azk = azk_in
    stat = 0
  end function

  integer function base_spln1d_clos() result(stat)
    stat = 0
  end function

  integer function base_spln2d_clos() result(stat)
    stat = 0
  end function

  integer function base_spln3d_clos() result(stat)
    stat = 0
  end function

  integer function base_spln1d_oper(adj, add, modl, data) result(stat)
    logical                :: adj, add
    real, dimension(axk%n) :: modl
    real, dimension(axd%n) :: data

    integer                :: ixd, ixk
    real                   :: xd, xk
    if (.not.add) then
       if (adj) then
          modl = 0.
       else 
          data = 0.
       end if
    end if

    if (adj) then
       do ixd=1, axd%n; xd = axd%o+(ixd-1)*axd%d
       do ixk=1, axk%n; xk = axk%o+(ixk-1)*axk%d
          modl(ixk) = modl(ixk) + data(ixd)*spline(xd-xk,axk%d)
       end do 
       end do
    else
       do ixd=1, axd%n; xd = axd%o+(ixd-1)*axd%d
       do ixk=1, axk%n; xk = axk%o+(ixk-1)*axk%d
          data(ixd) = data(ixd) + modl(ixk)*spline(xd-xk,axk%d)
       end do
       end do
    end if
    stat = 0
  end function


  integer function base_spln2d_oper(adj, add, modl, data) result(stat)
    logical            :: adj, add
    real, dimension(:) :: modl
    double precision, dimension(:) :: data

    if (.not.add) then
       if (adj) then
          modl = 0.
       else
          data = 0.
       end if
    end if
    call base_spln2d(adj, modl, data)
    stat = 0 
  end function

  subroutine base_spln2d(adj, modl, data)
    logical                      :: adj
    real, dimension(axk%n,ayk%n) :: modl
    double precision, dimension(axd%n,ayd%n) :: data

    integer                      :: ixd, ixk
    integer                      :: iyd, iyk
    real                         :: xd, yd
    real                         :: xk, yk

    if (adj) then
       do iyd=1, ayd%n; yd = ayd%o+(iyd-1)*ayd%d 
       do ixd=1, axd%n; xd = axd%o+(ixd-1)*axd%d
          do iyk=1, ayk%n; yk = ayk%o+(iyk-1)*ayk%d
          do ixk=1, axk%n; xk = axk%o+(ixk-1)*axk%d
             modl(ixk,iyk) = modl(ixk,iyk) + data(ixd,iyd)*spline(xd-xk,axk%d)*spline(yd-yk,ayk%d)
          end do
          end do
       end do
       end do
    else
       do iyd=1, ayd%n; yd = ayd%o+(iyd-1)*ayd%d
       do ixd=1, axd%n; xd = axd%o+(ixd-1)*axd%d
          do iyk=1, ayk%n; yk = ayk%o+(iyk-1)*ayk%d
          do ixk=1, axk%n; xk = axk%o+(ixk-1)*axk%d
             data(ixd,iyd) = data(ixd,iyd) + modl(ixk,iyk)*spline(xd-xk,axk%d)*spline(yd-yk,ayk%d)
          end do
          end do
       end do
       end do
    end if
  end subroutine

  integer function base_spln3d_oper(adj, add, modl, data) result(stat)
    logical            :: adj, add
    real, dimension(:) :: modl
!    real, dimension(:) :: data
    double precision, dimension(:) :: data

    if (.not.add) then
       if (adj) then
          modl = 0.
       else
          data = 0.
       end if
    end if

    call base_spln3d(adj, modl, data)

    stat = 0
  end function
  
  subroutine base_spln3d(adj, modl, data)
    logical                            :: adj
    real, dimension(axk%n,ayk%n,azk%n) :: modl
!    real, dimension(axd%n,ayd%n,azd%n) :: data
    double precision, dimension(axd%n,ayd%n,azd%n) :: data

    integer                            :: ixd, ixk
    integer                            :: iyd, iyk
    integer                            :: izd, izk
    real                               :: xd, yd, zd
    real                               :: xk, yk, zk

    if (adj) then
       do izd=1, azd%n; zd = azd%o+(izd-1)*azd%d
       do iyd=1, ayd%n; yd = ayd%o+(iyd-1)*ayd%d
       do ixd=1, axd%n; xd = axd%o+(ixd-1)*axd%d
          do izk=1, azk%n; zk = azk%o+(izk-1)*azk%d
          do iyk=1, ayk%n; yk = ayk%o+(iyk-1)*ayk%d
          do ixk=1, axk%n; xk = axk%o+(ixk-1)*axk%d
             modl(ixk,iyk,izk) = modl(ixk,iyk,izk) + data(ixd,iyd,izd)                       &
                               * spline(xd-xk,axk%d)*spline(yd-yk,ayk%d)*spline(zd-zk,azk%d)
          end do
          end do
          end do
       end do
       end do
       end do
    else
       do izd=1, azd%n; zd = azd%o+(izd-1)*azd%d
       do iyd=1, ayd%n; yd = ayd%o+(iyd-1)*ayd%d
       do ixd=1, axd%n; xd = axd%o+(ixd-1)*axd%d
          do izk=1, azk%n; zk = azk%o+(izk-1)*azk%d
          do iyk=1, ayk%n; yk = ayk%o+(iyk-1)*ayk%d
          do ixk=1, axk%n; xk = axk%o+(ixk-1)*axk%d
             data(ixd,iyd,izd) = data(ixd,iyd,izd) + modl(ixk,iyk,izk)                       &
                               * spline(xd-xk,axk%d)*spline(yd-yk,ayk%d)*spline(zd-zk,azk%d)
          end do
          end do
          end do
       end do
       end do
       end do
    end if
  end subroutine

  real function spline(x, dx)
    real :: x, dx
    if (x<-2*dx) then
       spline = 0.
    else if (x>=-2*dx .and. x<-dx) then
       spline = (2*dx+x)**3/6./dx**3
    else if (x>=-dx .and. x<0.) then
       spline = (-x)**3/2./dx**3-(-x)**2/dx**2+2./3.
    else if (x>=0 .and. x<dx) then
       spline = (+x)**3/2./dx**3-(+x)**2/dx**2+2./3.
    else if (x>=dx .and. x<=2*dx) then
       spline = (2*dx-x)**3/6./dx**3
    else
       spline = 0.
    end if
  end function

end module
