module derivatives

  use types, only: dp
  use data, only: Nmesh, dK1, Dk2, Dk3, dx, dn1, dn2, dn3, Dcoeffs1, Dcoeffs2, Dcoeffs3
  use data, only: r, Dr, DDr, DDDr
  implicit none

  private
  public genDcoeffs,xDiff1,xDiff2,xDiff3,rDiff1,rDiff2,rDiff3,Weizsacker

contains

  function facArray(array) result(fac)
    integer, intent(in) :: array(:)
    integer :: fac
    integer :: i

    fac=1
    do i=1,size(array)
       fac=fac*array(i)
    end do
  end function facArray

  function facInteger(number) result(fac)
    integer, intent(in) :: number
    integer :: fac
    integer :: i

    if (number<=0) then
       fac=1
    else
       fac=1
       do i=1,number
          fac=fac*i
       end do
    end if
  end function facInteger

  function genDcoeffs(order,stencil,h) result(d)
    integer, intent(in) :: stencil, order
    real(dp), intent(in) :: h
    real(dp), allocatable :: d(:,:)
    integer, allocatable :: a(:,:), array(:), facThese(:)
    integer :: n,m,i,j,k,l,o,p,index,index2

    n=stencil-1
    m=order
    allocate(a(n+1,n+1))
    allocate(d(n+1,n+1))
    a=0
    d=0.0_dp

    ! First, calculate a(m-1)_{n-1,i,j} = a_{m-1}(-i,...,k-i,...,n-i)
    ! k = 0,...,n and k != i or j

    do j=0,n
       do i=0,n
          if (i==j) then
             index=1
             allocate(array(n))
          else
             index=1
             allocate(array(n-1))
          end if
          do k=0,n
             if (k==i .or. k==j) then
                cycle
             else
                array(index)=-i+k
                index=index+1
             end if
          end do
          if (m-1==0) then
             a(i+1,j+1)=facArray(array)
          else if (m-1==1) then
             allocate(facThese(size(array)-1))
             do l=1,size(array)
                facThese=0
                index2=1
                do o=1,size(array)
                   if (o==l) then
                      cycle
                   else
                      facThese(index2)=array(o)
                      index2=index2+1
                   end if
                end do
                a(i+1,j+1)=a(i+1,j+1)+facArray(facThese)
             end do
             deallocate(facThese)
          else if (m-1==2) then
             allocate(facThese(size(array)-2))
             do l=1,size(array)
                do p=l+1,size(array)
                   facThese=0
                   index2=1
                   do o=1,size(array)
                      if (o==l .or. o==p) then
                         cycle
                      else
                         facThese(index2)=array(o)
                         index2=index2+1
                      end if
                   end do
                   a(i+1,j+1)=a(i+1,j+1)+facArray(facThese)
                end do
             end do
             deallocate(facThese)
          end if
          deallocate(array)
       end do
    end do

    ! Next, compute d(m)_{n+1,i,j}'s

    do j=0,n
       do i=0,n
          if (i==j) then
             cycle
          else
             d(i+1,j+1)=real((-1)**(m-j)*facInteger(m)*a(i+1,j+1),dp)/facInteger(j)/facInteger(n-j)/h**m
          end if
       end do
    end do

    do i=0,n
       do j=0,n
          if (j==i) then
             cycle
          else
             d(i+1,i+1)=d(i+1,i+1)-d(i+1,j+1)
          end if
       end do
    end do
  end function genDcoeffs

  function xDiff1(array) result(Darray)
    real(dp), intent(in) :: array(Nmesh)
    real(dp) :: Darray(Nmesh)
    integer :: i,j
    Darray=0.0_dp
    do i=1,dK1
       do j=1,dn1+1
          Darray(i)=sum(Dcoeffs1(i,1:dn1+1)*array(1:dn1+1))
       end do
    end do
    do i=dK1+1,Nmesh-dK1
       do j=1,dn1+1
          Darray(i)=sum(Dcoeffs1(dK1+1,1:dn1+1)*array(i-dK1:i-dK1+dn1))
       end do
    end do
    do i=Nmesh-dK1+1,Nmesh
       do j=1,dn1+1
          Darray(i)=sum(Dcoeffs1(dn1+i-Nmesh+1,1:dn1+1)*array(Nmesh-dn1:Nmesh))
       end do
    end do
  end function xDiff1

  function xDiff2(array) result(Darray)
    real(dp), intent(in) :: array(Nmesh)
    real(dp) :: Darray(Nmesh)
    integer :: i,j
    Darray=0.0_dp
    do i=1,dK2
       do j=1,dn2+1
          Darray(i)=sum(Dcoeffs2(i,1:dn2+1)*array(1:dn2+1))
       end do
    end do
    do i=dK2+1,Nmesh-dK2
       do j=1,dn2+1
          Darray(i)=sum(Dcoeffs2(dK2+1,1:dn2+1)*array(i-dK2:i-dK2+dn2))
       end do
    end do
    do i=Nmesh-dK2+1,Nmesh
       do j=1,dn2+1
          Darray(i)=sum(Dcoeffs2(dn2+i-Nmesh+1,1:dn2+1)*array(Nmesh-dn2:Nmesh))
       end do
    end do
  end function xDiff2

  function xDiff3(array) result(Darray)
    real(dp), intent(in) :: array(Nmesh)
    real(dp) :: Darray(Nmesh)
    integer :: i,j
    Darray=0.0_dp
    do i=1,dK3
       do j=1,dn3+1
          Darray(i)=sum(Dcoeffs3(i,1:dn3+1)*array(1:dn3+1))
       end do
    end do
    do i=dK3+1,Nmesh-dK3
       do j=1,dn3+1
          Darray(i)=sum(Dcoeffs3(dK3+1,1:dn3+1)*array(i-dK3:i-dK3+dn3))
       end do
    end do
    do i=Nmesh-dK3+1,Nmesh
       do j=1,dn3+1
          Darray(i)=sum(Dcoeffs3(dn3+i-Nmesh+1,1:dn3+1)*array(Nmesh-dn3:Nmesh))
       end do
    end do
  end function xDiff3

  function rDiff1(array) result(Darray)
    real(dp), intent(in) :: array(Nmesh)
    real(dp) :: Darray(Nmesh)
    Darray=xDiff1(array)/Dr
  end function rDiff1

  function rDiff2(array) result(Darray)
    real(dp), intent(in) :: array(Nmesh)
    real(dp) :: Darray(Nmesh)
    Darray=xDiff2(array)/Dr**2-xDiff1(array)*DDr/Dr**3
  end function rDiff2

  function rDiff3(array) result(Darray)
    real(dp), intent(in) :: array(Nmesh)
    real(dp) :: Darray(Nmesh)
    Darray=xDiff3(array)/Dr**3-3*xDiff2(array)*DDr/Dr**4+xDiff1(array)*(3*DDr**2/Dr**5-DDDr/Dr**4)
  end function rDiff3

function Weizsacker(Rho) result(tW)
  real(dp), intent(in) :: Rho(Nmesh)
  real(dp) :: tW(Nmesh)
  tW=-(0.5_dp*rDiff2(Rho)/Rho-0.25_dp*(rDiff1(Rho)/Rho)**2)/2.0_dp
end function Weizsacker

end module derivatives
