      PROGRAM FIELD3D

C     Generate file of replica vectors

C     Multiple receivers depths are handled inefficiently.

      INTEGER PRTFIL
      PARAMETER ( PRTFIL = 6,
     &   MXNR  = 5001, MXNRT  = 102000,  MXNT = 201,  MXNRD = 51,
     &   MXELT = 20000, MXNODE = 9000,   MAXM = 200, MXSET = 9000,
     &   MXNSD = 5 )

C     Some of the above parameters also occur in readin below
C     if you change MAXM, change it in eval3d.f, evalpdq.f too!

      INTEGER   NODE( 3, MXELT ), ADJELT( 3, MXELT ), M( MXSET ),
     &          ISET( MXNODE )
      REAL      R( MXNR ), THETA( MXNT ), SD( MXNSD ), RD( MXNRD ), 
     &          X( MXNODE ), Y( MXNODE )
      COMPLEX   CK( MAXM, MXSET ), P( MXNRT ), PHIS( MAXM, MXNSD, 3),
     &          PHIR( MAXM, MXSET ), PHIST( MAXM, 3 )
      CHARACTER OPT*50, TITLE*80, FILNAM( MXNODE )*50, TITENV*80

C     *** Read in all the user input ***

      CALL READIN( TITLE, OPT, MLIMIT, XS, YS, SD, NSD, RD, NRD,
     &   RMIN, RMAX, DELTAR, NR, R, THETA, NTHETA,
     &   X, Y, FILNAM, NNODES, NODE, NELTS )

C     *** Build ADJELT and identify the source element ***

      CALL GETADJ( NODE, NELTS, ADJELT )
      IELTSC = ISCELT( XS, YS, NODE, NELTS, X, Y )

C     Read modes at the source depths

      MSRC = MAXM
      DO ICOR = 1, 3
         INODE = NODE( ICOR, IELTSC )

         CALL GETMOD( 0, FILNAM( INODE ), MAXM, SD, NSD, 'N',
     &      CK( 1, ICOR ), PHIS( 1, 1, ICOR ), M( ICOR ), FREQ, TITENV )
         MSRC = MIN( MSRC, M( ICOR ) )
      END DO

C     *** Write header ***

      IF ( NRD .EQ. 1 ) THEN
         CALL HEADER( TITLE, SD, NSD, THETA, NTHETA,
     &      R, NR, FREQ, ATTEN, 'POLAR', XS, YS, RD(1)    )
      ELSE IF ( NSD .EQ. 1 ) THEN
         CALL HEADER( TITLE, RD, NRD, THETA, NTHETA,
     &      R, NR, FREQ, ATTEN, 'POLAR', XS, YS, SD(1)    )
      ELSE IF ( NTHETA .EQ. 1 ) THEN
         CALL HEADER( TITLE, SD, NSD, RD,    NRD,
     &      R, NR, FREQ, ATTEN, 'CARTE', XS, YS, THETA(1) )
      ENDIF

C     ****** MAIN Loop: over receiver depth ******

      DO IRD = 1, NRD
         WRITE( *, * ) 'IRD, RD( IRD ) = ', IRD, RD( IRD )

C        --- Get modes at the receiver depth
         CALL SLURP( NNODES, FILNAM, MXSET, MAXM, RD( IRD ),
     &      ISET, NSETS, M, CK, PHIR )

C        *** Loop over source depth ***

         DO ISD = 1, NSD
            WRITE( *, * ) 'ISD, SD( ISD ) = ', ISD, SD( ISD )

C           --- Get modes at the source depth
            DO ICOR = 1, 3
               DO MODE = 1, MSRC
                  PHIST( MODE, ICOR ) = PHIS( MODE, ISD, ICOR )
               END DO
            END DO

C           Call the appropriate routine to evaluate the field

            IF ( OPT(1:3) .EQ. 'STD' ) THEN
C              --- STANDARD (ignores hor. refraction)
               CALL EVAL3D( CK, PHIR, PHIST, M, IELTSC, ADJELT, 
     &            ISET, X, Y, NNODES, NODE, NELTS, XS, YS, 
     &            THETA, NTHETA, RMIN, RMAX, NR, MLIMIT, OPT, P )
            ELSE IF ( OPT(1:3) .EQ. 'PDQ' ) THEN
