function c_trapz2d_sample(func,dhx,dhy,xrange,yrange) result(int)
complex(8),dimension(:,:) :: func
real(8),optional :: dhx,dhy
real(8),dimension(2),optional :: xrange,yrange
real(8) :: hx,hy
real(8),dimension(2) :: xrange_,yrange_
integer :: Nx,Ny,i,j
complex(8) :: int
Nx=size(func,1)
Ny=size(func,2)
xrange_=[0d0,1d0];if(present(xrange))xrange_=xrange
yrange_=[0d0,1d0];if(present(yrange))yrange_=yrange
hx=xrange_(2)-xrange_(1);hx=hx/Nx
if(present(dhx))hx=dhx
hy=yrange_(2)-yrange_(1);hy=hy/Ny
if(present(dhy))hy=dhy
!
int = func(1,1)+func(1,Ny)+func(Nx,1)+func(Nx,Ny)
do i=2,Nx
do j=2,Ny
int=int+4d0*func(i,j)
enddo
enddo
do j=2,Ny
int=int+2d0*( func(1,j) + func(Nx,j) )
enddo
do i=2,Nx
int=int+2d0*( func(i,1) + func(i,Ny) )
enddo
int=int*hx*hy/4d0
end function c_trapz2d_sample