function c_simps2d_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/2
if(present(dhx))hx=dhx
hy=yrange_(2)-yrange_(1);hy=hy/Ny/2
if(present(dhy))hy=dhy
!
int = func(1,1)+func(1,Ny)+func(Nx,1)+func(Nx,Ny)
!
do j=1,Ny
int=int+4d0*( func(1,2*j-1) + func(Nx,2*j-1) )
enddo
do j=1,Ny-1
int=int+2d0*( func(1,2*j) + func(Nx,2*j) )
enddo
!
do i=1,Nx
int=int+4d0*( func(2*i-1,1) + func(2*i-1,Ny) )
enddo
do i=1,Nx-1
int=int+2d0*( func(2*i,1) + func(2*i,Ny) )
enddo
!
do j=1,Ny
do i=1,Nx
int=int+16d0*func(2*i-1,2*j-1)
enddo
enddo
do j=1,Ny-1
do i=1,Nx
int=int+8d0*func(2*i-1,2*j)
enddo
enddo
!
do j=1,Ny
do i=1,Nx-1
int=int+8d0*func(2*i,2*j-1)
enddo
enddo
do j=1,Ny-1
do i=1,Nx-1
int=int+4d0*func(2*i,2*j)
enddo
enddo
int=int*hx*hy/9d0
end function c_simps2d_sample