Program Potential
      Implicit NONE
c        real PI
c        parameter (PI=3.141593)
        real x_min /-3.0/, x_max /3.0/
        real y_min /-3.0/, y_max /3.0/
        real radius /0.2/      
        real step /0.1/
        real x1(2) / -1.0, 0.0/
        real x2(2) /  1.0, 0.0/
        Integer i,j,imax,jmax
c function:
c local:
        real x(2)
        real phi
        imax=(x_max-x_min)/step
        jmax=(y_max-y_min)/step
        Do i=1,imax
         Do j=1,jmax
           x(1)=x_min+real(i-1)*step
           x(2)=y_min+real(j-1)*step
           call calc_phi(x1,x2,x,radius,phi)
           write(*,' (3f10.4)') x(1), x(2), phi
         Enddo
         write(*,*) ' '
        Enddo
       end
       subroutine calc_phi(x1,x2,x,radius,phi)
        implicit NONE
        real x1(2)
        real x2(2)
        real x(2)
        real phi,radius
        real r1,r2
c        real r1,r2,r1radius,r2radius
          r1=sqrt((x(1)-x1(1))**2 + (x(2)-x1(2))**2)
          r2=sqrt((x(1)-x2(1))**2 + (x(2)-x2(2))**2)
c             r1radius=sqrt((-radius-x1(1))**2 + (0.-x1(2))**2)
c             r2radius=sqrt((radius-x2(1))**2 + (0.-x2(2))**2)
          if(r1.lt. radius )then
             phi=1./radius
          else if(r2.lt. radius )then
             phi=-1./radius
          else
             phi=1./r1 -1./r2
          Endif
        return
        end