C***************************************************************
C           Subroutine LIBEM2 by Stephen Kirkup                *        
C***************************************************************
C 
C  Version 2.2 (September 2015) Quadrature rule tidied up.
C   (August 2015). Linked to new guide. Versions 2
C   (July 2015) New GLS routine implemented and simplified 
C   integral equation formulation (first version 2001) 
C  School of Engineering 
C  University of Central Lancashire - www.uclan.ac.uk 
C  smkirkup@uclan.ac.uk
C  http://www.researchgate.net/profile/Stephen_Kirkup
C
C  This open source code can be found at
C   www.boundary-element-method.com/fortran/LIBEM2.FOR 
C  The guide to using the subroutine is found at
C   www.boundary-element-method.com/fortran/LIBEM2_FOR.pdf 
C
C  Issued under the GNU General Public License 2007, see gpl.txt
C
C  Part of the the author's open source BEM packages. 
C  All codes and manuals can be downloaded from 
C  www.boundary-element-method.com
C
C***************************************************************
C
C How to use the subroutine
C -------------------------
C
C The following diagram shows how the subroutine is to be used. A main
C program is required along with the definition of FNHANK, a routine
C that must return the spherical hankel function.
C
C                                   ....................................
C                                   :                                  :
C                                   :                                  :
C      ----------------------       :     --------------------------   :
C      |                    |       :     |                        |   :
C      |   MAIN PROGRAM     |------->-----|      LIBEM2            |   :
C      | (e.g. libem2_t.for)|       :     |                        |   :
C      |                    |       :     --------------------------   :
C      ----------------------       :                 |                :
C                                   :                 >                :
C                                   :                 |                :
C                                   :      ------------------------    :
C          Package ---------------->:      | subordinate routines |    :
C                                   :      ------------------------    :
C                                   :                                  :
C                                   :      (this file)                 :  
C                                   :..................................:
C                                  /         |                 |
C                               |_           >                 >
C                              /             |                 |
C             ................       ................   ................  
C             :              :       :   --------   :   :  --------    : 
C             : (geom2d.for) :---<---:   | L2LC |   :   :  | GLS  |    : 
C             :              :       :   --------   :   :  --------    :  
C             :..............:       : -------------:   : -------------:  
C                                    : |subordinate|:   : |subordinate|: 
C                                    : | routines  |:   : | routines  |:  
C                                    : -------------:   : -------------: 
C                                    :              :   :              : 
C                                    : (l2lc.for)   :   :  (gls.for)   :
C                                    :..............:   :..............:
C                                    
C
C The contents of the main program must be linked to LIBEM2.FOR, L2LC.FOR 
c  and GLS.FOR.
C
C Method of solution
C ------------------
C 
C In the main program, the boundary must be described as a set of
C  panels. The panels are defined by two indices (integers) which
C  label a node or vertex on the boundary. The data structure
C  VERTEX lists and enumerates the coordinates of the vertices, the
C  data structure SELV defines each panel by indicating the labels
C  for the two nodes it lies between and hence enumerates the panels.
C The boundary solution points (the points on the boundary at which 
C  {\phi} (SPHI) and d {\phi}/dn (SVEL) are returned) are at the centres
C  of the panels. The boundary functions {\alpha} (SALPHA), {\beta} 
C  (SBETA) and f (SF) are also defined by their values at the centres
C  of the panels.
C Normally a solution in the domain is required. By listing the 
C  coordinates of all the interior points in PINT, the subroutine
C  returns the value of {\phi} at these points in PDPHI.

C Notes on the geometric parameters
C ---------------------------------
C (1) Each of the vertices listed in VERTEX must be distinct points
C  with respect to EGEOM.
C (2) The boundary must be complete and closed. Thus 
C  SELV(i,2)=SELV(i+1,1) for i=1..NSE-1 and SELV(1,1)=SELV(NSE,2).
C (3) The indices of the nodes listed in SELV must be such that they
C  are ordered counter-clockwise around the boundary.
C (4) The largest panel must be no more than 10x the length of the
C  smallest panel.

C Notes on the interior points 
C ----------------------------
C (1) The points in PINT should lie within the boundary, as defined
C  by the parameters VERTEX and SELV. Any point lying outside the 
C  boundary will return a corresponding value in PDPHI that is near
C  zero.

C Notes on the boundary condition
C -------------------------------
C (1) For each i=1..NSE, it must not be the case that both of SALPHA(i)
C  and SBETA(i) are zero.