C              --- Quick (ignores and ignores)
               CALL EVALPDQ( CK, PHIR, PHIST, M, IELTSC, ADJELT, 
     &            ISET, X, Y, NNODES, NODE, NELTS, XS, YS, 
     &            THETA, NTHETA, RMIN, RMAX, NR, MLIMIT, OPT, P )
            ELSE IF ( OPT(1:3) .EQ. 'GBT' ) THEN
C              --- GAUSSIAN BEAM (include hor. refraction)
               CALL EVALGB( CK, PHIR, PHIST, M, MAXM, IELTSC, ADJELT, 
     &            ISET, X, Y, NNODES, NODE, NELTS, XS, YS, 
     &            THETA, NTHETA, RMIN, RMAX, NR, MLIMIT, OPT, P )
            ELSE
               CALL ERROUT( PRTFIL, 'F', 'FIELD3D', 'Unknown option' )
            ENDIF

C           Write out the field

            IF ( NRD .EQ. 1 ) THEN
               IREC = 6 + ( ISD - 1 ) * NTHETA
               CALL WRTFLD( P, NTHETA, NR, IREC )
            ELSE IF ( NSD .EQ. 1 ) THEN
               IREC = 6 + ( IRD - 1 ) * NTHETA
               CALL WRTFLD( P, NTHETA, NR, IREC )
            ELSE IF ( NTHETA .EQ. 1 ) THEN
               IREC = 6 + ( ISD - 1 ) * NRD + (IRD-1)
               CALL WRTFLD( P,  1, NR, IREC )
            ENDIF

         END DO     ! Next source   depth
      END DO     ! Next receiver depth

      STOP
      END
C**********************************************************************C
      SUBROUTINE READIN( TITLE, OPT, MLIMIT, XS, YS, SD, NSD, RD, NRD,
     &   RMIN, RMAX, DELTAR, NR, R, THETA, NTHETA, X, Y, FILNAM, NNODES,
     &   NODE, NELTS )

C     Reads in all the user input

      INTEGER PRTFIL, FLPFIL
      PARAMETER ( FLPFIL = 5, PRTFIL = 6,
     &   MXNR  = 5001, MXNRT  = 102000,  MXNT = 201,  MXNRD = 51,
     &   MXNODE = 9000,   MAXM = 200, MXNSD = 5 )

      INTEGER   NODE( 3, * )
      REAL      R( * ), THETA( * ), SD( * ), RD( * ), X( * ), Y( * )
      CHARACTER OPT*50, TITLE*80, FILNAM( * )*50

C     *** TITLE, OPT, MLIMIT ***

      READ(  FLPFIL, * ) TITLE
      WRITE( PRTFIL, * ) TITLE

      READ(  FLPFIL, * ) OPT
      WRITE( PRTFIL, * ) 'OPT = ', OPT
C      CALL STR$UPCASE( OPT, OPT )

      READ(  FLPFIL, * ) MLIMIT
      WRITE( PRTFIL, * ) 'Number of modes = ', MLIMIT
      MLIMIT = MIN( MLIMIT, MAXM )

C     *** Read source/receiver information ***

      READ(  FLPFIL, * ) XS, YS
      WRITE( PRTFIL, * ) 'Coords. of source = ', XS, YS

      XS = 1000.0 * XS
      YS = 1000.0 * YS

C     --- Read source/rcvr depths
      CALL SDRD( FLPFIL, PRTFIL, 0.0, 1.0E6,
     &   MXNSD, NSD, SD,
     &   MXNRD, NRD, RD )

C     --- Radii of receivers
      READ(  FLPFIL, * )  RMIN, RMAX, NR
      WRITE( PRTFIL, * ) 'RMIN, RMAX, NR = ', RMIN, RMAX, NR

      IF ( NR .GT. MXNR ) THEN
         WRITE( PRTFIL, * ) 'MXNR = ', MXNR
         CALL ERROUT( PRTFIL, 'F', 'FIELD3D-READIN', 'Too many ranges' )
      ENDIF

      RMIN   = 1000.0 * RMIN
      RMAX   = 1000.0 * RMAX
      DELTAR = ( RMAX - RMIN ) / ( NR - 1 )

      DO J = 1, NR
         R( J ) = RMIN + ( J - 1 ) * DELTAR
      END DO

