c     bem.f(program 1)

c     program for solution of two dimensional potential problems by the
c     b.i.e. method with constant elements

      common n,l,lec,imp
      dimension x(41),y(41),xm(40),ym(40),g(40,40),fi(40),dfi(40)
      dimension kode(40),cx(40),cy(40),sol(40),h(40,40)
      character*80 filea

c     initialization of program parameters
c     nx = maximum dimension of the system of equations.

      nx = 40

c     assign data set numbers for input, lec, and output, imp

      lec = 5
      imp = 6
      write(*,*) 'Enter the name of the input data file'
      read(*,'(a)') filea

      open(10,file=filea, status='old')

c     input

      call input(cx,cy,x,y,kode,fi)

c     form system of equations

      call fmat(x,y,xm,ym,g,h,fi,dfi,kode,nx)

c     solution of the system of equations

      call slnpd (g,dfi,d,n,nx)

c     compute the potential values in internal points

      call inter (fi,dfi,kode,cx,cy,x,y,sol)

c     output

      call outpt(xm,ym,fi,dfi,cx,cy,sol)
      stop
      end

      subroutine input (cx,cy,x,y,kode,fi)

c     program 2

      common n,l,lec,imp
      dimension cx(1),cy(1),x(1),y(1),kode(1),fi(1)
c      dimension title(18)

c     n = number of coundary elements
c     l = number of internal points where the function is calculated

      write (imp,100)
 100  format(' ',79('*'))

c     read name of the job

c      read (lec,150) title
c 150  format(18a4)
c      write(imp,250) title
c 250  format (25x, 18a4)

c     read basic parameters

c      write(*,*) 'enter n and l'
      read (10,*) n,l
c 200  format(2i5)
      write(imp,300) n,l
 300  format(//' data'//2x,'number of boundary elements =',
     1 i3/2x,'# of internal pts where the function is calculated =',i3)

c     read internal points coordinates

c	write(*,*) 'Enter internal point coordinates'
      do 1 i = 1,l
 1       read(10,*) cx(i),cy(i)
c 400     format(2f10.4)

