C***************************************************************
C           Subroutine AMBEMA by Stephen Kirkup                     
C***************************************************************
C 
C  Copyright 1998- Stephen Kirkup
C  John Tyndall Nuclear Research Institute
C  School of Computing Engineering and Physical Sciences
C  University of Central Lancashire - www.uclan.ac.uk 
C  Westlakes Campus
C  Samuel Lindow Building
C  West Lakes Science and Technology Park
C  Whitehaven
C  Cumbria CA24 3JY
C  United Kingdom
C  smkirkup@uclan.ac.uk
C
C  This open source code can be found at
C   www.boundary-element-method.com/fortran/AMBEMA.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 This subroutine computes the modal solutions to the axisymmetric
C Helmholtz equation
C                  __ 2                2
C                  \/    {\phi}   +   k  {\phi}   =  0   
C
C in the domain interior to a closed axisymmetric boundary.
C
C The boundary (S) is defined (approximated) by a set of straight line 
C elements. The domain of the equation is within the boundary.
C
C The nature of the boundary condition may be Dirichlet, Robin or 
C Neumann and it is specified by the boundary functions {\alpha} and
C {\beta}. It is assumed to have the following general form
C
C         {\alpha}(q) {\phi}(q) + {\beta}(q) v(q) = 0
C    
C where {\alpha} and {\beta} are complex-valued functions defined on S.
C Important examples are {\alpha}=1, {\beta}=0 which is equivalent to a
C Dirichlet  boundary condition and {\alpha}=0, {\beta}=1 which is 
C equivalent to a Neumann boundary condition. 
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     |------->-----|      AMBEMA            |   :
C      |(e.g. ambema_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) :---<---:   | H3ALC |  :   :  | CGLS |    : 
C             :              :       :   --------   :   :  --------    :  
C             :..............:       : -------------:   : -------------:  
C                                    : |subordinate|:   : |subordinate|: 
C             ................       : | routines  |:   : | routines  |:  
C             :              :       : -------------:   : -------------: 
C             : (geom3d.for) :---<---:              :   :              : 
C             :              :       : (h3alc.for)  :   : (cgls.for)   :
C             :..............:       :..............:   :..............:
C                                    
C
C The contents of the main program must be linked to AMBEMA.FOR, 
C H3ALC.FOR, CGLS.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  elements. The elements 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 element by indicating the labels for the two edge nodes on the
C  generator and hence enumerates the elements.
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 elements. The boundary functions {\alpha} 
C  (SALPHA), {\beta} (SBETA) and f (SF) are also defined by their values
C  at the centres of the elements.
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 PIPHI.
C
C

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 element must be no more than 10x 
C  the length of the generator of the smallest element.

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 PIPHI 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 H3ALC: Returns the individual discrete Helmholtz integral
C  operators. (in file H3ALC.FOR)
C subroutine INTEIG: Solves a general linear system of equations. 
C  (in file INTEIG.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 FNEXP(Z): real X : Returns the complex exponential of Z.


C The subroutine

      SUBROUTINE AMBEMA(KA,KB,MAXNK,NK,
     *                 MAXNEIG,
     *                 MAXNV,NV,VERTEX,MAXNSE,NSE,SELV,
     *                 MAXNPI,NPI,PINT,
     *                 SALPHA,SBETA,
     *                 LVALID,EGEOM,
     *                 NEIG,EIGVAL,SPHI,SVEL,PIPHI,
     *                 WKSPC1,WKSPC2,WKSPC3,WKSPC4,WKSPC5,WKSPC6,WKSPC7,
     *                 WKSP00,WKSP01,WKSP02,WKSP03,WKSP04,WKSP05,
     *                 WKSP06,WKSP07,WKSP08,WKSP09,WKSP10,WKSP11,
     *                 WKSP12)

      PARAMETER (MAXNGQ=32)
      PARAMETER (MAXNTQ=10000)

C  Wavenumber interpolation information
C   Lower limit of k-range
      REAL*8     KA
C   Upper limit of k-range
      REAL*8     KB
C   Limit on the number of interpolation points in the k-range
      INTEGER    MAXNK
C   Number of interpolation points in the k-range
      INTEGER    NK

C  Eigenvalue information
C   Limit on the number of eigenvalues determined in the k-range
      INTEGER    MAXNEIG

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 elements describing S
      INTEGER    MAXNSE
C   The number of elements describing S
      INTEGER    NSE
C   The indices of the vertices describing each element
      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 = 0
C  where alpha, beta are complex valued functions over S.
C  The functions are set values at the collocation points.
C   function alpha
      COMPLEX*16 SALPHA(MAXNSE)
C   function beta
      COMPLEX*16 SBETA(MAXNSE)

C  Validation and control parameters
      LOGICAL    LVALID
      REAL*8     EGEOM

C  Solution 
C   Number of eigenvalues determined in the k-range
      INTEGER    NEIG
C   eigenvalues
      COMPLEX*16     EIGVAL(MAXNEIG)
C   function phi
      COMPLEX*16 SPHI(MAXNEIG,MAXNSE)
C   function vel
      COMPLEX*16 SVEL(MAXNEIG,MAXNSE)
C   domain solution
      COMPLEX*16 PIPHI(MAXNEIG,MAXNPI)

C  Working space 
      COMPLEX*16 WKSPC1(MAXNK,MAXNSE,MAXNSE)
      COMPLEX*16 WKSPC2(MAXNK,MAXNSE,MAXNSE)
      COMPLEX*16 WKSPC3(MAXNK,MAXNSE,MAXNSE)
      LOGICAL    WKSPC4(MAXNSE)

C  External function
      REAL*8     DIST2

      REAL*8     K
      COMPLEX*16 CK

C  Constants
C   Real scalars: 0, 1, 2, half, pi
      REAL*8 ZERO,ONE,TWO,HALF,PI
C   Complex scalars: (0,0), (1,0), (0,1)
      COMPLEX*16 CZERO,CONE,CIMAG

C  Geometrical description of the boundary
C   Elements counter
      INTEGER    ISE,JSE
C   The points interior to the boundary where the solution is sought 
      INTEGER    IIP
C   Parameters for H2LC
      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 element (LPONEL=.TRUE.) and
C    one for the case when P does not lie on the element.]
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 H3ALC
C   Counter through the quadrature points
      INTEGER    IGQ