C     --- Read angles for radials
      THETA( 3 ) = -999.9
      READ( FLPFIL, * ) NTHETA, ( THETA( ITH ), ITH = 1, NTHETA )
      CALL SUBTAB( THETA, NTHETA )

      IF ( NTHETA .GT. MXNT ) THEN
         WRITE( PRTFIL, * ) 'MXNT = ', MXNT
         CALL ERROUT( PRTFIL, 'F', 'FIELD3D-READIN', 'Too many radials')
      ENDIF

      IF ( NR * NTHETA .GT. MXNRT ) THEN
         WRITE( PRTFIL, * ) 'MXNRT = ', MXNRT
         CALL ERROUT( PRTFIL, 'F', 'FIELD3D-READIN',
     &                                      'Too many receiver points' )
      ENDIF

      IF ( NSD .NE. 1 .AND. NRD .NE. 1 .AND. NTHETA .NE. 1 ) THEN
         CALL ERROUT( PRTFIL, 'F', 'FIELD3D-READIN',
     &      'Multiple sources, rcvrs, and bearings in a single run' )
      END IF

C     *** Read nodal coordinates ***

      READ(  FLPFIL, * ) NNODES
      WRITE( PRTFIL, * ) 'NNODES = ', NNODES

      IF ( NNODES .GT. MXNODE ) THEN
         WRITE( PRTFIL, * ) 'MXNODE = ', MXNODE
         CALL ERROUT( PRTFIL, 'F', 'FIELD3D-READIN', 'Too many nodes' )
      ENDIF

      DO I = 1, NNODES

         READ( FLPFIL, * ) X( I ), Y( I ), FILNAM( I )

         X( I ) = 1000.0 * X( I )
         Y( I ) = 1000.0 * Y( I )

C         CALL STR$UPCASE( FILNAM( I ), FILNAM( I ) )

C        --- Append the extension '.MOD' to FILNAM
         JJ = INDEX( FILNAM( I ), ' ' )
         FILNAM( I )( JJ:JJ+3 ) = '.mod'

      END DO

C     *** Read in element definitions ***

      READ(  FLPFIL, * ) NELTS
      WRITE( PRTFIL, * ) 'NELTS  = ', NELTS

      DO IELT = 1, NELTS
         READ( FLPFIL, * ) ( NODE( ICOR, IELT ), ICOR=1, 3 )
      END DO

      IF ( OPT(4:4) .EQ. 'T' ) THEN
         WRITE( PRTFIL, * ) 'Performing a Tesselation check'
         CALL TESCHK( X, Y, NODE, NELTS, PRTFIL )
         WRITE( PRTFIL, * ) 'Passed the Tesselation check'
      ENDIF

      RETURN
      END
C**********************************************************************C
      SUBROUTINE GETADJ( NODE, NELTS, ADJELT )

C     Constructs a table ADJELT(IELT, ISIDE) which gives the
C     element number which shares ISIDE with element number IELT

      INTEGER NODE( 3, * ), ADJELT( 3, * ), ICOR( 3, 2 )
      DATA ((ICOR(ISIDE, J), ISIDE=1, 3), J=1, 2) /1, 2, 3, 2, 3, 1/

C     ICOR maps a side (1, 2 or 3) and a local node (1 or 2) to a
C     corner (1, 2, OR 3) of the triangle

      DO IELT = 1, NELTS
         DO ISIDE = 1, 3
            ADJELT( ISIDE, IELT ) = 0
         END DO
      END DO

C     *** Loop over each triangle ***

      DO IELT = 1, NELTS

C        *** Loop over triangle sides ***

         DO ISIDE = 1, 3

            IF ( ADJELT( ISIDE, IELT ) .EQ. 0 ) THEN

            NODE1 = NODE( ICOR( ISIDE, 1 ), IELT )
            NODE2 = NODE( ICOR( ISIDE, 2 ), IELT )