c     read coordinates of extreme points of the boundary elements 
c     in array x and y

         write(imp,500)
 500     format(//2x,'coordinates of the extreme points of the 
     1 boundary elements',//4x,'point',10x,'x',18x,'y')

c	write(*,*) 'Enter extreme points of boundary elements'
         do 10 i = 1,n
            read (10,*) x(i),y(i)
c 600        format(2f10.4)
 10         write(imp,700) i,x(i),y(i)
 700        format (5x,i3,2(5x,e14.7))

c     read boundary conditions
c     fi(i) = value of the potential in the node i if kode = 0,
c     value of the potential derivative if kode = 1.

            write(imp,800)
 800        format(//2x,'boundary conditions'//5x,'node',6x,'code',
     1           5x,'prescribed value')

c	write(*,*) 'Enter boundary conditions'
            do 20 i = 1,n
               read(10,*) kode(i),fi(i)
c 900           format(i5,f10.4)
 20            write (imp,950) i,kode(i),fi(i)
 950           format(5x,i3,8x,i1,8x,e14.7)
               return
               end
      subroutine fmat (x,y,xm,ym,g,h,fi,dfi,kode,nx)

c     program 3

c     this subroutine computes g and h matrices and form the system ax = f

      common n,l,lec,imp
      dimension x(1), y(1),xm(1),ym(1),g(nx,nx),h(nx,nx),fi(1),kode(1)
      dimension dfi(1)

c     compute the mid-point coordinates and store in array xm and ym

      x(n+1) = x(1)
      y(n+1) = y(1)
      do 10 i = 1,n
         xm(i) = (x(i) + x(i+1))/2
 10      ym(i) = (y(i) + y(i+1))/2

c     compute g and h matrices

         do 30 i = 1,n
            do 30 j = 1,n
               if(i-j) 20,25,20
 20            call inte(xm(i),ym(i),x(j),y(j),x(j+1),
     1             y(j+1),h(i,j),g(i,j))
               go to 30
 25            call inlo(x(j),y(j),x(j+1),y(j+1),g(i,j))
               h(i,j) = 3.1415926
 30            continue

c     arrange the system of equations ready to be solved

c               do 50 j = 1,n
c                  if(kode(j)) 50,50,40
c 40               do 50 i = 1,n
c                     ch = g(i,j)
c                     g(i,j) = -h(i,j)
c                     h(i,j) = -ch
c 50               continue

               do 50 j = 1,n
                  if(kode(j).gt.0)then
                     do 40 i = 1,n
                        ch = g(i,j)
                        g(i,j) = -h(i,j)
                        h(i,j) = -ch
 40                  continue
                  else
                     continue
                  endif
 50            continue

c     dfi originally contains the independent coefficients 
c     after solution it will contain the values of the system unknowns

                     do 60 i = 1,n
                        dfi(i) = 0
                        do 60 j = 1,n
                           dfi(i) = dfi(i) + h(i,j) * fi(j)
 60                         continue
                            return
                            end

      subroutine inte(xp,yp,x1,y1,x2,y2,h,g)

c     program 4

c     this subroutine computes the values of the h and g matrix
c     off diagonal elements by means of numberical integration
c     along the boundary elements

c     dist = distance from the point under consideration to the 
c     boundary elements
c     ra = distance from the point under consideration to the 
c     integration points in the boundary elements

      dimension xco(4), yco(4),gi(4),ome(4)
      gi(1) = 0.86113631
      gi(2) = -gi(1)
      gi(3) = 0.33998104
      gi(4) = -gi(3)
      ome(1) = 0.34785485
      ome(2) = ome(1)
      ome(3) = 0.65214515
      ome(4) = ome(3)
      ax = (x2 - x1)/2
      bx = (x2 + x1)/2
      ay = (y2 - y1)/2
      by = (y2 + y1)/2
      if(ax)10,20,10
 10   ta = ay/ax
      dist = abs((ta*xp - yp+y1 - ta*x1)/sqrt(ta**2+1))
      go to 30
 20   dist = abs(xp - x1)
 30   sig = (x1 - xp)*(y2-yp) - (x2 - xp)*(y1 - yp)
      if(sig)31,32,32
 31   dist = -dist
 32   g = 0.
      h = 0.
      do 40 i = 1,4
         xco(i) = ax*gi(i) + bx
         yco(i) = ay*gi(i) + by
         ra = sqrt((xp - xco(i))**2 + (yp - yco(i))**2)
         g = g + alog(1/ra)*ome(i)*sqrt(ax**2 + ay**2)
 40       h = h - (dist*ome(i)*sqrt(ax**2 + ay**2)/ra**2)
          return
          end



      subroutine inlo(x1,y1,x2,y2,g)

c     program 5

c     this subroutine computes the values of the diagonal
c     elements of the g matrix

      ax = (x2 - x1)/2
      ay = (y2 - y1)/2
      sr = sqrt(ax**2 + ay**2)
      g = 2*sr*(alog(1/sr)+1)
      return
      end


      subroutine slnpd(a,b,d,n,nx)
c     program 6

c     solution of linear systems of equations
c     by the gauss elimination method providing
c     for interchanging rows when encountering a 
c     zero diagonal coefficient

c     a : system matrix
c     b : originally it contains the independent coefficients.
c         after solution it contains the values of the system unknowns.

c     n : actual number of unknowns
c     nx : row and column dimension of a

      dimension a(nx,nx),b(nx)
      n1 = n-1
      do 100 k = 1,n1
         k1 = k+1
         c = a(k,k)
         if(abs(c)-0.000001)1,1,3
 1       do 7 j = k1,n

c     try to interchange rows to get non zero diagonal coefficient

            if(abs(a(j,k)) - 0.000001)7,7,5
 5          do 6 l = k,n
               c = a(k,l)
               a(k,l) = a(j,l)
 6             a(j,l) = c
               c = b(k)
               b(k) = b(j)
               b(j) = c
               c = a(k,k)
               go to 3
 7             continue
 8             write(6,2)k
 2             format('*** singularity in row',i5)
               d = 0.
               go to 300

c     divide row by diagonal coefficient

 3              c = a(k,k)
                do 4 j = k1,n
 4                 a(k,j) = a(k,j)/c
                   b(k) = b(k)/c

c     eliminate unknown x(k) from row i
                   do 10 i = k1,n
                      c = a(i,k)
                      do 9 j = k1,n
 9                       a(i,j) = a(i,j) - c*a(k,j)
 10                      b(i) = b(i) - c*b(k)
 100                     continue

c     compute last unknown

                         if ((abs(a(n,n)) - 0.000001).lt.0) then
                            write(6,22)k
 22                      format('*** singularity in row',i5)
                      else
                            continue
                         endif
 101                     b(n) = b(n)/a(n,n)

c     apply backsubstitution process to compute remaining unknowns

                         do 200 l = 1,n1
                            k = n - l
                            k1 = k+1
                            do 200 j = k1,n
 200                        b(k) = b(k) - a(k,j)*b(j)

c     compute value of determinant

                            d = 1.
                            do 250 i = 1,n
 250                           d = d*a(i,i)
 300                           return
                               end

      subroutine inter(fi,dfi,kode,cx,cy,x,y,sol)
c     program 7

c     this subroutine computes the potential value for internal points.

      common n,l,lec,imp
      dimension fi(1),dfi(1),kode(1),cx(1),cy(1),x(1),y(1),sol(1)

c     reorder fi and dfi array to put all the values of the potential
c     in fi and all the values of the derivative in dfi

      do 20 i = 1,n
         if (kode(i)) 20,20,10
 10      ch = fi(i)
         fi(i) = dfi(i)
         dfi(i) = ch
 20      continue

c     compute the potential values for internal points

         do 40 k = 1,l
            sol (k) = 0
            do 30 j = 1,n
               call inte(cx(k),cy(k),x(j),y(j),x(j+1),y(j+1),a,b)
 30            sol(k) = sol(k) + dfi(j)*b - fi(j)*a
 40            sol(k) = sol(k)/(2*3.1415926)
               return
               end

      subroutine outpt(xm,ym,fi,dfi,cx,cy,sol)
c     program 8

      common n,l,lec,imp
      dimension xm(1),ym(1),fi(1),dfi(1),cx(1),cy(1),sol(1)
      write(imp,100)
 100  format(' ',79('*')//1x,'results'//2x,'boundary nodes'//9x,
     1 'x',16x,'y',12x,'potential',10x,'potential derivative'/)
      do 10 i = 1,n
 10      write (imp,200) xm(i),ym(i),fi(i),dfi(i)
 200     format(2x,e14.7,3(6x,e14.7))
         write(imp,300)
 300     format(//,2x,'internal points',//11x,'x',18x,'y',14x,
     1   'potential',/)
       do 20 k = 1,l
 20       write(imp,400) cx(k),cy(k),sol(k)
 400      format(3(5x,e14.7))
          write(imp,500)
 500      format(' ',79('*'))
          return
          end