      PROGRAM HELM4Q
C=======================================================================
C           2-DIM FEM PROGRAM FOR SOLVING HELMHOLZ EQUATION
C                  USING UPPER HALF BANDED MATRIX
C              EQUATION: DP/DX + DP/DY + ALPHA * P = 0.
C                ELEMENT : 4-NODED ISO-PARAMETERIC
C         NUMBERING ORDER: (-1,-1),(+1,-1),(+1,+1),(-1,+1)
C   ORIGINAL:1984 EIJI FUKUMORI BUFFALO NY & REVISED: 1994 ACHI
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=4,MXE=30000,MXN=30000,MXB=3000,INTEPT=2,MXW=150 )
      PARAMETER ( ND3=(ND*ND-ND)/2+ND )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN), IBND(MXB),
     * ITYPE(MXB), BVALUE(MXB),RK(MXN,MXW), RHS(MXN),SK(ND,ND), 
     * X(ND), Y(ND),BX(2,ND),SAI(INTEPT), W(INTEPT),SS(ND,ND),
     * BPP(2,ND,INTEPT,INTEPT), SF(ND,INTEPT,INTEPT)
      DIMENSION GSTIFF(MXE,ND3), SOURCE(MXE,ND3)
      CHARACTER INPFILE*12
      WRITE (*,*)' HELM4Q.FOR SOLVER'
C=======================================================================
      CALL GRULE ( INTEPT, SAI, W   )
      CALL DERIV ( ND, INTEPT, X, Y, SAI, BPP )
      CALL SHAPEF( ND, INTEPT, X, SAI, SF )
C=======================================================================
      INPFILE = 'HEL4DATA.QQQ'
      WRITE(*,*)' INPUT FILE NAME =', INPFILE
      CALL INPUT ( ND,MXE,MXN,MXB,NE,NNODE,NB,ALPHA,NODEX,
     * XCOORD,YCOORD,IBND,ITYPE,BVALUE, INPFILE )
C=======================================================================
      CALL BANDWID ( MXE, ND, NE,NODEX, NBW )
      IF ( NBW .GT. MXW ) STOP 'Error ... NBW > MXW'
C=======================================================================
      CALL GSM ( INTEPT,ND,ND3,MXE,MXN,NE,BPP,W,BX,SK,NODEX,
     * XCOORD,YCOORD,GSTIFF )
C=======================================================================
      CALL GDM ( MXE,MXN,INTEPT,ND,ND3,BPP,W,NE,SF,
     * XCOORD,YCOORD,NODEX,SS, SOURCE )
C=======================================================================
      CALL SUMUP ( ND,ND3,MXE,MXN,MXW,NE,NBW,NNODE,ALPHA,
     * NODEX,GSTIFF,SOURCE,SS,RK )
C=======================================================================
      CALL FORM (MXN,MXB,MXW,NNODE,NB,NBW,RK,RHS,ITYPE,BVALUE,IBND)
C=======================================================================
      CALL SYSTEMQ ( MXN, MXW, NNODE, NBW, RK, RHS )
C=======================================================================
      CALL PRINT ( MXN, NNODE , RHS )
      CALL PRINTG ( MXN, NNODE, RHS, XCOORD, YCOORD )
      STOP
      END
C
C
      SUBROUTINE GSM ( INTEPT,ND,ND3,MXE,MXN,NE,BPP,W,BX,SS,NODEX,
     * XCOORD,YCOORD,GSTIFF )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),GSTIFF(MXE,ND3),
     * SS(ND,ND), BX(2,ND),W(INTEPT), BPP(2,ND,INTEPT,INTEPT)