C           --- Search through other elements to find common side
            DO IELTT = 1, NELTS
               IF ( IELTT .NE. IELT ) THEN
                  DO ISIDET = 1, 3
                     NODE1T = NODE( ICOR( ISIDET, 1 ), IELTT )
                     NODE2T = NODE( ICOR( ISIDET, 2 ), IELTT )

C                    --- Do IELT and IELTT share this side?
                     IF ( ( NODE1 .EQ. NODE1T .AND.
     &                             NODE2 .EQ. NODE2T ) .OR.
     &                    ( NODE1 .EQ. NODE2T .AND.
     &                             NODE2 .EQ. NODE1T ) ) THEN
                         ADJELT( ISIDE,  IELT  ) = IELTT
                         ADJELT( ISIDET, IELTT ) = IELT
                         GOTO 5000
                     ENDIF

                  END DO  ! next side
               ENDIF
            END DO   ! next element
            ENDIF
 5000    ENDDO   ! next side

      END DO   ! next element

      RETURN
      END
C**********************************************************************C
      FUNCTION ISCELT( XS, YS, NODE, NELTS, X, Y )

C     Identifies the element containing the source at (XS, YS)

C     We define a function ENCL( XS, YS ) which is 1 when (XS , YS)
C     is inside a given triangle and decreases from 1 the further the
C     source moves away from the triangle.

C     The element with the highest value of ENCL is identified as the source
C     element.

C     If several elements enclose, the highest numbered
C     element is given posession

      INTEGER NODE( 3, * )
      REAL    X( * ), Y( * )

      ENCLMX = 0.0

      DO IELT = 1, NELTS
         NODE1 = NODE( 1, IELT )
         NODE2 = NODE( 2, IELT )
         NODE3 = NODE( 3, IELT )

         X1 = X( NODE1 )
         Y1 = Y( NODE1 )
         X2 = X( NODE2 )
         Y2 = Y( NODE2 )
         X3 = X( NODE3 )
         Y3 = Y( NODE3 )

C        --- Compute areas of triangles
         DELTA =
     &        ( X2*Y3 - Y2*X3 ) - ( X1*Y3 - Y1*X3 ) + ( X1*Y2 - Y1*X2 )

         A1 = ( X2*Y3 - Y2*X3 ) - ( XS*Y3 - YS*X3 ) + ( XS*Y2 - YS*X2 )
         A2 = ( XS*Y3 - YS*X3 ) - ( X1*Y3 - Y1*X3 ) + ( X1*YS - Y1*XS )
         A3 = ( X2*YS - Y2*XS ) - ( X1*YS - Y1*XS ) + ( X1*Y2 - Y1*X2 )

         ENCL = ABS( DELTA ) / ( ABS( A1 ) + ABS( A2 ) + ABS( A3 ) )

         IF ( ENCL .GT. ENCLMX ) THEN
            ISCELT = IELT
            ENCLMX = ENCL
         ENDIF

      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE TESCHK( X, Y, NODE, NELTS, PRTFIL )

C     Checks to see that triangulation is a tesselation
C     that is, that there are no overlapping triangles
C     (holes or triangles inside triangles are still possible)

      INTEGER NODE( 3, * ), ICOR( 3, 2 ), PRTFIL
      REAL    X( * ), Y( * )

      DATA ((ICOR(ISIDE, J), ISIDE=1, 3), J=1, 2) /1, 2, 3, 2, 3, 1/