C External modules in external files
C ==================================
C subroutine L2LC: Returns the individual discrete Laplace integral
C  operators (in file L2LC.FOR).
C subroutine GLS: Solves a general linear system of equations. 
C  (in file GLS.FOR)
C real function DIST2: Returns the distance between two 2-vectors. (in
C  file GEOM2D.FOR)

C External modules provided in the package (this file)
C ====================================================
C subroutine GL8: Returns the points and weights of the 8-point Gauss-
C  Legendre quadrature rule.
C real function FNSQRT(X): real X : Returns the square root of X.
C real function FNLOG(X): real X : Returns the natural logarithm of X.


C The subroutine

      SUBROUTINE LIBEM2(MAXNV,NV,VERTEX,MAXNSE,NSE,SELV,
     *                 MAXNPI,NPI,PINT,
     *                 SALPHA,SBETA,SF,SIPHI,PDIPHI,
     *                 LSOL,LVALID,EGEOM,
     *                 SPHI,SVEL,PDPHI,
     *                 L_SS,M_SSPLUS,L_DS,M_DS,
     *                 C,PERM,XORY,WKSPC1,WKSPC2,WKSPC3)
      PARAMETER (MAXNQ=32)

C  Boundary geometry
C   Limit on the number of vertices on S
      INTEGER    MAXNV
C   The number of vertices on S
      INTEGER    NV
C   The coordinates of the vertices on S
      REAL*8     VERTEX(MAXNV,2)
C   Limit on the number of panels describing S
      INTEGER    MAXNSE
C   The number of panels describing S
      INTEGER    NSE
C   The indices of the vertices describing each panel
      INTEGER    SELV(MAXNSE,2)
      
C  Interior points at which the solution is to be observed
C   Limit on the number of points interior to the boundary where 
C    solution is sought
      INTEGER    MAXNPI
C   The number of interior points
      INTEGER    NPI
C   Coordinates of the interior points
      REAL*8     PINT(MAXNPI,2)

C  The boundary condition is such that {\alpha} {\phi} + {\beta} v = f
C  where alpha, beta and f are real valued functions over S.
C  The functions are set values at the central points.
C   function alpha
      REAL*8     SALPHA(MAXNSE)
C   function beta
      REAL*8     SBETA(MAXNSE)
C   function f
      REAL*8     SF(MAXNSE)

C  The incident potential on S
      REAL*8     SIPHI(MAXNSE)
C  The incident potential at the chosen interior points
      REAL*8     PDIPHI(MAXNPI)

C  Validation and control parameters
      LOGICAL    LSOL
      LOGICAL    LVALID
      REAL*8     EGEOM

C  Solution 
C   function phi
      REAL*8     SPHI(MAXNSE)
C   function vel
      REAL*8     SVEL(MAXNSE)
C   domain solution
      REAL*8     PDPHI(MAXNPI)

C  Other information 
      REAL*8     L_SS(MAXNSE,MAXNSE)
      REAL*8     M_SSPLUS(MAXNSE,MAXNSE)
      REAL*8     L_DS(MAXNPI,MAXNSE)
      REAL*8     M_DS(MAXNPI,MAXNSE)
      REAL*8     C(MAXNSE)
      INTEGER*4  PERM(MAXNSE)
      LOGICAL    XORY(MAXNSE)

C  Work space
      REAL*8     WKSPC1(MAXNSE)
      REAL*8     WKSPC2(MAXNSE)
      REAL*8     WKSPC3(MAXNSE)

C  External function
      REAL*8     DIST2

C  Constants
C   Real scalars: 0, 1, 2, half, pi
      REAL*8     ZERO,ONE,TWO,HALF,PI

C  Geometrical description of the boundary
C   panels counter
      INTEGER    ISE,JSE
C   The points interior to the boundary where the solution is sought 
      INTEGER    IPI
C   Parameters for L2LC
      REAL*8     P(2),PA(2),PB(2),QA(2),QB(2),VECP(2)
      LOGICAL    LPONEL

C   Quadrature rule parameters for L2LC
C    Actual number of quadrature points
      INTEGER    NQ
C    Abscissae of the actual quadrature rule
      REAL*8     AQ(MAXNQ)
C    Weights of the actual quadrature rule
      REAL*8     WQ(MAXNQ)

C  Validation and control parameters for subroutine L2LC
      LOGICAL    LVAL
      REAL*8     EQRULE
      LOGICAL    LFAIL1
      LOGICAL    LL
      LOGICAL    LM
      LOGICAL    LMT
      LOGICAL    LN