C
      DO  IEL = 1 ,NE
      DO I = 1 , ND
      DO J = I , ND
      SS(I,J) = 0.
      END DO
      END DO
      DO  K = 1 , INTEPT
      DO  L = 1 , INTEPT
      WEIGHT = W(K) * W(L)
      YAC11 = 0.
      YAC12 = 0.
      YAC21 = 0.
      YAC22 = 0.
      DO I = 1 , ND
      YAC11 = YAC11 + BPP(1,I,K,L) * XCOORD(NODEX(IEL,I))
      YAC12 = YAC12 + BPP(1,I,K,L) * YCOORD(NODEX(IEL,I))
      YAC21 = YAC21 + BPP(2,I,K,L) * XCOORD(NODEX(IEL,I))
      YAC22 = YAC22 + BPP(2,I,K,L) * YCOORD(NODEX(IEL,I))
      END DO
      DETJ = YAC11 * YAC22 - YAC12 * YAC21
      Y11 =   YAC22 / DETJ
      Y12 = - YAC12 / DETJ
      Y21 = - YAC21 / DETJ
      Y22 =   YAC11 / DETJ
      DO J = 1 , ND
      BX(1,J) = Y11 * BPP(1,J,K,L) + Y12 * BPP(2,J,K,L)
      BX(2,J) = Y21 * BPP(1,J,K,L) + Y22 * BPP(2,J,K,L)
      END DO
      BETA = WEIGHT * DETJ
      DO I = 1 , ND
      DO J = I , ND
      SS(I,J) = SS(I,J) + ( BX(1,I)*BX(1,J)+BX(2,I)*BX(2,J) )*BETA
      END DO
      END DO
      END DO
      END DO
      M = 0
      DO K = 1 , ND
      DO L = K , ND
      M = M + 1
      GSTIFF(IEL,M) = SS(K,L)
      END DO
      END DO
      END DO
      RETURN
      END
C
C
      SUBROUTINE INPUT ( ND,MXE,MXN,MXB,NE,NNODE,NB,ALPHA,NODEX,
     * XCOORD,YCOORD,IBND,ITYPE,BVALUE, INPFILE )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBND(MXB),
     * ITYPE(MXB), BVALUE(MXB)
      LOGICAL YES
      CHARACTER INPFILE*12
      IR = 1
      INQUIRE ( FILE=INPFILE, EXIST=YES )
      IF ( YES ) OPEN ( IR, FILE=INPFILE, STATUS='OLD' )
      IF ( .NOT. YES ) STOP' NO INPUT FILE'
C========> PARAMETER
      READ (1,*) ALPHA
C========> ELEMENT
      READ (1,*) NE
      DO I = 1 , NE
      READ (1,*) IEL, ( NODEX(IEL,J), J = 1 , ND )
      END DO
C========> COORDINATES
      READ (1,*) NNODE
      DO I = 1 , NNODE
      READ (1,*) NODE, XCOORD(NODE) , YCOORD(NODE)
      END DO
C========> BOUNDARY CONDITIONS
      READ(1,*) NB
      DO I = 1 , NB
      READ (1,*) IBND(I) , ITYPE(I), BVALUE(I)
      END DO
      CLOSE (1)
      RETURN
      END
C
C
      SUBROUTINE PRINT ( MXN, NNODE, RHS )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION RHS(MXN)
      CHARACTER OUTFILE*12, EXFILE*3
      LOGICAL YES
      OUTFILE ='SOLUTION.BIN'
      EXFILE ='NEW'
      IW = 1
      INQUIRE ( FILE=OUTFILE, EXIST=YES )
      IF ( YES ) EXFILE='OLD'
      OPEN ( IW, FILE=OUTFILE, STATUS=EXFILE,FORM='UNFORMATTED' )
      WRITE (IW) (RHS(I),I=1,NNODE)
      CLOSE (IW)
      RETURN
      END
