implicit real*8 (a-h,o-z)

        parameter(a= -1.0d0,b=1.0d0,c=-1.0d0,d=1.0d0,
     1      m=80,n=80,n1=160)

        dimension x1(n1+1),y1(n1+1),f(0:m,0:n),x(0:m),y(0:n)
     8        ,work(2*(m+1)*(n+1)),tru(m,0:n)

        h= (b -a)/m
        do i=0,m
          x(i) = a + float(i)*h
	enddo

        do j=0,n
           y(j) = c + float(j)*h
	enddo

C--------- prescribe the sourse term

	do i=0,m
	  do j=0,n
	     f(i,j) = 0.0
	  enddo
	enddo

C--------- prescribe the boundary condition

	do i=0,m
	   f(i,0) = 0.0
	   f(i,n) = 0.0
	enddo
	
	do j=0,m
           f(0,j) = 0.0
           f(m,j) = 0.0
        enddo


C---------- Get the control points on the interface ------------------

        pi = datan(1.0d0)*4.0
        dh = 2.0d0*pi/n1

        do i=1,n1
          t = dh*(i-1.0)
          x1(i) = 0.4 * dcos(t)
          y1(i) = 0.5 * dsin(t)
	enddo
	x1(n1+1) = x1(1)
	y1(n1+1) = y1(1)

C--------- Peskin discrete delta function

	do j=1,n-1
	   do i=1,m-1
	      do k=1,n1
		dsk = dsqrt( (x1(k+1)-x1(k))**2 + (y1(k+1)-y1(k))**2 )
		f(i,j) = f(i,j) + dsk * delcos(x(i)-x1(k),h)*
     1                   delcos(y(j)-y1(k),h)
	      enddo	      
	   enddo
	enddo

C------ Call fast Poisson solver from Fishpack

        idimf = m+1
        mbdcnd = 1
        nbdcnd = 1
	elmbda = 0.0

        call hwscrt(a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,
     1                   elmbda,f,idimf,pertrb,ierror,work)

C------- Output solution 

	open(412,file='do.m',status='unknown')
        do j=0,n
          write(412,100)(f(i,j),i=0,m)
        enddo
        close(412)

 100    format(161e13.4)

	stop
	end