PROGRAM HELM1DIM
C********************************************************************
C                EXACT SOLUTION TO 1DIM HELMHOLTZ EQUATION
C                        D2U/DXDX + ALPHASQ*U = 0
C                       -- VARIABLE DEFINITION --
C                   IBC1 = BOUNDARY CONDITION AT X=0
C                    BVALUE1 = BOUNDARY VALUE AT X=0
C                   IBC2 = BOUNDARY CONDITION AT X=L
C                    BVALUE2 = BOUDARY VALUE AT X=L
C              EIJI FUKUMORI AT AICHI JAPAN   21 DEC 1998
C********************************************************************
      PARAMETER ( NSEG=20 )
      COMPLEX R1,R2,A11,A12,A21,A22, Z1,Z2
C------- A*R**2 + B*R + C = 0
      CALL INPUT ( IBC1,BVALUE1,IBC2,BVALUE2,DL,ALPHASQ )
      A = 1.
      B = 0.
      C = ALPHASQ
      D = B*B-4.*A*C
      R1 = ( -B-CSQRT(CMPLX(D,0.)) ) / (2.*A)
      R2 = ( -B+CSQRT(CMPLX(D,0.)) ) / (2.*A)
C
C-------- IMPLEMENTATION OF BOUNDARY CONDITIONS
C-- AT X=0.
      CALL BCIMPLE ( 0.,IBC1,R1,R2,A11,A12 )
C-- AT X=DL
      CALL BCIMPLE ( DL,IBC2,R1,R2,A21,A22 )
C
C------- DETERMINATION OF Z1 AND Z2
      Z1 = ( BVALUE1*A22 - A12*BVALUE2 ) / ( A11*A22 - A12*A21 )
      Z2 = ( BVALUE2*A11 - A21*BVALUE1 ) / ( A11*A22 - A12*A21 )
C
      CALL OUTPUT ( NSEG,DL,Z1,Z2,R1,R2 )
      STOP
      END
C
C
      SUBROUTINE BCIMPLE ( X,IBC,R1,R2,A1,A2 )
      COMPLEX R1,R2,A1,A2
C-------DIRICHLET
      IF ( IBC .EQ. 1 ) THEN
      A1 = CEXP(R1*X)
      A2 = CEXP(R2*X)
      ENDIF
      IF ( IBC .EQ. 2 ) THEN
C-------NEUMANN
      A1 = R1*CEXP(R1*X)
      A2 = R2*CEXP(R2*X)
      END IF
      RETURN
      END
C
C
      SUBROUTINE INPUT ( IBC1,BVALUE1,IBC2,BVALUE2,DL,ALPHASQ )
      WRITE (*,110)
  110 FORMAT( ' TYPE IN LENGTH OF DOMAIN = ' $ )      
      READ (*,*) DL
      IF ( DL .LE. 0. ) STOP 'EXECUTION TERMINATED'
C
      WRITE (*,120)
  120 FORMAT( ' TYPE IN SQUARE OF ALPHA = ' $ )      
      READ (*,*) ALPHASQ
      IF (ALPHASQ .LE. 0. ) STOP 'EXECUTION TERMINATED'
C
      WRITE (*,130)
  130 FORMAT( ' TYPE IN BOUNDARY CONDITION AT X=0 = ' $ )
      READ (*,*) IBC1
      WRITE (*,140)
  140 FORMAT( ' TYPE IN BOUNDARY VALUE AT X=0 = ' $ )
      READ (*,*) BVALUE1
C
      WRITE (*,150)
  150 FORMAT( ' TYPE IN BOUNDARY CONDITION AT X=L = ' $ )
      READ (*,*) IBC2
      WRITE (*,160)
  160 FORMAT( ' TYPE IN BOUNDARY VALUE AT X=L = ' $ )
      READ (*,*) BVALUE2
      RETURN
      END
C
C
      SUBROUTINE OUTPUT ( NSEG,DL,Z1,Z2,R1,R2 )
      COMPLEX R1,R2, Z1,Z2, U
      OPEN ( 1, FILE='EXACT.DAT', STATUS='UNKNOWN')
      WRITE(1,*)'X U(REAL) U(IMAGINARY)'
      DO I = 1 , NSEG+1
      X = (I-1)*DL/NSEG
      U = Z1*CEXP(R1*X) + Z2*CEXP(R2*X)
      WRITE(1,'(3F13.8)') X, U
      END DO
      CLOSE (1)
      RETURN
      END