C
C
      SUBROUTINE PRINTG ( MXN, NNODE, RHS, XCOORD, YCOORD )
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( IY =50, LAYER=36 )
      DIMENSION RHS(MXN), XCOORD(MXN), YCOORD(MXN), IWAVE(110,110)
      CHARACTER OUTFILE*12, EXFILE*3, IWAVE*1, SYMBOL(LAYER)*1
      LOGICAL YES
      DATA SYMBOL /'0', '1', '2', '3', '4', '5', '6', '7', '8', '9',
     *             'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J',
     *             'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T',
     *             'U', 'V', 'W', 'X', 'Y', 'Z' /
      OUTFILE ='SOLUTION.FEM'
      EXFILE ='NEW'
      IW = 1
      INQUIRE ( FILE=OUTFILE, EXIST=YES )
      IF ( YES ) EXFILE='OLD'
      OPEN ( IW, FILE=OUTFILE, STATUS=EXFILE,FORM='FORMATTED' )
      CALL MINMAX ( MXN, NNODE, XCOORD, XMIN, XMAX )
      CALL MINMAX ( MXN, NNODE, YCOORD, YMIN, YMAX )
      CALL MINMAX ( MXN, NNODE, RHS, RMIN, RMAX )
      SPANX = XMAX - XMIN
      SPANY = YMAX - YMIN
      SPANR = RMAX - RMIN
      DO I = 1 , 110
      DO J = 1 , 110
      IWAVE(I,J) = ' '
      END DO
      END DO
C
      IX = SPANX/SPANY*IY
      DO NODE = 1 , NNODE
      I = (( XCOORD(NODE) - XMIN ) / SPANX ) * IX + 1
      J = (( YCOORD(NODE) - YMIN ) / SPANY ) * IY + 1
      K = (( RHS(NODE) - RMIN ) / SPANR)*LAYER + 0.5
      IF ( K .LE. 0 ) IWAVE(I,J) = SYMBOL(1)
      IF ( K .EQ. 1 ) IWAVE(I,J) = SYMBOL(2)
      IF ( K .EQ. 2 ) IWAVE(I,J) = SYMBOL(3)
      IF ( K .EQ. 3 ) IWAVE(I,J) = SYMBOL(4)
      IF ( K .EQ. 4 ) IWAVE(I,J) = SYMBOL(5)
      IF ( K .EQ. 5 ) IWAVE(I,J) = SYMBOL(6)
      IF ( K .EQ. 6 ) IWAVE(I,J) = SYMBOL(7)
      IF ( K .EQ. 7 ) IWAVE(I,J) = SYMBOL(8)
      IF ( K .EQ. 8 ) IWAVE(I,J) = SYMBOL(9)
      IF ( K .EQ. 9 ) IWAVE(I,J) = SYMBOL(10)
      IF ( K .EQ. 10) IWAVE(I,J) = SYMBOL(11)
      IF ( K .EQ. 11) IWAVE(I,J) = SYMBOL(12)
      IF ( K .EQ. 12) IWAVE(I,J) = SYMBOL(13)
      IF ( K .EQ. 13) IWAVE(I,J) = SYMBOL(14)
      IF ( K .EQ. 14) IWAVE(I,J) = SYMBOL(15)
      IF ( K .EQ. 15) IWAVE(I,J) = SYMBOL(16)
      IF ( K .EQ. 16) IWAVE(I,J) = SYMBOL(17)
      IF ( K .EQ. 17) IWAVE(I,J) = SYMBOL(18)
      IF ( K .EQ. 18) IWAVE(I,J) = SYMBOL(19)
      IF ( K .EQ. 19) IWAVE(I,J) = SYMBOL(20)
      IF ( K .EQ. 20) IWAVE(I,J) = SYMBOL(21)
      IF ( K .EQ. 21) IWAVE(I,J) = SYMBOL(22)
      IF ( K .EQ. 22) IWAVE(I,J) = SYMBOL(23)
      IF ( K .EQ. 23) IWAVE(I,J) = SYMBOL(24)
      IF ( K .EQ. 24) IWAVE(I,J) = SYMBOL(25)
      IF ( K .EQ. 25) IWAVE(I,J) = SYMBOL(26)
      IF ( K .EQ. 26) IWAVE(I,J) = SYMBOL(27)
      IF ( K .EQ. 27) IWAVE(I,J) = SYMBOL(28)
      IF ( K .EQ. 28) IWAVE(I,J) = SYMBOL(29)
      IF ( K .EQ. 29) IWAVE(I,J) = SYMBOL(30)
      IF ( K .EQ. 30) IWAVE(I,J) = SYMBOL(31)
      IF ( K .EQ. 31) IWAVE(I,J) = SYMBOL(32)
      IF ( K .EQ. 32) IWAVE(I,J) = SYMBOL(33)
      IF ( K .EQ. 33) IWAVE(I,J) = SYMBOL(34)
      IF ( K .EQ. 34) IWAVE(I,J) = SYMBOL(35)
      IF ( K .GE. 35) IWAVE(I,J) = SYMBOL(36)
      END DO
