subroutine CELAPPER(N,N1,N2,xm,ym,xb,yb,nx,ny,lg,BCT,BCV,phi,dphi)

c     BVP governed by 2D Laplace's equation with periodic BCs
c     Determining unknowns on the boundary by constant elements      
c     N: number of elements
c     N1: no. of elements on each vertical side
c     N2: no. of elements on each horizontal side      
c     xm,ym: real arrays containing midpoints of elements     
c     xb,yb: real arrays containing starting points of elements
c     nx,ny: real arrays containing normal vectors to elements
c     lg: real arrays containing lengths of elements
c     BCT: integer arrays containing type of boundary conditions
c     BCT(k)=0 if phi is specified on k-th element
c     BCT(k)=1 if normal derivative of phi is specified on k-th element
c     BCV: real arrays containing values specified on elements
c     phi: real arrays containing values of phi on elements (returned)
c     dphi: real arrays containing values of normal derivative of phi on elements (returned)

c     For more details, read document "2D potential problems with periodic BCs"
      
      integer m,k,N,N1,N2,BCT(1000),k1
      
      double precision xm(1000),ym(1000),xb(1000),yb(1000),
     & nx(1000),ny(1000),lg(1000),BCV(1000),A(1000,1000),
     & B(1000),pi,PF1,PF2,del,phi(1000),dphi(1000),F1,F2,
     & Z(1000)
     
      pi=4d0*datan(1d0)

	do 2 m=1,N+2*N1
	B(m)=0d0
	do 2 k=1,N+2*N1
	A(m,k)=0d0
2	continue
     
      do 10 m=1,N
	k1=0
      B(m)=0d0
      do 5 k=1,N
      call CPF(xm(m),ym(m),xb(k),yb(k),nx(k),ny(k),lg(k),PF1,PF2)
      F1=PF1/pi
      F2=PF2/pi
      if (k.eq.m) then
      del=1d0
      else
      del=0d0
      endif
      if (BCT(k).eq.0) then
      A(m,k)=-F1
      B(m)=B(m)+BCV(k)*(-F2+0.5d0*del)
      else if (BCT(k).eq.1) then
      A(m,k)=F2-0.5d0*del
      B(m)=B(m)+BCV(k)*F1
      else
	k1=k1+1
      A(m,k)=F2-0.5d0*del
      A(m,N+k1)=-F1
      endif      
5     continue      
10    continue

	do 12 k=1,N1
	A(N+k,k)=1d0
	A(N+k,2*N1+N2+1-k)=-1d0
	A(N+N1+k,N+k)=1d0
	A(N+N1+k,N+2*N1+N2+1-k-N2)=1d0
12    continue


      call solver(A,B,N+2*N1,1,Z)


	k1=0
      
      do 15 m=1,N
      if (BCT(m).eq.0) then
      phi(m)=BCV(m)
      dphi(m)=Z(m)
      else if (BCT(m).eq.1) then
      phi(m)=Z(m)
      dphi(m)=BCV(m)
      else
	k1=k1+1
	phi(m)=Z(m)
      dphi(m)=Z(N+k1)
      endif
15    continue      

      return
      end      
      
c     Copyright 2006-2007 WT Ang, A Beginner's Course in Boundary Element Methods