PROGRAM WRM1X1S1
C=======================================================================
C SOLUTION OF D2U/DXDX + ALPHASQ*U = 0 USING WEIGHTED RESIDUAL METHOD
C         WITH AN APPROXIMATING FUNCTION OF U(X)=F0(X)+A1*F1(X)
C              AND BOUNDARY CONDITIONS OF U(0)=U0. & DU/DX(L)=S;
C -------------- VARIABLE DEFNITION ----------- NOV/2023  EIJI FUKUMORI 
C   XST & XEN: INTEGRATION LIMITS; NSEG: NUMBER OF SEGMENTS IN LIMITS;
C   UNKNOWN COEFFICENT (A1) IN THE APPROXIMATING FUNCTION IS EVALUATED
C    BY THE FOLLOWING EQUATION:           B1 * A1 + C1 = 0.
C--------------------      F0 = U0 + S*RL*(X/RL)
C--------------------      F1 = (X/RL)*(1.D0-X/RL) + (X/RL)
C=======================================================================
      IMPLICIT REAL * 8 ( A-H , O-Z )
      PARAMETER ( INTEPT = 3, NSEG=100, MULTI=10 )
      DIMENSION SAI(INTEPT) , W(INTEPT)
      COMMON / DEL / DELTAX     /DOMAIN/ RL      /BORDER/ U0, S
      COMMON / COFF / ALPHASQ
      EXTERNAL F0, F1
C=======================================================================
C            THREE-SAMPLING-POINT GAUSS INTEGRATION METHOD
C            N: NUMBER OF SAMPLING POINTS IN EACH SEGMENET
C   SAI(I) & W(I): NON-DIMENSIONALIZED COORDINATE & WEIGHTING FACTOR
      DATA SAI/-0.7745966692415D0,0.0000000000000D0, 0.7745966692415D0/
      DATA W  / 0.5555555555555D0,0.8888888888888D0, 0.5555555555555D0/
C=======================================================================
C------------------- MATERIAL DATA AND BOUNDARY VALUES
      ALPHASQ=1.D0
      XST=0.D0
      XEN=0.5D0
      U0 = 1.D0
      S  = 0.D0
C=======================================================================
      WRITE(*,fmt='(a)', advance='no') 'TYPE IN S='
      READ(*,*) S
C=======================================================================
      OPEN ( 1, FILE='WRM1X1S1.FEM',STATUS='UNKNOWN' )
      WRITE(1,*)'==== ONE DIMENSIONAL HELMHOLTZ EQUATION ===='
      WRITE(1,*)'==== DIRICHLET ------- NEUMANN PROBLEM ===='
      WRITE(1,*)'---- GALERKINS WEIGHTING FUNCTION----'
      WRITE(1,*)'# OF GL INTEGRATION SAMPLING POINTS =',INTEPT
      WRITE(1,*)'APPROXIMATING FUNCTION: U(X) = F0(X) + A1*F1(X)'
      WRITE(1,*)'WHERE F0(X) = U0+SL(X/L)'
      WRITE(1,*)'AND   F1(X) = (X/L)*(1-X/L) + (X/L)'
C=======================================================================
C   DELTAX: SPACIAL DEFERENTIAL LENGTH FOR DERIVATIVE EVALUATION.
      RL = XEN - XST
      DELTAX = RL / ( MULTI * NSEG )
C=======================================================================
      WRITE(1,*)'LENGTH OF DOMAIN =',RL
      WRITE(1,*)'NUMBER OF SEGMENTS FOR INTEGRATION =', NSEG
      WRITE(1,*)'DX FOR DERIVATIVE EVALUATION =', DELTAX
C=======================================================================
      WRITE(1,*)'X-COORDINATE OF LEFT  END BOUNDARY =',XST
      WRITE(1,*)'X-COORDINATE OF RIGHT END BOUNDARY =',XEN
      WRITE(1,*)'ALPHASQ =', ALPHASQ
      WRITE(1,*)'NUMBER OF SEGMENTS =', NSEG
      WRITE(1,*)'DX FOR DERIVATIVE EVALUATION =', DELTAX
