C***************************************************************           
C              Subroutine LIBEMA by Stephen Kirkup                      
C*************************************************************** 
C  Version 2 (July 2015) (previous version 2001) Stephen Kirkup
C  New GLS routine implemented and simplified integral equation
C   formulation
C  School of Engineering, University of Central Lancashire
C  www.uclan.ac.uk 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/LIBEMA.FOR 
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
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.
C
C                                   ....................................
C                                   :                                  :
C                                   :                                  :
C      ----------------------       :     --------------------------   :
C      |                    |       :     |                        |   :
C      |   MAIN PROGRAM     |------->-----|      LIBEMA            |   :
C      | (e.g. libema_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) :---<---:   | L3ALC |  :   :  | GLS |     : 
C             :              :       :   --------   :   :  --------    :  
C             :..............:       : -------------:   : -------------:  
C                                    : |subordinate|:   : |subordinate|: 
C             ................       : | routines  |:   : | routines  |:  
C             :              :       : -------------:   : -------------: 
C             : (geom3d.for) :---<---:              :   :              : 
C             :              :       : (l3alc.for)  :   : (gls.for)    :
C             :..............:       :..............:   :..............:
C                                    
C
C The contents of the main program must be linked to LIBEMA.FOR, 
C L3ALC.FOR, GLS.FOR, GEOM2D.FOR and GEOM3D.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 generator (the (r,z)
C  coordinate). The data structure VERTEX lists and enumerates the 
C  (r,z) coordinates of the vertices, the data structure SELV defines 
C  each panel by indicating the labels for the two edge nodes on the
C  generator 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 generator of the panels. The boundary functions {\alpha} 
C  (SALPHA), {\beta} (SBETA) and f (SF) are also defined by their values
C  at the centres of the panels.
C Normally a solution in the domain is required. By listing the (r,z)
C  coordinates of all the interior points in PINT, the subroutine
C  returns the value of {\phi} at these points in PDPHI.`
C
C
C Format and parameters of the subroutine
C ---------------------------------------
C
C The subroutine has the form
C
C      SUBROUTINE LIBEMA(MAXNV,NV,VERTEX,MAXNSE,NSE,SELV,
C     *                 MAXNPI,NPI,PINT,
C     *                 SALPHA,SBETA,SF,
C     *                 LSOL,LVALID,EGEOM,
C     *                 SPHI,SVEL,PDPHI,
C     *                 L_SS,M_SSPLUS,L_DS,M_DS,
C     *                 PERM,XORY,C,WKSPC1,M_SSPLUS,L_DS)

C The parameters to the subroutine
C ================================

C Geometry of the boundary S (input)
C integer MAXNV: The limit on the number of vertices of the generator
C  that defines (approximates) S. MAXNV>=3.
C integer NV: The number of vertices on S. 3<=NV<=MAXNV.
C real VERTEX(MAXNV,2): The coordinates of the vertices of the 
C  generator. VERTEX(i,1),VERTEX(i,2) are the r,z coordinates of the
C  i-th vertex.
C integer MAXNSE: The limit on the number of panels describing S.
C  MAXNSE>=3.
C integer NSE: The number of panels describing S. 3<=NSE<=MAXNSE.
C integer SELV(MAXNSE,2): The indices of the two vertices defining 
C  each panel. The generator of the i-th panel has vertices with
C  r,z coordinates (VERTEX(SELV(i,1),1),VERTEX(SELV(i,1),2)) and
C  (VERTEX(SELV(i,2),1),VERTEX(SELV(i,2),2)).

C Interior points at which the solution is to be observed (input)
C integer MAXNPI: Limit on the number of points interior to the
C  boundary. MAXNPI>=1.
C integer NPI: The number of interior points. 0<=NPI<=MAXNPI.
C real PINT(MAXNPI,2). The coordinates of the interior point.
C  PINT(i,1),PINT(i,2) are the r,z coordinates of the i-th point.

C The boundary condition ("{\alpha} {\phi} + {\beta} v = f") (input)
C real SALPHA(MAXNSE): The values of "{\alpha}" at the centres
C  of the generator of the panels.
C real SBETA(MAXNSE): The values of "{\beta}" at the centres
C  of the generator of the panels.
C real SF(MAXNSE): The values of "f" at the centres of the 
C  of the generator of the panels.

C Incident field
C real SFFPHI(MAXNSE): The incident potential at the centres of the 
C  generator of the panels
C real SFFVEL(MAXNSE): The derivative of the incident centres of the
C  generator of the panels
C real SIPHI(MAXNPI): The incident potential at the chosen interior 
C  points

C Validation and control parameters (input) 
C logical LSOL: A switch to control whether a particular solution
C  is required.
C logical LVALID: A switch to enable the choice of checking of 
C  subroutine parameters.
C real EGEOM: The maximum absolute error in the parameters that
C  describe the geometry.

C  Solution (output)
C  real SPHI(MAXNSE): The potential ("{\phi}") at the centres of the 
C   boundary panels.
C  real SVEL(MAXNSE):  ("v" or "d{\phi}/dn" where n is the outward 
C   normal to the boundary) at the centres of the boundary panels.
C  real PDPHI(MAXNPI): The potential ("{\phi}") at the interior points.

C Further information from system solution GLS
C      INTEGER*4  PERM(MAXNSE)
C      LOGICAL    XORY(MAXNSE)
C      REAL*8     C(MAXNSE)

C Work Space
C      REAL*8     WKSPC1(MAXNSE)
C      REAL*8     M_SSPLUS(MAXNSE)
C      REAL*8     L_DS(MAXNSE)

 
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 VERTEX(SELV(1,1),1)=0
C  and VERTEX(SELV(NSE,2),1)=0.
C (3) The indices of the nodes listed in SELV must be such that they
C  are ordered counter-clockwise around the generator of the boundary.
C (4) The generator of the largest panel must be no more than 10x 
C  the length of the generator of the 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 L3ALC: Returns the individual discrete Laplace integral
C  operators. (in file L3ALC.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 LIBEMA(MAXNV,NV,VERTEX,MAXNSE,NSE,SELV,
     *                 MAXNPI,NPI,PINT,
     *                 SALPHA,SBETA,SF,SIPHI,SIVEL,PDIPHI,
     *                 LSOL,LVALID,EGEOM,
     *                 SPHI,SVEL,PDPHI,
     *                 L_SS,M_SSPLUS,L_DS,M_DS,
     *                 PERM,XORY,C,ACOLNORMS,BCOLNORMS,WKSPC)
      PARAMETER (MAXNGQ=32)
      PARAMETER (MAXNTQ=10000)

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 derivative of the incident potential on S
      REAL*8 SIVEL(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  Working space 
      REAL*8 L_SS(MAXNSE,MAXNSE)
      REAL*8 M_SSPLUS(MAXNSE,MAXNSE)
      REAL*8 L_DS(MAXNPI,MAXNSE)
      REAL*8 M_DS(MAXNPI,MAXNSE)

C Further information from system solution GLS
      INTEGER*4  PERM(MAXNSE)
      LOGICAL    XORY(MAXNSE)
      REAL*8     C(MAXNSE)

C Work Space
      REAL*8     ACOLNORMS(MAXNSE)
      REAL*8     BCOLNORMS(MAXNSE)
      REAL*8     WKSPC(MAXNSE)


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 L3ALC
      REAL*8     P(2),PA(2),PB(2),QA(2),QB(2),VECP(2)
      LOGICAL    LPONEL


C  Quadrature rule information
C   [Note that in this program two quadrature rules are used: one for
C    the case when the point P lies on the panel (LPONEL=.TRUE.) and
C    one for the case when P does not lie on the panel.]
C   Quadrature rule used when LPONEL=.TRUE.
C    Number of quadrature points
      INTEGER    NGQON
C    Abscissae of the actual quadrature rule
      REAL*8     AGQON(MAXNGQ)
C    Weights of the actual quadrature rule
      REAL*8     WGQON(MAXNGQ)
C   Quadrature rule used when LPONEL=.FALSE.
C    Number of quadrature points
      INTEGER    NGQOFF
C    Abscissae of the actual quadrature rule
      REAL*8     AGQOFF(MAXNGQ)
C    Weights of the actual quadrature rule
      REAL*8     WGQOFF(MAXNGQ)
C   Quadrature rule parameters for L3ALC
C    Actual number of quadrature points
      INTEGER    NGQ
C    Abscissae of the actual quadrature rule in the generator direction
      REAL*8     AGQ(MAXNGQ)
C    Weights of the actual quadrature rule
      REAL*8     WGQ(MAXNGQ)
C   Counter through the quadrature points
      INTEGER    IGQ
C    Abscissae of the actual quadrature rule in the theta direction
      REAL*8     ATQ(MAXNTQ)
C    Weights of the actual quadrature rule
      REAL*8     WTQ(MAXNTQ)

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

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

      REAL*8     WKSPCE(2*MAXNTQ+MAXNGQ)


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
      REAL*8     RADMID,SGLEN,GLEN,CIRMID,TDIV
      INTEGER    NDIV
      REAL*8     DIST2


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

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

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


C Validation of parameters of LIBEMA
C ---------------------------------

      IF (LVALID) THEN

C Validation of main paramters
        LERROR=.FALSE.
        IF (MAXNV.LT.3) THEN
          WRITE(*,*) 'MAXNV = ',MAXNV
          WRITE(*,*) 'ERROR(LIBEMA) - must have MAXNV>=3'
          LERROR=.TRUE.
        END IF
        IF (NV.LT.3.OR.NV.GT.MAXNV) THEN
          WRITE(*,*) 'NV = ',NV
          WRITE(*,*) 'ERROR(LIBEMA) - must have 3<=NV<=MAXNV'
          LERROR=.TRUE.
        END IF
        IF (MAXNSE.LT.3) THEN
          WRITE(*,*) 'MAXNSE = ',MAXNSE
          WRITE(*,*) 'ERROR(LIBEMA) - must have MAXNSE>=3'
          LERROR=.TRUE.
        END IF
        IF (NSE.LT.3.OR.NSE.GT.MAXNSE) THEN
          WRITE(*,*) 'NSE = ',NSE
          WRITE(*,*) 'ERROR(LIBEMA) - must have 3<=NSE<=MAXNSE'
          LERROR=.TRUE.
        END IF
        IF (MAXNPI.LT.1) THEN
          WRITE(*,*) 'MAXNPI = ',MAXNPI
          WRITE(*,*) 'ERROR(LIBEMA) - must have MAXNPI>=1'
          LERROR=.TRUE.
        END IF
        IF (NPI.LT.0.OR.NPI.GT.MAXNPI) THEN
          WRITE(*,*) 'NPI = ',NPI
          WRITE(*,*) 'ERROR(LIBEMA) - must have 3<=NPI<=MAXNPI'
          LERROR=.TRUE.
        END IF
        IF (EGEOM.LE.ZERO) THEN
          WRITE(*,*) 'NPI = ',NPI
          WRITE(*,*) 'ERROR(LIBEMA) - EGEOM must be positive'
          LERROR=.TRUE.
        END IF
        IF (LERROR) THEN
          LFAIL=.TRUE.
          WRITE(*,*)
          WRITE(*,*) 'Error(s) found in the main parameters of LIBEMA'
          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 (VERTEX(SELV(1,1),1).GT.EGEOM) LERROR=.TRUE.
        IF (VERTEX(SELV(NSE,2),1).GT.EGEOM) LERROR=.TRUE.
        IF (LERROR) THEN
          WRITE(*,*) 'ERROR(LIBEMA) - 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(LIBEMA) - 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(LIBEMA) - 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(LIBEMA) - 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=2,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)
          QC(1)=VERTEX(SELV(ISE-1,1),1)
          QC(2)=VERTEX(SELV(ISE-1,1),2)
          IF (DIST2(QC,QB).LT.MAX(DIST2(QA,QB),DIST2(QA,QC))) THEN
            WRITE(*,*) 'Sharp angle at node ',SELV(ISE,1)
            LERROR=.TRUE.
          END IF
160     CONTINUE
        IF (LERROR) THEN
          WRITE(*,*)
          WRITE(*,*) 'WARNING(LIBEMA) - 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(LIBEMA) - 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 validation and control parameters
C  Switch off the validation of L3ALC
      LVAL=.FALSE.
C  Set EQRULE
      EQRULE=1.0D-6

C Set up the quadrature rule(s).
C  Set up quadrature rule for the case when P is not on the panel.
C   Set up 8 point Gauss-Legendre rules
      CALL GL8(MAXNGQ,NGQOFF,WGQOFF,AGQOFF)
C  Set up quadrature rule for the case when P is on the panel.
C   This is done by splitting the quadrature rule at the centre.
      NGQON=2*NGQOFF
      DO 330 IGQ=1,NGQOFF
        AGQON(IGQ)=AGQOFF(IGQ)/TWO
        AGQON(NGQOFF+IGQ)=0.5D0+AGQOFF(IGQ)/TWO
        WGQON(IGQ)=WGQOFF(IGQ)/TWO
        WGQON(NGQOFF+IGQ)=WGQOFF(IGQ)/TWO
330   CONTINUE



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 Quadrature rule in the theta direction is constructed out of individual
C Gauss rules so that the length of each is approximately equal to the
C length of the panel at the generator.
          NGQ=NGQON
          RADMID=(QA(1)+QB(1))/TWO
          SGLEN=(QA(1)-QB(1))*(QA(1)-QB(1))+
     *     (QA(2)-QB(2))*(QA(2)-QB(2))
          GLEN=SQRT(SGLEN)
          CIRMID=PI*RADMID
          NDIV=1+CIRMID/GLEN
          TDIV=ONE/DBLE(NDIV)
          NTQ=NDIV*NGQ
          IF (NTQ.GT.MAXNTQ) THEN
            WRITE(*,*) 'ERROR(LIBEMA) - MAXNTQ is set too small'
            STOP
          END IF
          DO 146 IDIV=1,NDIV
            DO 156 IGQ=1,NGQ
              WTQ((IDIV-1)*NGQ+IGQ)=WGQON(IGQ)/DBLE(NDIV)
              ATQ((IDIV-1)*NGQ+IGQ)=AGQON(IGQ)/DBLE(NDIV)+
     *         TDIV*DBLE(IDIV-1)
156         CONTINUE
146       CONTINUE



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

C     Call L3ALC.
          CALL L3ALC(P,VECP,QA,QB,LPONEL,
     *     MAXNGQ,NGQON,AGQON,WGQON,MAXNTQ,NTQ,ATQ,WTQ,
     *     LVAL,EGEOM,EQRULE,LFAIL,
     *     LL,LM,LMT,LN,DISL,DISM,DISMT,DISN,
     *     WKSPCE)
          
          SUMM=SUMM+DISM
180     CONTINUE
        IF (ABS(SUMM-0.5D0).LT.0.01) THEN
          WRITE(*,*) 
          WRITE(*,*) 'ERROR(LIBEMA) - 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(LIBEMA) - in geometry'
          WRITE(*,*) ' The boundary panels may be arranged incorrectly'
        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 Quadrature rule in the theta direction is constructed out of individual
C Gauss rules so that the length of each is approximately equal to the
C length of the panel at the generator.
          NGQ=NGQON
          RADMID=(QA(1)+QB(1))/TWO
          SGLEN=(QA(1)-QB(1))*(QA(1)-QB(1))+
     *     (QA(2)-QB(2))*(QA(2)-QB(2))
          GLEN=SQRT(SGLEN)
          CIRMID=PI*RADMID
          NDIV=1+CIRMID/GLEN
          TDIV=ONE/DBLE(NDIV)
          NTQ=NDIV*NGQ
          IF (NTQ.GT.MAXNTQ) THEN
            WRITE(*,*) 'ERROR(LIBEMA) - MAXNTQ is set too small'
            STOP
          END IF
          DO 147 IDIV=1,NDIV
            DO 157 IGQ=1,NGQ
              WTQ((IDIV-1)*NGQ+IGQ)=WGQON(IGQ)/DBLE(NDIV)
              ATQ((IDIV-1)*NGQ+IGQ)=AGQON(IGQ)/DBLE(NDIV)+
     *         TDIV*DBLE(IDIV-1)
157         CONTINUE
147       CONTINUE

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 L3ALC.
          CALL L3ALC(P,VECP,QA,QB,LPONEL,
     *     MAXNGQ,NGQON,AGQON,WGQON,MAXNTQ,NTQ,ATQ,WTQ,
     *     LVAL,EGEOM,EQRULE,LFAIL,
     *     LL,LM,LMT,LN,DISL,DISM,DISMT,DISN,
     *     WKSPCE)
          
          SUMM=SUMM+DISM
210     CONTINUE

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


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
          IF (ISE.EQ.JSE) THEN
            LPONEL=.TRUE.
          ELSE
            LPONEL=.FALSE.
          END IF

C     Select quadrature rule for L3ALC
C   :  Select the quadrature rule AGQON-WGQON in the case when the
C   :   point p lies on the panel, otherwise select AGQOFF-WGQOFF
C      [Note that the overall method would benefit from selecting from
C       a wider set of quadrature rules, and an appropriate method
C       of selection]
          IF (LPONEL) THEN
          NGQ=NGQON
          DO 600 IGQ=1,NGQ
            AGQ(IGQ)=AGQON(IGQ)
            WGQ(IGQ)=WGQON(IGQ)
600       CONTINUE
          ELSE 
          NGQ=NGQOFF
          DO 610 IGQ=1,NGQ
            AGQ(IGQ)=AGQOFF(IGQ)
            WGQ(IGQ)=WGQOFF(IGQ)
610       CONTINUE
          END IF



C Quadrature rule in the theta direction is constructed out of individual
C Gauss rules so that the length of each is approximately equal to the
C length of the panel at the generator.
          RADMID=(QA(1)+QB(1))/TWO
          SGLEN=(QA(1)-QB(1))*(QA(1)-QB(1))+
     *     (QA(2)-QB(2))*(QA(2)-QB(2))
          GLEN=SQRT(SGLEN)
          CIRMID=PI*RADMID
          NDIV=1+CIRMID/GLEN
          TDIV=ONE/DBLE(NDIV)
          NTQ=NDIV*NGQ
          IF (NTQ.GT.MAXNTQ) THEN
            WRITE(*,*) 'ERROR(LIBEMA) - MAXNTQ is set too small'
            STOP
          END IF
          DO 145 IDIV=1,NDIV
            DO 155 IGQ=1,NGQ
              WTQ((IDIV-1)*NGQ+IGQ)=WGQ(IGQ)/DBLE(NDIV)
              ATQ((IDIV-1)*NGQ+IGQ)=AGQ(IGQ)/DBLE(NDIV)+
     *         TDIV*DBLE(IDIV-1)
155         CONTINUE
145       CONTINUE


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


C    Call of L3ALC routine to compute [L], [M], [Mt], [N]
          CALL L3ALC(P,VECP,QA,QB,LPONEL,
     *     MAXNGQ,NGQ,AGQ,WGQ,MAXNTQ,NTQ,ATQ,WTQ,
     *     LVAL,EGEOM,EQRULE,LFAIL,
     *     LL,LM,LMT,LN,DISL,DISM,DISMT,DISN,
     *     WKSPCE)

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

        L_SS(ISE,ISE)=L_SS(ISE,ISE)
        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, ACOLNORMS, BCOLNORMS, WKSPC)
      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 Quadrature rule in the generator direction
          NGQ=NGQOFF
          DO 1610 IGQ=1,NGQ
            AGQ(IGQ)=AGQOFF(IGQ)
            WGQ(IGQ)=WGQOFF(IGQ)
1610      CONTINUE

C Quadrature rule in the theta direction is constructed out of individual
C Gauss rules so that the length of each is approximately equal to the
C length of the panel at the generator.
          RADMID=(QA(1)+QB(1))/TWO
          SGLEN=(QA(1)-QB(1))*(QA(1)-QB(1))+
     *     (QA(2)-QB(2))*(QA(2)-QB(2))
          GLEN=SQRT(SGLEN)
          CIRMID=PI*RADMID
          NDIV=1+CIRMID/GLEN
          TDIV=ONE/DBLE(NDIV)
          NTQ=NDIV*NGQ
          IF (NTQ.GT.MAXNTQ) THEN
            WRITE(*,*) 'ERROR(LIBEMA) - MAXNTQ is set too small'
            STOP
          END IF
          DO 1145 IDIV=1,NDIV
            DO 1155 IGQ=1,NGQ
              WTQ((IDIV-1)*NGQ+IGQ)=WGQ(IGQ)/DBLE(NDIV)
              ATQ((IDIV-1)*NGQ+IGQ)=AGQ(IGQ)/DBLE(NDIV)+
     *         TDIV*DBLE(IDIV-1)
1155        CONTINUE
1145      CONTINUE
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 L3ALC.
             CALL L3ALC(P,VECP,QA,QB,LPONEL,
     *        MAXNGQ,NGQ,AGQ,WGQ,MAXNTQ,NTQ,ATQ,WTQ,
     *        LVAL,EGEOM,EQRULE,LFAIL,
     *        LL,LM,LMT,LN,DISL,DISM,DISMT,DISN,
     *        WKSPCE)

          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 ----------------------------------------------------------------------

C Subordinate routines for LIBEMA
C ===============================

C ----------------------------------------------------------------------
C            Subroutine GL8.FOR by www.numerical-methods.com           |
C ----------------------------------------------------------------------
C
C Subroutine GL8 assigns the weights and points of a 8 point Gaussian
C (Gauss-Legendre) quadrature rule defined on the interval [0,1].
C
C SUBROUTINE GL8(MAXN, N, WTS, PTS)
C integer  maxn: the maximimum number of weights/points
C integer     n: the number of weights/pointsr
C real      wts: the weights
C real      pts: the points
C
C Source of the code: http://www.numerical-methods.com/fortran/GL8.FOR
C Source of the user-guide: http://www.numerical-methods.com/fortran/
C  gl8.htm
C
C Licence: This is 'open source'; the software may be used and applied
C  within other systems as long as its provenance is appropriately
C  acknowledged. See the GNU Licence http://www.gnu.org/licenses/lgpl.txt
C  for more information or contact webmaster@numerical-methods.com

C Original code 1998. Documentation enhanced 2014

       SUBROUTINE GL8(MAXN,N,WTS,PTS)
       INTEGER MAXN
       INTEGER N
       REAL*8 WTS(MAXN)
       REAL*8 PTS(MAXN)
       N=8
       WTS(1)=       5.061426814519E-02
       PTS(1)=       1.985507175123E-02
       WTS(2)=         0.111190517227
       PTS(2)=         0.101666761293
       WTS(3)=         0.156853322939
       PTS(3)=         0.237233795042
       WTS(4)=         0.181341891689
       PTS(4)=         0.408282678752
       WTS(5)=         0.181341891689
       PTS(5)=         0.591717321248
       WTS(6)=         0.156853322939
       PTS(6)=         0.762766204958
       WTS(7)=         0.111190517227
       PTS(7)=         0.898333238707
       WTS(8)=       5.061426814519E-02
       PTS(8)=         0.980144928249
       END


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


