! check_test_4d.f90  check every derivative function
!                    from package test_4d

module global
  double precision :: xmin, xmax, ymin, ymax, zmin, zmax, tmin, tmax
  double precision :: hx, hy, hz, ht
  integer :: nx, ny, nz, nt
  double precision, dimension(10) :: xg, yg, zg, tg
  double precision :: gmaxerror
  interface
    subroutine rderiv1(order, npoints, point, h, c)
      integer, intent(in) :: order    ! derivative order
      integer, intent(in) :: npoints  ! number of points to be used
      integer, intent(in) :: point    ! compute derivative at this point
      double precision, intent(in) :: h  ! divisor for use h^order
      double precision, intent(out), dimension(*)::c ! returned coefficients
    end subroutine rderiv1
  end interface
end module global

subroutine check(name, F, points, xord, yord, zord, tord)
  use global
  use test_4d
  implicit none
  character(len=*), intent(in) :: name
  integer, intent(in) :: points, xord, yord, zord, tord
  interface
    function F(x, y, z, t) result(val)
      double precision :: x, y, z, t, val
    end function F
  end interface
  double precision, dimension(10) ::  cx, cy, cz, ct
  double precision, dimension(10) ::  p1, p2, p3
  double precision :: analytic, discrete, error, maxerror
  double precision ::  x, y, z, t
  integer :: i, ii, iii, iiii
  integer :: j, jj, jjj, jjjj

  maxerror = 0.0
  do i = 1,points 
    call rderiv1(xord, points, i, hx, cx)
    x = xg(i)
    do ii = 1,points 
      call rderiv1(yord, points, ii, hy, cy)
      y = yg(ii)
      do iii = 1,points 
        call rderiv1(zord, points, iii, hz, cz)
        z = zg(iii)
        do iiii = 1,points 
        call rderiv1(tord, points, iiii, ht, ct)
        t = tg(iiii)
        analytic = F(x, y, z, t)
        discrete = 0.0

        ! four cases of one partial x, y, z, t to any order
        if(xord.ne.0.and.yord.eq.0.and.zord.eq.0.and.tord.eq.0) then
          do j = 1,points 
            discrete = discrete + cx(j)*u(xg(j),y,z,t)
          end do ! j

        else if(xord.eq.0.and.yord.ne.0.and.zord.eq.0.and.tord.eq.0) then
          do jj = 1,points 
            discrete = discrete + cy(jj)*u(x,yg(jj),z,t)
          end do ! jj

        else if(xord.eq.0.and.yord.eq.0.and.zord.ne.0.and.tord.eq.0) then
          do jjj = 1,points 
            discrete = discrete + cz(jjj)*u(x,y,zg(jjj),t)
          end do ! jjj

        else if(xord.eq.0.and.yord.eq.0.and.zord.eq.0.and.tord.ne.0) then
          do jjjj = 1,points 
            discrete = discrete + ct(jjjj)*u(x,y,z,tg(jjjj))
          end do ! jjjj

        ! six cases of two partials xy, xz, xt, yz, yt, zt to any order
        else if(xord.ne.0.and.yord.ne.0.and.zord.eq.0.and.tord.eq.0) then
          do j = 1,points 
            p1(j) = 0.0
            do jj = 1,points 
              p1(j) = p1(j) + cy(jj)*u(xg(j),yg(jj),z,t)
            end do ! jj
            discrete = discrete + cx(j)*p1(j)
          end do ! j

        else if(xord.ne.0.and.yord.eq.0.and.zord.ne.0.and.tord.eq.0) then
          do j = 1,points 
            p1(j) = 0.0
            do jjj = 1,points 
              p1(j) = p1(j) + cz(jjj)*u(xg(j),y,zg(jjj),t)
            end do ! jjj
            discrete = discrete + cx(j)*p1(j)
          end do ! j

        else if(xord.ne.0.and.yord.eq.0.and.zord.eq.0.and.tord.ne.0) then
        do j = 1,points 
          p1(j) = 0.0
          do jjjj = 1,points 
            p1(j) = p1(j) + ct(jjjj)*u(xg(j),y,z,tg(jjjj))
          end do ! jjjj
          discrete = discrete + cx(j)*p1(j)
        end do ! j

        else if(xord.eq.0.and.yord.ne.0.and.zord.ne.0.and.tord.eq.0) then
          do jj = 1,points 
            p1(jj) = 0.0
            do jjj = 1,points 
              p1(jj) = p1(jj) + cz(jjj)*u(x,yg(jj),zg(jjj),t)
            end do ! jjj
            discrete = discrete + cy(jj)*p1(jj)
          end do ! jj

        else if(xord.eq.0.and.yord.ne.0.and.zord.eq.0.and.tord.ne.0) then
          do jj = 1,points 
            p1(jj) = 0.0
            do jjjj = 1,points 
              p1(jj) = p1(jj) + ct(jjjj)*u(x,yg(jj),z,tg(jjjj))
            end do ! jjjj
            discrete = discrete + cy(jj)*p1(jj)
          end do ! jj

        else if(xord.eq.0.and.yord.eq.0.and.zord.ne.0.and.tord.ne.0) then
          do jjj = 1,points 
            p1(jjj) = 0.0
            do jjjj = 1,points 
              p1(jjj) = p1(jjj) + ct(jjjj)*u(x,y,zg(jjj),tg(jjjj))
            end do ! jjjj
            discrete = discrete + cz(jjj)*p1(jjj)
          end do ! jjj

        ! four cases of three partials xyz, xyt, xzt, yzt to any order
        else if(xord.ne.0.and.yord.ne.0.and.zord.ne.0.and.tord.eq.0) then
          do j = 1,points 
            p1(j) = 0.0
            do jj = 1,points 
              p2(jj) = 0.0
              do jjj = 1,points 
                p2(jj) = p2(jj) + cz(jjj)*u(xg(j),yg(jj),zg(jjj),t)
              end do ! jjj
              p1(j) = p1(j) + cy(jj)*p2(jj)
            end do ! jj
            discrete = discrete + cx(j)*p1(j)
          end do ! j

        else if(xord.ne.0.and.yord.ne.0.and.zord.eq.0.and.tord.ne.0) then
          do j = 1,points 
            p1(j) = 0.0
            do jj = 1,points 
              p2(jj) = 0.0
              do jjjj = 1,points 
                p2(jj) = p2(jj) + ct(jjjj)*u(xg(j),yg(jj),z,tg(jjjj))
              end do ! jjjj
              p1(j) = p1(j) + cy(jj)*p2(jj)
            end do ! jj
            discrete = discrete + cx(j)*p1(j)
          end do ! j

        else if(xord.ne.0.and.yord.eq.0.and.zord.ne.0.and.tord.ne.0) then
          do j = 1,points 
            p1(j) = 0.0
            do jjj = 1,points 
              p2(jjj) = 0.0
              do jjjj = 1,points 
                p2(jjj) = p2(jjj) + ct(jjjj)*u(xg(j),y,zg(jjj),tg(jjjj))
              end do ! jjjj
              p1(j) = p1(j) + cz(jjj)*p2(jjj)
            end do ! jjj
            discrete = Discrete + cx(j)*p1(j)
          end do ! j

        else if(xord.eq.0.and.yord.ne.0.and.zord.ne.0.and.tord.ne.0) then
          do jj = 1,points 
            p1(jj) = 0.0
            do jjj = 1,points 
              p2(jjj) = 0.0
              do jjjj = 1,points 
                p2(jjj) = p2(jjj) + ct(jjjj)*u(x,yg(jj),zg(jjj),tg(jjjj))
              end do ! jjjj
              p1(jj) = p1(jj) + cz(jjj)*p2(jjj)
            end do ! jjj
            discrete = discrete + cy(jj)*p1(jj)
          end do ! jj

        ! one case of four partials
        else
          do j = 1,points 
            p1(j) = 0.0
            do jj = 1,points 
              p2(jj) = 0.0
              do jjj = 1,points 
                p3(jjj) = 0.0
                do jjjj = 1,points 
                  p3(jjj) = p3(jjj) +ct(jjjj)*u(xg(j),yg(jj),zg(jjj),tg(jjjj))
                end do ! jjjj
                p2(jj) = p2(jj) + cz(jjj)*p3(jjj)
                end do ! jjj
              p1(j) = p1(j) + cy(jj)*p2(jj)
              end do ! jj
            discrete = Discrete + cx(j)*p1(j)
            end do ! j
          end if
          error = analytic-discrete
          if(abs(Error).gt.maxerror) then
            Maxerror = abs(error)
          end if
          if(maxerror.gt.0.001.and.i.eq.0.and.ii.eq.0.and.iii.eq.0.and.iiii.eq.0) then
            print *,'analytic=',analytic,', discrete=',discrete, &
                    ', error=',error
          end if
        end do ! iiii  t
      end do ! iii  z
    end do ! ii  y
  end do ! i  x
  print *, '                                 ',name,' maxerror=',maxerror
  if(maxerror.gt.gmaxerror) then
    gmaxerror = maxerror
  end if