C  Parameters for subroutine L2LC. 
      REAL*8     DISL
      REAL*8     DISM
      REAL*8     DISMT
      REAL*8     DISN

C  Other variables
C   Error flag
      LOGICAL    LERROR
C   Failure flag
      LOGICAL    LFAIL
C   Accumulation of solution {\phi}
      REAL*8     SUMPHI
C   Stores a vertex
      REAL*8     QC(2)
C   Maximum,minimum sizes of panels
      REAL*8     SIZMAX,SIZMIN
C   The `diameter' of the boundary or the maximum distance between any
C    two vertices
      REAL*8     DIAM
      REAL*8     SUMM     

C INITIALISATION
C --------------

C Set constants
      ZERO=0.0D0
      ONE=1.0D0
      TWO=2.0D0
      HALF=ONE/TWO
      PI=4.0D0*ATAN(ONE)

      EQRULE=1.0D-6

C Validation
C ==========

C Validation of parameters of LIBEM2
C ---------------------------------

      IF (LVALID) THEN

C Validation of main paramters
        LERROR=.FALSE.
        IF (MAXNV.LT.3) THEN
          WRITE(*,*) 'MAXNV = ',MAXNV
          WRITE(*,*) 'ERROR(LIBEM2) - must have MAXNV>=3'
          LERROR=.TRUE.
        END IF
        IF (NV.LT.3.OR.NV.GT.MAXNV) THEN
          WRITE(*,*) 'NV = ',NV
          WRITE(*,*) 'ERROR(LIBEM2) - must have 3<=NV<=MAXNV'
          LERROR=.TRUE.
        END IF
        IF (MAXNSE.LT.3) THEN
          WRITE(*,*) 'MAXNSE = ',MAXNSE
          WRITE(*,*) 'ERROR(LIBEM2) - must have MAXNSE>=3'
          LERROR=.TRUE.
        END IF
        IF (NSE.LT.3.OR.NSE.GT.MAXNSE) THEN
          WRITE(*,*) 'NSE = ',NSE
          WRITE(*,*) 'ERROR(LIBEM2) - must have 3<=NSE<=MAXNSE'
          LERROR=.TRUE.
        END IF
        IF (MAXNPI.LT.1) THEN
          WRITE(*,*) 'MAXNPI = ',MAXNPI
          WRITE(*,*) 'ERROR(LIBEM2) - must have MAXNPI>=1'
          LERROR=.TRUE.
        END IF
        IF (NPI.LT.0.OR.NPI.GT.MAXNPI) THEN
          WRITE(*,*) 'NPI = ',NPI
          WRITE(*,*) 'ERROR(LIBEM2) - must have 3<=NPI<=MAXNPI'
          LERROR=.TRUE.
        END IF
        IF (EGEOM.LE.ZERO) THEN
          WRITE(*,*) 'NPI = ',NPI
          WRITE(*,*) 'ERROR(LIBEM2) - EGEOM must be positive'
          LERROR=.TRUE.
        END IF
        IF (LERROR) THEN
          LFAIL=.TRUE.
          WRITE(*,*)
          WRITE(*,*) 'Error(s) found in the main parameters of LIBEM2'
          WRITE(*,*) 'Execution terminated'
          STOP
        END IF
      END IF

C Find the diameter DIAM of the boundary
      DIAM=0.0
      DO 100 IV=1,NV-1
        PA(1)=VERTEX(IV,1)
        PA(2)=VERTEX(IV,2)
        DO 110 JV=IV+1,NV
          PB(1)=VERTEX(JV,1)
          PB(2)=VERTEX(JV,2)
          DIAM=MAX(DIAM,DIST2(PA,PB))
110     CONTINUE
100   CONTINUE

      IF (LVALID) THEN
        LERROR=.FALSE.
C Check that the boundary defined by SELV is complete and closed
        DO 120 ISE=1,NSE-1
          IF (SELV(ISE,2).NE.SELV(ISE+1,1)) LERROR=.TRUE.
120     CONTINUE
        IF (LERROR) THEN
          WRITE(*,*) 'ERROR(LIBEM2) - boundary defined by SELVis not'
          WRITE(*,*) ' complete and closed'
        END IF
C Check that EGEOM is not too large
        IF (EGEOM.GT.DIAM/100.0D0) THEN
          WRITE(*,*) 'EGEOM = ',EGEOM
          WRITE(*,*) 'ERROR(LIBEM2) - EGEOM is set too large'
          LERROR=.TRUE.
        END IF
        IF (LERROR) THEN
          LFAIL=.TRUE.
          WRITE(*,*)
          WRITE(*,*) 'Error in boundary geometry or EGEOM'
          WRITE(*,*) 'Execution terminated'
        END IF
      END IF                  

      IF (LVALID) THEN
C Check that the vertices are distinct with respect to EGEOM
        LERROR=.FALSE.
        DO 130 IV=1,NV-1
          PA(1)=VERTEX(IV,1)
          PA(2)=VERTEX(IV,2)
          DO 140 JV=IV+1,NV
            PB(1)=VERTEX(JV,1)
            PB(2)=VERTEX(JV,2)
            IF (ABS(PA(1)-PB(1)).LT.EGEOM) THEN
              IF (ABS(PA(2)-PB(2)).LT.EGEOM) THEN
                WRITE(*,*) 'Vertices ',IV,JV,' are not distinct'
                LERROR=.TRUE.
              END IF
            END IF
140       CONTINUE
130     CONTINUE
        IF (LERROR) THEN
          WRITE(*,*) 
          WRITE(*,*) 'ERROR(LIBEM2) - Vertices (see above) coincide'
          WRITE(*,*) 'Execution terminated'
          STOP
        END IF
      END IF          


C Check that the panels are not of disproportionate sizes
      IF (LVALID) THEN
        SIZMAX=ZERO
        SIZMIN=DIAM
        DO 150 ISE=1,NSE
          QA(1)=VERTEX(SELV(ISE,1),1)
          QA(2)=VERTEX(SELV(ISE,1),2)
          QB(1)=VERTEX(SELV(ISE,2),1)
          QB(2)=VERTEX(SELV(ISE,2),2)
          SIZMAX=MAX(SIZMAX,DIST2(QA,QB))
          SIZMIN=MIN(SIZMIN,DIST2(QA,QB))
150     CONTINUE
        IF (SIZMAX.GT.10.0D0*SIZMIN) THEN
          WRITE(*,*) 'WARNING(LIBEM2) - panels of disproportionate'
          WRITE(*,*) ' sizes'
        END IF
      END IF          
          
C Check that the boundary does not contain sharp angles
      IF (LVALID) THEN
        LERROR=.FALSE.
        DO 160 ISE=1,NSE
          QA(1)=VERTEX(SELV(ISE,1),1)
          QA(2)=VERTEX(SELV(ISE,1),2)
          QB(1)=VERTEX(SELV(ISE,2),1)
          QB(2)=VERTEX(SELV(ISE,2),2)
          IF (ISE.GE.2) THEN
            QC(1)=VERTEX(SELV(ISE-1,1),1)
            QC(2)=VERTEX(SELV(ISE-1,1),2)
          ELSE
            QC(1)=VERTEX(SELV(NSE,1),1)
            QC(2)=VERTEX(SELV(NSE,1),2)
          END IF
          IF (DIST2(QC,QB).LT.MAX(DIST2(QA,QB),DIST2(QA,QC))/5) THEN
            WRITE(*,*) 'Sharp angle at node ',SELV(ISE,1)
            LERROR=.TRUE.
          END IF
160     CONTINUE
        IF (LERROR) THEN
          WRITE(*,*)
          WRITE(*,*) 'WARNING(LIBEM2) - Boundary has sharp angles'
        END IF
      END IF          
     

C Validation of the surface functions
      IF (LVALID.AND.LSOL) THEN
        LERROR=.FALSE.
        DO 170 ISE=1,NSE
          IF (MAX(ABS(SALPHA(ISE)),ABS(SBETA(ISE))).LT.1.0D-6) 
     *     LERROR=.TRUE.
170     CONTINUE
        IF (LERROR) THEN
          WRITE(*,*) 
          WRITE(*,*) 'ERROR(LIBEM2) - at most one of SALPHA(i),SBETA(i)'
          WRITE(*,*) ' may be zero for all i'
          WRITE(*,*) 'Execution terminated'
          STOP
        END IF
      END IF

C Set up the quadrature rule.
C   Set up 8 point Gauss-Legendre rules
      CALL GL8(MAXNQ,NQ,WQ,AQ)

C Validation that the surface is closed
      IF (LVALID) THEN
        PA(1)=VERTEX(SELV(1,1),1)
        PA(2)=VERTEX(SELV(1,1),2)
        PB(1)=VERTEX(SELV(1,2),1)
        PB(2)=VERTEX(SELV(1,2),2)
        P(1)=(PA(1)+PB(1))/TWO
        P(2)=(PA(2)+PB(2))/TWO
        VECP(1)=0.0D0
        VECP(2)=1.0D0
        SUMM=0.0D0
        DO 180 JSE=1,NSE
C  Set QA and QB, the coordinates of the edges of the JSEth panel
          QA(1)=VERTEX(SELV(JSE,1),1)
          QA(2)=VERTEX(SELV(JSE,1),2)
          QB(1)=VERTEX(SELV(JSE,2),1)
          QB(2)=VERTEX(SELV(JSE,2),2)
C     Set LPONEL
          LPONEL=(JSE.EQ.1)

C     Switch off the validation parameter for L2LC
          LVAL=.FALSE.


C     Only M operators is required. Set LM true, 
C      LL, LMT,LN false. 
          LL=.FALSE.
          LM=.TRUE.
          LMT=.FALSE.
          LN=.FALSE.

C     Call L2LC.
          CALL L2LC(P,VECP,QA,QB,LPONEL,
     *     MAXNQ,NQ,AQ,WQ,
     *     LVAL,EGEOM,EQRULE,LFAIL1,
     *     LL,LM,LMT,LN,DISL,DISM,DISMT,DISN)
          
          SUMM=SUMM+DISM
180     CONTINUE
        IF (ABS(SUMM+0.5D0).GT.0.01) THEN
          WRITE(*,*) 
          WRITE(*,*) 'ERROR(LIBEM2) - in geometry'
          WRITE(*,*) ' The boundary could be oriented wrongly'
          WRITE(*,*) '  On the outer boundary arrange panels'
     *     // ' in clockwise order'
          WRITE(*,*) '  If there are inner boundaries arrange the'
     *     // ' panels in anticlockwise order'
          STOP
        END IF
        IF (ABS(SUMM+0.5D0).GT.0.01) THEN
          WRITE(*,*) 
          WRITE(*,*) 'WARNING(LIBEM2) - in geometry'
          WRITE(*,*) ' The boundary panels may be arranged incorrectly'
          STOP
        END IF
      END IF  


C Validation that the points in PINT are interior points
      IF (LVALID) THEN
        DO IPI=1,NPI
        P(1)=PINT(IPI,1)
        P(2)=PINT(IPI,2)
        VECP(1)=0.0D0
        VECP(2)=1.0D0
        SUMM=0.0D0
        DO 210 JSE=1,NSE
C  Set QA and QB, the coordinates of the edges of the JSEth panel
          QA(1)=VERTEX(SELV(JSE,1),1)
          QA(2)=VERTEX(SELV(JSE,1),2)
          QB(1)=VERTEX(SELV(JSE,2),1)
          QB(2)=VERTEX(SELV(JSE,2),2)
C     Set LPONEL
          LPONEL=.FALSE.

C     Only the M operator is required. Set LM true, 
C      LL,LMT,LN false. 
          LL=.FALSE.
          LM=.TRUE.
          LMT=.FALSE.
          LN=.FALSE.

C     Call L2LC.
          CALL L2LC(P,VECP,QA,QB,LPONEL,
     *     MAXNQ,NQ,AQ,WQ,
     *     LVAL,EGEOM,EQRULE,LFAIL1,
     *     LL,LM,LMT,LN,DISL,DISM,DISMT,DISN)
          
          SUMM=SUMM+DISM
210     CONTINUE

        IF (ABS(SUMM+1.0D0).GT.0.01) THEN
          WRITE(*,*) 
          WRITE(*,*) 'WARNING(LIBEM2) - The observation point'
          WRITE(*,*) ' (',P(1),',',P(2),')'
          WRITE(*,*) ' is may not be interior to the boundary'
        END IF
      END DO
      END IF



C Set up validation and control parameters
C  Switch on the validation of L2LC
      LVAL=.FALSE.
C  Set EQRULE
      EQRULE=1.0D-6


C  Compute the discrete L, M, Mt and N matrices
C   Loop(ISE) through the points on the boundary
      DO 510 ISE=1,NSE
C    Set P
        PA(1)=VERTEX(SELV(ISE,1),1)
        PA(2)=VERTEX(SELV(ISE,1),2)
        PB(1)=VERTEX(SELV(ISE,2),1)
        PB(2)=VERTEX(SELV(ISE,2),2)
        P(1)=(PA(1)+PB(1))/TWO
        P(2)=(PA(2)+PB(2))/TWO
C    Set VECP to the normal on the boundary of the panel at P
        CALL NORM2(PA,PB,VECP)
C    Loop(ISE) through the panels
        DO 520 JSE=1,NSE
C     Set QA and QB, the coordinates of the edges of the JSEth panel
          QA(1)=VERTEX(SELV(JSE,1),1)
          QA(2)=VERTEX(SELV(JSE,1),2)
          QB(1)=VERTEX(SELV(JSE,2),1)
          QB(2)=VERTEX(SELV(JSE,2),2)

C     Set LPONEL
          LPONEL=(ISE.EQ.JSE)


C     All operators are required
          LL=.TRUE.
          LM=.TRUE.
          LMT=.FALSE.
          LN=.FALSE.

C     Call L2LC.
          CALL L2LC(P,VECP,QA,QB,LPONEL,
     *     MAXNQ,NQ,AQ,WQ,
     *     LVAL,EGEOM,EQRULE,LFAIL1,
     *     LL,LM,LMT,LN,DISL,DISM,DISMT,DISN)

          L_SS(ISE,JSE)=DISL
          M_SSPLUS(ISE,JSE)=DISM
         
C    Close loop(JSE) 
520     CONTINUE

        M_SSPLUS(ISE,ISE)=M_SSPLUS(ISE,ISE)+HALF
        IF (LSOL) C(ISE)=SIPHI(ISE)
            
C   Close loop(ISE) 
510   CONTINUE

      IF (LSOL) THEN
        CALL GLS(MAXNSE,NSE,M_SSPLUS,L_SS,C,SALPHA,SBETA,SF,
     *   SPHI,SVEL,LFAIL,PERM,XORY,WKSPC1,WKSPC2,WKSPC3)
      END IF  

C  SOLUTION IN THE DOMAIN

C   Compute sound pressures at the selected interior points.
C    Loop through the the points in the interior region
      DO 800 IPI=1,NPI
C    Set P
        P(1)=PINT(IPI,1)
        P(2)=PINT(IPI,2)
C    Set VECP, this is arbitrary as the derivative of the potential
C     at P is not sought.
        VECP(1)=ONE
        VECP(2)=ZERO

C    Initialise SUMPHI to the incident potential
        SUMPHI=PDIPHI(IPI)

C    Loop(ISE) through the panels
        DO 850 JSE=1,NSE
C     Compute the discrete L and M integral operators. 
            
C     Set QA and QB, the coordinates of the edges of the JSEth panel
          QA(1)=VERTEX(SELV(JSE,1),1)
          QA(2)=VERTEX(SELV(JSE,1),2)
          QB(1)=VERTEX(SELV(JSE,2),1)
          QB(2)=VERTEX(SELV(JSE,2),2)

C     All the points do not lie on the boundary hence LPONEL=.FALSE.
          LPONEL=.FALSE.              

C     Only L, M operators are required. Set LL,LM true, 
C      LMT,LN false. 
          LL=.TRUE.
          LM=.TRUE.
          LMT=.FALSE.
          LN=.FALSE.
                
C     Call L2LC.
          CALL L2LC(P,VECP,QA,QB,LPONEL,
     *     MAXNQ,NQ,AQ,WQ,
     *     LVAL,EGEOM,EQRULE,LFAIL,
     *     LL,LM,LMT,LN,DISL,DISM,DISMT,DISN)

          IF (.NOT.LSOL) THEN
            L_DS(IPI,JSE)=DISL
            M_DS(IPI,JSE)=DISM
          END IF

C     Accumulate phi 
          IF (LSOL) SUMPHI=SUMPHI+DISL*SVEL(JSE)-DISM*SPHI(JSE)

C      Close loop (JSE) through the panels
850     CONTINUE

        PDPHI(IPI)=SUMPHI

C     Close loop(IPI) through the interior points
800   CONTINUE

      END



C Subroutines required for L2LC (not in file L2LC.FOR) 
C  Subroutine for returning the square root.
       REAL*8 FUNCTION FNSQRT(X)
       REAL*8 X
       FNSQRT=SQRT(X)
       END

C  Subroutine for returning the natural logarithm.
       REAL*8 FUNCTION FNLOG(X)
       REAL*8 X
       FNLOG=LOG(X)
       END

C ----------------------------------------------------------------------