
	program test

	parameter(n=50)

	real, dimension(n+1,n+1) :: Dfd, Dfo, Dch
	real, dimension(n+1,n+1)   :: f,dfa,dfn,x,y
	real, dimension(n+1) :: c,xc

	pi=4*atan(1.)

c Initialise coordinates x and y

	do i=0,n
	do j=0,n
	x(i+1,j+1)=cos(i*pi/n)
	xc(i+1)   =cos(i*pi/n)          ! for diff matrix
	y(i+1,j+1)=cos(j*pi/n)
	enddo
	enddo

c Chebyshev differentiation matrix

	c=1
	c(1)=2
	c(n+1)=2

	do i=0,n
	do j=0,n
	if(i==0 .and. j==0) Dch(i+1,j+1)= (2.*n**2+1)/6.
	if(i==n .and. j==n) Dch(i+1,j+1)=-(2.*n**2+1)/6.
	if(i==j .and. j.ne.0 .and.  j.ne.n) 
     &Dch(i+1,j+1)=-xc(j+1)/(2.*(1.-xc(j+1)**2))
	if(i.ne.j) 
     &Dch(i+1,j+1)=(c(i+1)*(-1.)**(i+j))/(c(j+1)*(xc(i+1)-xc(j+1)))
	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(Dch,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(Dch))

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


