program EXPRD

c     Example for 2D potential problem with periodic conditions
      
      integer N1,N2,BCT(1000),N,i,ians
      
      double precision xb(1000),yb(1000),xm(1000),ym(1000),
     & nx(1000),ny(1000),lg(1000),BCV(1000),
     & phi(1000),dphi(1000),pint,dl,xi,eta,pi
     
      print*,'Enter no. of elements on a vertical side:'
      read*,N1
      print*,'Enter no. of elements on a horizontal side:'
      read*,N2
      N=2*N1+2*N2

c     Putting N well spaced out points in anticlockwise direction on the square boundary
      
	pi=4d0*datan(1d0)

      dl=0.5d0/dfloat(N1)
      
      do 10 i=1,N1
      xb(i)=0d0
      yb(i)=0.50d0-dfloat(i-1)*dl
      xb(N1+N2+i)=2d0*pi
      yb(N1+N2+i)=dfloat(i-1)*dl 
10    continue

      dl=2d0*pi/dfloat(N2)
      
      do 12 i=1,N2
      xb(N1+i)=dfloat(i-1)*dl
      yb(N1+i)=0d0
      xb(N1+N2+N1+i)=2d0*pi-dfloat(i-1)*dl
      yb(N1+N2+N1+i)=0.5d0
12    continue

      xb(N+1)=xb(1)
      yb(N+1)=yb(1)
      
c     Working out the midpoints, lengths and normal vectors of the elements      
      
      do 20 i=1,N
      xm(i)=0.5d0*(xb(i)+xb(i+1))
      ym(i)=0.5d0*(yb(i)+yb(i+1))
      lg(i)=dsqrt((xb(i+1)-xb(i))**2d0+(yb(i+1)-yb(i))**2d0)
      nx(i)=(yb(i+1)-yb(i))/lg(i)
      ny(i)=(xb(i)-xb(i+1))/lg(i)
20    continue
 
c     Specifying the boundary conditions on the sides of the square  
 
      do 30 i=1,N
      if ((i.gt.N1).and.(i.le.(N1+N2))) then
      BCT(i)=0
      BCV(i)=dsin(xm(i))
      else if (i.gt.(2*N1+N2)) then
      BCT(i)=1
      BCV(i)=dexp(0.5d0)*dsin(xm(i))
      else            
      BCT(i)=2
      endif
30    continue

      call CELAPPER(N,N1,N2,xm,ym,xb,yb,nx,ny,lg,BCT,BCV,phi,dphi) 
      
50    print*,'Enter coordinates xi and eta of an interior point:'
      
      read*,xi,eta
      
      call CELAP2(N,xi,eta,xb,yb,nx,ny,lg,phi,dphi,pint)               
 
      write(*,60)pint,dexp(eta)*dsin(xi)

60    format('Numerical and exact values are:',
     & F14.6,' and',F14.6,' respectively')
     
      print*,'To continue with another point enter 1:'
      read*,ians
      
      if (ians.eq.1) goto 50
        
      end
     
c     Copyright 2006-2007 WT Ang, A Beginner's Course in Boundary Element Methods