C----------------BOUNDARY CONDITIONS
      WRITE(1,*) 'U(X) AT X=0 =',U0
      WRITE(1,*) 'S    AT X=L =', S
C=======================================================================
C                COMPUTATION OF H(F0,F1) AND H(F1,F1)
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, INTEPT, SAI, W, F0, F1, C1 )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, INTEPT, SAI, W, F1, F1, B1 )
C=======================================================================
C      EVALUATION OF UNKNOWN A1 IN THE APPROXIMATING FUNCTION U(X)
      A1 = -(C1+S) / B1
C=======================================================================
C                     PRINTING RESULTS
      WRITE (1,*) 'H(F1,F1)=',B1
      WRITE (1,*) 'H(F0,F1)=',C1
      WRITE (1,*) 'A1 = -(H(F0,F1)+S)/H(F1,F1) =', A1
C=======================================================================
      CALL OUTPUT ( XST, XEN, A1 )
      CLOSE (1)
      STOP 'NORMAL TERMINATION'
      END
C
C
      SUBROUTINE INTE ( ALPHASQ,XST,XEN,NSEG,INTEPT,SAI,W, G1,G2,TOTAL )
      IMPLICIT REAL * 8 ( A-H , O-Z )
      DIMENSION SAI(INTEPT) , W(INTEPT)
      EXTERNAL G1, G2
      TOTAL = 0.D0
      DX = ( XEN - XST ) / NSEG
      DO I = 1 , NSEG
      X1 = DX*(I-1) + XST
      X2 = X1 + DX
      SUM = 0.D0
      SH  = ( X2 - X1 ) / 2.D0
      AVE = ( X1 + X2 ) / 2.D0
      DO J = 1 , INTEPT
      X = SH * SAI(J) + AVE
      SUM = SUM + (-DERIV(G1,X)*DERIV(G2,X)+ALPHASQ*G1(X)*G2(X)) * W(J)
      END DO
      TOTAL = TOTAL + SH * SUM
      END DO
      RETURN
      END
C
C
      FUNCTION F0(X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON /DOMAIN/ RL
      COMMON /BORDER/ U0, S
      F0 = U0 + S*RL*(X/RL)
      RETURN
      END
C
C
      FUNCTION F1(X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON /DOMAIN/ RL
      F1 = (X/RL)*(1.D0-X/RL) + (X/RL)
      RETURN
      END
C
C
      FUNCTION DERIV(F,X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON / DEL / DELTAX
      EXTERNAL F
      DERIV = ( F(X+DELTAX) - F(X-DELTAX) ) / ( 2.D0*DELTAX )
      RETURN
      END
C
C
      SUBROUTINE OUTPUT ( XST,XEN,A1 )
      IMPLICIT REAL * 8 ( A-H , O-Z )
      EXTERNAL F0, F1
      NDIV = 10
      DX = ( XEN - XST ) / NDIV
      WRITE(1,*)'X-COORDINATE U(X)-RIGHT-NEUMANN DU/DX EXACT'
      DO I = 1 , NDIV+1
      X = DX*(I-1) + XST
      UX = F0(X) + A1*F1(X)
      DUDX = DERIV(F0,X)+A1*DERIV(F1,X)
      WRITE(1,*) X,  UX, DUDX, EXACT(X)
      END DO
      RETURN
      END
C
C
      FUNCTION EXACT (X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON / DEL / DELTAX     /DOMAIN/ RL      /BORDER/ U0, S
      COMMON / COFF / ALPHASQ
C
      A = U0
      AL = DSQRT (ALPHASQ)
      B = (AL*U0*DSIN(AL*RL)+S)/(AL*DCOS(AL*RL))
      EXACT = A*DCOS(AL*X) + B*DSIN(AL*X)
      RETURN
      END