end subroutine check


program check_test_4d
  use test_4d
  use global
  implicit none
  integer :: i, ii, iii, iiii

  print *, 'check_test_4d.f90 running'
  ! a test case, all do not need to be same
  xmin = -1.0
  xmax = 1.0
  nx = 5
  hx = (xmax-xmin)/(nx-1)
  do i = 1,nx
    xg(i) = xmin + (i-1)*hx
  end do
  ymin = -1.0
  ymax = 1.0
  ny = 5
  hy = (ymax-ymin)/(ny-1)
  do ii = 1,ny 
    yg(ii) = ymin + (ii-1)*hy
  end do
  zmin = -1.0
  zmax = 1.0
  nz = 5
  hz = (zmax-zmin)/(nz-1)
  do iii = 1,nz 
    zg(iii) = zmin + (iii-1)*hz
  end do
  tmin = -1.0
  tmax = 1.0
  nt = 5
  ht = (tmax-tmin)/(nt-1)
  do iiii = 1,nt 
    tg(iiii) = tmin + (iiii-1)*ht
  end do
  gmaxerror = 0.0
  call check("Ux   ", Ux,    5, 1, 0, 0, 0)
  call check("Uy   ", Uy,    5, 0, 1, 0, 0)
  call check("Uz   ", Uz,    5, 0, 0, 1, 0)
  call check("Ut   ", Ut,    5, 0, 0, 0, 1)
  call check("Uxx  ", Uxx,   5, 2, 0, 0, 0)
  call check("Uxy  ", Uxy,   5, 1, 1, 0, 0)
  call check("Uxz  ", Uxz,   5, 1, 0, 1, 0)
  call check("Uxt  ", Uxt,   5, 1, 0, 0, 1)
  call check("Uyy  ", Uyy,   5, 0, 2, 0, 0)
  call check("Uyz  ", Uyz,   5, 0, 1, 1, 0)
  call check("Uyt  ", Uyt,   5, 0, 1, 0, 1)
  call check("Uzz  ", Uzz,   5, 0, 0, 2, 0)
  call check("Uzt  ", Uzt,   5, 0, 0, 1, 1)
  call check("Utt  ", Utt,   5, 0, 0, 0, 2)
  call check("Uxxx ", Uxxx,  5, 3, 0, 0, 0)
  call check("Uxxy ", Uxxy,  5, 2, 1, 0, 0)
  call check("Uxxz ", Uxxz,  5, 2, 0, 1, 0)
  call check("Uxxt ", Uxxt,  5, 2, 0, 0, 1)
  call check("Uxyy ", Uxyy,  5, 1, 2, 0, 0)
  call check("Uxyz ", Uxyz,  5, 1, 1, 1, 0)
  call check("Uxyt ", Uxyt,  5, 1, 1, 0, 1)
  call check("Uxzz ", Uxzz,  5, 1, 0, 2, 0)
  call check("Uxzt ", Uxzt,  5, 1, 0, 1, 1)
  call check("Uxtt ", Uxtt,  5, 1, 0, 0, 2)
  call check("Uyyy ", Uyyy,  5, 0, 3, 0, 0)
  call check("Uyyz ", Uyyz,  5, 0, 2, 1, 0)
  call check("Uyyt ", Uyyt,  5, 0, 2, 0, 1)
  call check("Uyzz ", Uyzz,  5, 0, 1, 2, 0)
  call check("Uyzt ", Uyzt,  5, 0, 1, 1, 1)
  call check("Uytt ", Uytt,  5, 0, 1, 0, 2)
  call check("Uzzz ", Uzzz,  5, 0, 0, 3, 0)
  call check("Uzzt ", Uzzt,  5, 0, 0, 2, 1)
  call check("Uztt ", Uztt,  5, 0, 0, 1, 2)
  call check("Uttt ", Uttt,  5, 0, 0, 0, 3)

  call check("Uxxxx", Uxxxx, 5, 4, 0, 0, 0)
  call check("Uxxxy", Uxxxy, 5, 3, 1, 0, 0)
  call check("Uxxxz", Uxxxz, 5, 3, 0, 1, 0)
  call check("Uxxxt", Uxxxt, 5, 3, 0, 0, 1)
  call check("Uxxyy", Uxxyy, 5, 2, 2, 0, 0)
  call check("Uxxyz", Uxxyz, 5, 2, 1, 1, 0)
  call check("Uxxyt", Uxxyt, 5, 2, 1, 0, 1)
  call check("Uxxzz", Uxxzz, 5, 2, 0, 2, 0)
  call check("Uxxzt", Uxxzt, 5, 2, 0, 1, 1)
  call check("Uxxtt", Uxxtt, 5, 2, 0, 0, 2)
  call check("Uxyyy", Uxyyy, 5, 1, 3, 0, 0)
  call check("Uxyyz", Uxyyz, 5, 1, 2, 1, 0)
  call check("Uxyyt", Uxyyt, 5, 1, 2, 0, 1)
  call check("Uxyzz", Uxyzz, 5, 1, 1, 2, 0)
  call check("Uxyzt", Uxyzt, 5, 1, 1, 1, 1)
  call check("Uxytt", Uxytt, 5, 1, 1, 0, 2)
  call check("Uxzzz", Uxzzz, 5, 1, 0, 3, 0)
  call check("Uxzzt", Uxzzt, 5, 1, 0, 2, 1)
  call check("Uxztt", Uxztt, 5, 1, 0, 1, 2)
  call check("Uxttt", Uxttt, 5, 1, 0, 0, 3)
  call check("Uyyyy", Uyyyy, 5, 0, 4, 0, 0)
  call check("Uyyyz", Uyyyz, 5, 0, 3, 1, 0)
  call check("Uyyyt", Uyyyt, 5, 0, 3, 0, 1)
  call check("Uyyzz", Uyyzz, 5, 0, 2, 2, 0)
  call check("Uyyzt", Uyyzt, 5, 0, 2, 1, 1)
  call check("Uyytt", Uyytt, 5, 0, 2, 0, 2)
  call check("Uyzzz", Uyzzz, 5, 0, 1, 3, 0)
  call check("Uyzzt", Uyzzt, 5, 0, 1, 2, 1)
  call check("Uyztt", Uyztt, 5, 0, 1, 1, 2)
  call check("Uyttt", Uyttt, 5, 0, 1, 0, 3)
  call check("Uzzzz", Uzzzz, 5, 0, 0, 4, 0)
  call check("Uzzzt", Uzzzt, 5, 0, 0, 3, 1)
  call check("Uzztt", Uzztt, 5, 0, 0, 2, 2)
  call check("Uzttt", Uzttt, 5, 0, 0, 1, 3)
  call check("Utttt", Utttt, 5, 0, 0, 0, 4)
  print *, 'check_test_4d.f90 finished, max maxerror=',gmaxerror
end program check_test_4d