common/aaa/h
        write (*,*) 'Enter xl,yl,nx,ny:'
        read (*,*) xl,yl,nx,ny
        dx=xl/nx
        dy=yl/ny
        h=dx
        write (*,*) 'Enter jx,jy for source point:'
        read (*,*) jx,jy
        sum=0.
        do ix=1,nx+1
        do iy=1,ny+1
        factor=h*h
        if(ix.eq.1 .or. ix.eq.(nx+1))factor=factor/2.
        if(iy.eq.1 .or. iy.eq.(ny+1))factor=factor/2.
        xx=(ix-jx)*h
        yy=(iy-jy)*h
        sum=sum+factor/sqrt(xx**2+yy**2+h**2)
        enddo
        enddo
        write (*,*) sum
        a1=(1-jx)*h
        a2=((nx+1)-jx)*h
        b1=(1-jy)*h
        b2=((ny+1)-jy)*h
        call exact(a1,b1,a2,b2)
        end

        subroutine exact(a1,b1,a2,b2)
        common/aaa/h
        ans=-q(a1,b1)+q(a2,b1)+q(a1,b2)-q(a2,b2)
        write (*,*) ans
        end

        real function q(r1,r2)
        common/aaa/h
        d=sqrt(r1**2+r2**2+h**2)
        q=h*atan((r1*r2)/(h*d))-r2*alog(r1+d)-r1*alog(r2+d)
        return
        end