C
      PROGRAM MHANK
C
      PARAMETER(MOM=200,MOB=21)
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VS,SM,WB,SS,FHANK,VV(MOM,MOB),CJ0,CN0,VH(MOM,MOB)
      DIMENSION XOB(MOB)
      COMMON /DATA/ OM,X,WB,SM
      CHARACTER FNAME*15
      EXTERNAL FHANK
C
C *** HEAD ***
C
      WRITE(*,*)
      WRITE(*,1000)
 1000 FORMAT(/' From MHANK : SOURCE POINT is located at (0,0)'
     * /' Z OF OB IS ASSUMED TO BE 0'
     * /' EXACT GREEN FUNTION IN A 2-D SH MEDIUM (= HANKEL FUNC)')
C
      Q=100.D0
      DN=2.5D0
      V=3000.D0
      WRITE(*,1010)Q,DN,V
 1010 FORMAT( ' Q,DN,V    =',4F8.1)
C
C *** DATA INPUT ***
C
   10 WRITE(*,*)'OMEGA NUMBER?'
      READ(*,*)NOM
C
      IF(MOM.LT.NOM)THEN
        WRITE(*,*)'FROM MAIN MHANK'
        WRITE(*,*)'MOM < NOM RE-ENTRER !!!'
        WRITE(*,*)'MOM,NOM =',MOM,NOM
        GOTO 10
      END IF
C
      WRITE(*,*)'FIRST OM?'
      READ(*,*)FOM
      WRITE(*,*)'DELTA OM?'
      READ(*,*)DOM
C
   20 WRITE(*,1110)
      READ(*,*)NOB
      IF(NOB.GT.MOB)THEN
        WRITE(*,*)'NOB > MOB !!! : Re-Enter NOB : MOB=',MOB
        GOTO 20
      END IF
      IF(NOB.EQ.0)GOTO 20
C
      DO 100 I=1,NOB
        WRITE(*,1120)I
        READ(*,*)XOB(I)
  100 CONTINUE
C
 1110 FORMAT(/' Enter Number of Observation Point (Z=0) : ',\)
 1120 FORMAT( ' Enter X-axis of Observation Point of Num. =',I3,' : ',\)
C
      WRITE(*,1210)
      READ(*,*)SCC
      WRITE(*,1220)
      READ(*,*)NS
C
 1210 FORMAT(/' Enter the lowest phase velocity : ',\)
 1220 FORMAT( ' Enter div. num. for intergral   : ',\)
C
C *** DATA CAL ***
C
      VS=V*DCMPLX(1.D0,0.5D0/Q)
      SM=DN*VS*VS
C
C *** GREEN FUNC ***
C
      DO 200 IOB=1,NOB
C
        X=XOB(IOB)
C
        DO 210 IOM=1,NOM
C
          OM=FOM+DOM*DFLOAT(IOM-1)
          WB=OM*OM/(VS*VS)
C
          WNA=0.D0
          WNB=OM/SCC
          IF(OM.EQ.0.D0)WRITE(*,*)'OM=0 !!!!!!!!!!!!!!!!!!!!'
C
          CALL NEWTON(FHANK,WNA,WNB,NS,SS,IER)
          IF(IER.NE.0)WRITE(*,*)' IOB,IOM,IER',IOB,IOM,IER,'!!!!!!!!!!'
C
          VV(IOM,IOB)=-DCMPLX(0.D0,1.D0)*SS/(2.D0*3.1415927D0*SM)
C
          CALL CJN0(X*WB,CJ0,CN0)
          VH(IOM,IOB)=-DCMPLX(CN0,CJ0)/(4.D0*SM)
C
  210   CONTINUE
C
  200 CONTINUE
C
C *** OUTPUT OF DATA ***
C
        WRITE(*,1400)DOM,NOM
        WRITE(*,1401)
        DO 300 IOB=1,NOB
        WRITE(*,1410)XOB(IOB)
          DO 310 K=1,NOM,4
            WRITE(*,1420)(VV(K+J,IOB),J=0,3)
  310     CONTINUE
  300   CONTINUE
C
        WRITE(*,1402)
        DO 400 IOB=1,NOB
        WRITE(*,1410)XOB(IOB)
          DO 410 K=1,NOM,4
            WRITE(*,1420)(VH(K+J,IOB),J=0,3)
  410     CONTINUE
  400   CONTINUE
C
 1400 FORMAT(/' DOM,NOM=',F6.2,I5)
 1401 FORMAT(/' V(INTERGRAL):'/)
 1402 FORMAT(/' V(EXACT):'/)
 1410 FORMAT(/' X=',F10.2)
 1420 FORMAT(8D10.3)
C
      STOP
      END
C
C
C
      COMPLEX FUNCTION FHANK*16 (WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 SM,WB,CM
      COMMON /DATA/ OM,X,WB,SM
C
      W2=WN*WN
      CM=CDSQRT(WB-W2)
      IF(IMAG(CM).GT.0.D0)CM=-CM
C
      FHANK=DCOS(WN*X)/CM
C
      RETURN
      END
C
C
C
      SUBROUTINE NEWTON (FUNC,A,B,N,S,IER)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 FUNC,Y,FX,Q,S
C
      IER = 0
      IF( N-1 ) 100,110,110
  100 IER=2000
      S=(0.D0,0.D0)
      GO TO 9999
      IF(A.EQ.B)THEN
        IER=2000
        GO TO 9999
      END IF
C
  110 H=(B-A)/N
      Y   = FUNC(A)
      Q=(0.D0,0.D0)
C
      DO  120 I=1,N
      X=A+H*I
      FX = FUNC(X)
      Q   = Q+H*(Y+FX)/2.D0
      Y   = FX
  120 CONTINUE
C
      S = Q
 9999 CONTINUE
      IER=0
      RETURN
C
      END