C
      WRITE(IW,*) ' MAX=',RMAX , '     MIN=',RMIN
      DO I = 1 , IX
      WRITE (IW,'(1X,110A1)') ( IWAVE(I,J), J=1, IY )
      END DO
C
      WRITE (IW,*) ' SYMBOL     CENTRAL VALUE'
      DO K = 1 , LAYER
      VALUE = K*SPANR/LAYER + RMIN
      WRITE(IW,*) SYMBOL(K) , VALUE
      END DO
C
      WRITE(IW,*)' POTENTIAL AT NODAL POINTS'
      WRITE (IW,'(I6,G17.7)') (I,RHS(I),I=1,NNODE)
      CLOSE (IW)
      RETURN
      END
C
C
      SUBROUTINE MINMAX ( MXN, NNODE, X, XMIN, XMAX )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION X(MXN)
      XMAX = X(1)
      XMIN = X(1)
      DO I = 1 , NNODE
      IF ( X(I) .GT. XMAX ) XMAX = X(I)
      IF ( X(I) .LT. XMIN ) XMIN = X(I)
      END DO
      RETURN
      END
C
C
      SUBROUTINE BANDWID ( MXE, ND, NE, NODEX , NBW  )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND)
      NBW = 0
      DO  I = 1 , NE
      DO  J = 1 , ND-1
      DO  K = J+1 , ND
      NBW = MAX0(NBW,IABS(NODEX(I,J)-NODEX(I,K)))
      END DO
      END DO
      END DO
      NBW = NBW + 1
      write(6,*) ' bandwidth =', nbw
      RETURN
      END

C
C
      SUBROUTINE SYSTEMQ ( MXN, MXW, NNODE, NBW, A, B )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION A(MXN,MXW) , B(MXN)
C---------- GAUSS ELIMINATION --------------
      DO N = 1 , NNODE-1
      DO L = 2 , MIN((NNODE-N+1),NBW)
      C = A(N,L) / A(N,1)
      IF ( C .NE. 0. ) THEN
      I = N + L - 1 
      DO K = L , NBW
      J = K-L+1
      A(I,J) = A(I,J) - C * A(N,K)
      END DO
      A(N,L) = C
      B(I) = B(I) - C * B(N)
      END IF
      END DO
      B(N) = B(N) / A(N,1)
      END DO
C------ SOLUTION OF LAST UNKNOWN ---------
      B(NNODE) = B(NNODE) / A(NNODE,1)
C---------- BACK SUBSTITUTION -------------
      DO I = NNODE-1 , 1 , -1
      DO J = 2 , MIN((NNODE-I+1),NBW)
      B(I) = B(I)-A(I,J)*B(I+J-1) 
      END DO
      END DO
      RETURN
      END 
C
C
      SUBROUTINE FORM (MXN,MXB,MXW,NNODE,NB,NBW,A,RHS,IBC,BV,IBND)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION RHS(MXN),A(MXN,MXW),IBC(MXB),BV(MXB),IBND(MXB)
C----- REFORMATION OF VECTOR RHS DUE TO FIRST KIND B.C. 
      DO K = 1 , NB
      IF ( IBC(K) .EQ. 1 ) THEN
