
	program test

	parameter(n=51)

	real, dimension(n,n) :: Dfd, Dfo, Dch
	real, dimension(n,n)   :: f,dfa,dfn,x,y


c Fourier differentiation matrix

	pi=4*atan(1.)

      do i=1,n
	do j=1,n
	if(i.ne.j)Dfo(i,j)=pi*(-1)**(i-j)/(2*sin(pi*(i-j)/n))
	enddo
	enddo



c Initialise coordinates x and y


	dx=2./(n-1)
	do i=0,n-1
	do j=0,n-1
	x(i+1,j+1)=-1+i*dx
	y(i+1,j+1)=-1+j*dx
	enddo
	enddo

c Read half width of Gaussian test function
	
	write(*,*)' Give half-width of Gauss function '
	read(*,*)a

	
c calculate test function f 	

	f = exp(-1./a**2*(x*x+y*y))

c ------------------------------------------------------

c calculate derivative of test function - d/dx

	dfa=-2./a**2*x*exp(-1/a**2*(x*x+y*y))
	dfn=matmul(Dfo,f)

c ... and compare min/max and some points with the analytical solution

	write(*,*)minval(f),' < f < ',maxval(f)
	write(*,*)minval(dfa),' < df (ana) < ',maxval(dfa)
	write(*,*)minval(dfn),' < dfn (num) < ',maxval(dfn)

c compare some central points

	write(*,*)' Analytical v.s. Numerical Derivative (x)'
	do i=n/2-5,n/2+5
	write(*,*)dfa(i,n/2),'   ',dfn(i,n/2)
	enddo
	write(*,*)' '

c -------------------------------------------------------


c calculate derivative of test function - d/dy

	dfa=-2./a**2*y*exp(-1/a**2*(x*x+y*y))
	dfn=matmul(f,transpose(Dfo))

c ... and compare min/max and some points with the analytical solution

	write(*,*)minval(f),' < f < ',maxval(f)
	write(*,*)minval(dfa),' < df (ana) < ',maxval(dfa)
	write(*,*)minval(dfn),' < dfn (num) < ',maxval(dfn)

c compare some central points

	write(*,*)' Analytical v.s. Numerical Derivative (y)'
	do i=n/2-5,n/2+5
	write(*,*)dfa(n/2,i),'   ',dfn(n/2,i)
	enddo
	write(*,*)' '
	end program test