C    Abscissae of the actual quadrature rule in the theta direction
      REAL*8     AGQ(MAXNGQ)
C    Weights of the actual quadrature rule
      REAL*8     WGQ(MAXNGQ)
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 H2LC
      LOGICAL    LVAL
      REAL*8     EK
      REAL*8     EQRULE
      LOGICAL    LLK
      LOGICAL    LMK
      LOGICAL    LMKT
      LOGICAL    LNK

C  Parameters for subroutine H2LC. 
      COMPLEX*16 DISLK
      COMPLEX*16 DISMK
      COMPLEX*16 DISMKT
      COMPLEX*16 DISNK

      REAL*8     WKSPCE(2*MAXNTQ+MAXNGQ)

C  Parameters for subroutine INTEIG
C   The k-coordiate of the interpolation matrices
      REAL*8 WKSPC5(MAXNK)

C The list of eigenvectors. EIGVCT(i,j) gives the value of the 
C  j-th component of the i-th eigenvector. The eigenvectors are
C  normalised
      COMPLEX*16 WKSPC6((MAXNK-1)*MAXNSE,MAXNSE)
      LOGICAL    WKSPC7(MAXNSE)

C Work space
C ----------
      COMPLEX*16 WKSP00((MAXNK-1)*MAXNSE,(MAXNK-1)*MAXNSE)
      COMPLEX*16 WKSP01((MAXNK-1)*MAXNSE,(MAXNK-1)*MAXNSE)
      COMPLEX*16 WKSP02((MAXNK-1)*MAXNSE,(MAXNK-1)*MAXNSE)
      REAL*8     WKSP03((MAXNK-1)*MAXNSE,(MAXNK-1)*MAXNSE)
      REAL*8     WKSP04((MAXNK-1)*MAXNSE,(MAXNK-1)*MAXNSE)
      REAL*8     WKSP05((MAXNK-1)*MAXNSE,(MAXNK-1)*MAXNSE)
      REAL*8     WKSP06((MAXNK-1)*MAXNSE,(MAXNK-1)*MAXNSE)
      REAL*8     WKSP07((MAXNK-1)*MAXNSE)
      REAL*8     WKSP08((MAXNK-1)*MAXNSE)
      REAL*8     WKSP09(MAXNK)
      INTEGER    WKSP10((MAXNK-1)*MAXNSE)
      COMPLEX*16 WKSP11((MAXNK-1)*MAXNSE,(MAXNK-1)*MAXNSE)
      COMPLEX*16 WKSP12((MAXNK-1)*MAXNSE)

C  Other variables
C   Failure flag
      LOGICAL    LFAIL
C   Accumulation of solution {\phi}
      COMPLEX*16 SUMPHI
C   The `diameter' of the boundary or the maximum distance between any
C    two vertices
      REAL*8     DIAM
     
      REAL*8     RADMID,SGLEN,GLEN,CIRMID,TDIV
      INTEGER    NDIV
      REAL*8     SUMMK


C Code omitted.
C In order to purchase the full code please see the author's website
C www.boundary-element-method.com
C Code forms part of the ABEMFULL package.