C          CASE OF DIRICHLET BOUNDARY CONDITION
      I = IBND(K) 
      DO J = 2 , NBW 
      I = I - 1 
      IF ( I.GT. 0 ) RHS(I) = RHS(I) - BV(K) * A(I,J)
      END DO
      I = IBND(K) 
      DO J = 2 , NBW 
      I = I + 1 
      IF ( I .LE. NNODE ) RHS(I) = RHS(I) - BV(K) * A(IBND(K),J)
      END DO
      END IF
      END DO
C-----REFORMATION OF MATRIX A.
      DO  K = 1 , NB
      I = IBND (K)
      IF ( IBC(K) .EQ. 1 ) THEN
C          CASE OF DIRICHLET BOUNDARY CONDITION
      RHS(I) = BV(K) 
      A(I,1) = 1. 
      DO J = 2 , NBW 
      L = I - J + 1 
      A(I,J) = 0. 
      IF ( L  .GT. 0 ) A(L,J) = 0. 
      END DO
      END IF
C          CASE OF NUEMANN BOUNDARY CONDITION
      IF ( IBC(K) .EQ. 2 ) RHS (I) = RHS(I) - BV(K)
      END DO
      RETURN
      END 
C
C
      SUBROUTINE GRULE ( N , SAI , W )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(N) , W(N)
      IF ( N .LT. 2 ) STOP'N<2'
      IF ( N .GT. 6 ) STOP'N>6'
      IF ( N .EQ. 1 ) THEN
      SAI(1) = 0.
      W(1) = 2.
      RETURN
      END IF
C
      IF ( N .EQ. 2 ) THEN
      SAI(1) = 0.5773502692
      W(1) = 1.
      END IF
C
      IF ( N .EQ. 3 ) THEN
      SAI(1) = 0.7745966692
      SAI(2) = 0.
      W(1) = 5./ 9.
      W(2) = 8./ 9.
      END IF
C
      IF ( N .EQ. 4 ) THEN
      SAI(1) = 0.3399810436
      SAI(2) = 0.8611363116
        W(1) = 0.6521451549
        W(2) = 0.3478548451
      END IF
C
      IF ( N .EQ. 5 ) THEN
      SAI(1) = 0.9061798459
      SAI(2) = 0.5384693101
      SAI(3) = 0.
        W(1) = 0.2369268850
        W(2) = 0.4786286705
        W(3) = 5.12 / 9.
      END IF
C
      IF ( N .EQ. 6 ) THEN
      SAI(1) = 0.2386191861
      SAI(2) = 0.6612093865
      SAI(3) = 0.9324695142
        W(1) = 0.4679139346
        W(2) = 0.3607615730
        W(3) = 0.1713244924
      END IF
C
      NN = N / 2
      DO  I = 1 , NN
      J = N - I + 1
      SAI(J) = - SAI(I)
      W(J) = W(I)
      END DO
      RETURN
      END
C
C
      SUBROUTINE DERIV ( ND, INTEPT, F0, F1, SAI, BPP )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(INTEPT),BPP(2,ND,INTEPT,INTEPT),F0(ND),F1(ND)
      DO  K = 1 , INTEPT
      E1 = SAI (K)
      DO  L = 1 , INTEPT
      E2 = SAI (L)
COMPUTATION OF BP(J) = DN(J) / DETA1
      CALL ISOPARA ( ND, E1+0.5 , E2 , F1 )
      CALL ISOPARA ( ND, E1-0.5 , E2 , F0 )
      DO I = 1 , ND
      BPP(1,I,K,L) = F1(I) - F0(I)
      END DO
COMPUTATION OF BP(J)= DN(J)/DETA2
      CALL ISOPARA ( ND, E1 , E2+0.5 , F1 )
      CALL ISOPARA ( ND, E1 , E2-0.5 , F0 )
      DO I = 1 , ND
      BPP(2,I,K,L) = F1(I) - F0(I)
      END DO
      END DO
      END DO
      RETURN
      END