C     ICOR maps a side (1, 2 or 3) and a local node (1 or 2) to a
C     corner (1, 2, OR 3) of the triangle

      DO IELT1 = 1, NELTS-1
         DO ISIDE1 = 1, 3
            X1 = X( NODE( ICOR( ISIDE1, 1 ), IELT1 ) )
            Y1 = Y( NODE( ICOR( ISIDE1, 1 ), IELT1 ) ) 
            X2 = X( NODE( ICOR( ISIDE1, 2 ), IELT1 ) )
            Y2 = Y( NODE( ICOR( ISIDE1, 2 ), IELT1 ) ) 
            UX = X2 - X1
            UY = Y2 - Y1

            DO IELT2 = IELT1 + 1, NELTS
               DO ISIDE2 = 1, 3
                  X3 = X( NODE( ICOR( ISIDE2, 1 ), IELT2 ) )
                  Y3 = Y( NODE( ICOR( ISIDE2, 1 ), IELT2 ) ) 
                  X4 = X( NODE( ICOR( ISIDE2, 2 ), IELT2 ) )
                  Y4 = Y( NODE( ICOR( ISIDE2, 2 ), IELT2 ) )
                  VX = X4 - X3
                  VY = Y4 - Y3
                  WX = X4 - X2
                  WY = Y4 - Y2
                  DEL = UX * VY - UY * VX

                  IF ( ABS( DEL ) .GT.
     &               MAX( ABS( UX*UY ), ABS( VX*VY ) ) / 10000.0 ) THEN
                     S1 = ( UX * WY - UY * WX ) / DEL
                     S2 = ( VX * WY - VY * WX ) / DEL
                     IF ( 0.001 .LT. S1 .AND. S1 .LT. 0.999 .AND.
     &                    0.001 .LT. S2 .AND. S2 .LT. 0.999 ) THEN
                        WRITE( PRTFIL, * )
     &                     'Sides cross for elts = ', IELT1, IELT2
                        WRITE( PRTFIL, * ) 'SIDE 1:'
                        WRITE( PRTFIL, * ) '   (', X1, Y1, ')'
                        WRITE( PRTFIL, * ) '   (', X2, Y2, ')'
                        WRITE( PRTFIL, * ) '   S = ', S1

                        WRITE( PRTFIL, * ) 'SIDE 2:'
                        WRITE( PRTFIL, * ) '   (', X3, Y3, ')'
                        WRITE( PRTFIL, * ) '   (', X4, Y4, ')'
                        WRITE( PRTFIL, * ) '   S = ', S2
                        STOP
                     ENDIF
                  ENDIF

               END DO
            END DO

         END DO
      END DO

      RETURN
      END                   
C**********************************************************************C
      SUBROUTINE SLURP( NNODES, FILNAM, MXSET, MAXM, RD,
     &   ISET, NSETS, M, CK, PHIR )

C     Reads in the values of the modes at the given receiver depth
C     for every node in the triangulation

      INTEGER PRTFIL
      PARAMETER ( PRTFIL = 6 )

      INTEGER   M( * ), ISET( * )
      COMPLEX   CK( MAXM, * ), PHIR( MAXM, * )
      CHARACTER FILNAM( * )*50, TITENV*80

      NSETS = 0

      DO 2000 INODE = 1, NNODES

C        --- Check if the modes have already been read
         IF ( INODE .GE. 2 ) THEN
            DO JNODE = 1, INODE-1
               IF ( FILNAM( INODE ) .EQ. FILNAM( JNODE ) ) THEN
C                 --- Copy previously read modes
                  ISET( INODE ) = ISET( JNODE )
                  GOTO 2000   ! process next node
               ENDIF
            END DO
         ENDIF

         NSETS = NSETS + 1
         IF ( NSETS .GT. MXSET ) THEN
         WRITE( PRTFIL, * ) 'MXSET = ', MXSET
            CALL ERROUT( PRTFIL, 'F', 'FIELD3D-SLURP',
     &                             'Too many mode sets' )
         ENDIF
         ISET( INODE ) = NSETS

C        --- Check for 'DUMMY' elts (acoustic absorbers)
         IF ( FILNAM( INODE )(1:5) .EQ. 'DUMMY' ) THEN
            M( NSETS ) = 0
         ELSE
C           --- Necessary to read in modes
            WRITE( *, * ) 'Reading INODE', INODE
            CALL GETMOD( 0, FILNAM( INODE ), MAXM, RD, 1, 'N',
     &         CK( 1, NSETS ), PHIR( 1, NSETS ), M( NSETS ),
     &         FREQ, TITENV )
         ENDIF
 2000 ENDDO   ! next node

      RETURN
      END