C
C
      SUBROUTINE ISOPARA ( ND, E1 , E2 , F )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION F(ND)
      F(1) = 0.25*(1.- E1)*(1.- E2)
      F(2) = 0.25*(1.+ E1)*(1.- E2)
      F(3) = 0.25*(1.+ E1)*(1.+ E2)
      F(4) = 0.25*(1.- E1)*(1.+ E2)
      RETURN
      END
C
C
      SUBROUTINE SHAPEF ( ND, INTEPT, F, SAI, SF )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION F(ND) , SAI(INTEPT) , SF(ND,INTEPT,INTEPT)
      DO  K = 1 , INTEPT
      E1 = SAI (K)
      DO  L = 1 , INTEPT
      E2 = SAI( L )
      CALL ISOPARA ( ND, E1 , E2 , F )
      DO  I = 1 , ND
      SF(I,K,L) = F(I)
      END DO
      END DO
      END DO
      RETURN
      END
C
C
      SUBROUTINE GDM ( MXE,MXN,INTEPT,ND,ND3,BPP,W,NE,SF,
     * XCOORD,YCOORD,NODEX,SS, SOURCE )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),SOURCE(MXE,ND3),
     * SS(ND,ND),BPP(2,ND,INTEPT,INTEPT),W(INTEPT),SF(ND,INTEPT,INTEPT)
C
      DO  IEL = 1 ,NE
      DO I = 1 , ND
      DO J = I , ND
      SS(I,J) = 0.
      END DO
      END DO
C
      DO  K = 1 , INTEPT
      DO  L = 1 , INTEPT
      WEIGHT = W(K) * W(L)
      YAC11 = 0.
      YAC12 = 0.
      YAC21 = 0.
      YAC22 = 0.
      DO I = 1 , ND
      YAC11 = YAC11 + BPP(1,I,K,L) * XCOORD(NODEX(IEL,I))
      YAC12 = YAC12 + BPP(1,I,K,L) * YCOORD(NODEX(IEL,I))
      YAC21 = YAC21 + BPP(2,I,K,L) * XCOORD(NODEX(IEL,I))
      YAC22 = YAC22 + BPP(2,I,K,L) * YCOORD(NODEX(IEL,I))
      END DO
      DETJ = YAC11 * YAC22 - YAC12 * YAC21
      BETA = WEIGHT * DETJ
C
      DO I = 1 , ND
      DO J = I , ND
      SS(I,J) = SS(I,J) + SF(I,K,L) * SF(J,K,L) * BETA
      END DO
      END DO
C
      END DO
      END DO
C
      M = 0
      DO K = 1 , ND
      DO L = K , ND
      M = M + 1
      SOURCE(IEL,M) = SS(K,L)
      END DO
      END DO
C
      END DO
      RETURN
      END
C
C
      SUBROUTINE SUMUP ( ND,ND3,MXE,MXN,MXW,NE,NBW,NNODE,ALPHA,
     * NODEX,GSTIFF,SOURCE,SS,RK )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SS(ND,ND),RK(MXN,MXW),GSTIFF(MXE,ND3),SOURCE(MXE,ND3),
     * NODEX(MXE,ND)
C
      DO I = 1 , NNODE
      DO J = 1 , NBW
      RK(I,J) = 0.
      END DO
      END DO
C
      DO  IEL = 1 , NE
      DO K = 1 , ND
      DO L = 1 , ND
      SS(K,L) = 0.
      END DO
      END DO
      M = 0
      DO K = 1 , ND
      DO L = K , ND
      M = M + 1
      SS(K,L) = - GSTIFF(IEL,M) + ALPHA*SOURCE(IEL,M)
      SS(L,K) = SS(K,L)
      END DO
      END DO
C----ASSEMBLY
      DO K = 1 , ND
      I = NODEX(IEL,K)
      DO L = 1 , ND
      J = NODEX(IEL,L) - I + 1
      IF ( J .GT. 0 ) RK(I,J) = RK(I,J) + SS(K,L)
      END DO
      END DO
      END DO
      RETURN
      END
