c *** Integration of Bessel Function J0 using the original and Hisada's methods

      program integp

      implicit real*8(a-h,o-z)
      complex*16 s1(100),s2,s3,cc,mu,im,oc,fa,fb,fc,OM
      real*4 a

      write(*,*)'Enter iom'
      read(*,*)iom
      write(*,*)'Enter nom'
      read(*,*)nom
      write(*,*)'Enter dom'
      read(*,*)dom
      write(*,*)'Enter omi'
      read(*,*)omi
      write(*,*)'Enter Velocity'
      read(*,*)c
      write(*,*)'Enter Q'
      read(*,*)q
      cc=c*dcmplx(1.d0,-0.5d0/q)
      write(*,*)'Enter Horizontal Distance'
      read(*,*)r

      im=dcmplx(0.d0,1.d0)
      DO 100 IO=IOM,NOM
        OM=DCMPLX(DOM*IO,OMI)
        oc=om/cc
        s1(IO)=cdexp(im*oc*r)/r
 100  CONTINUE
C      
      OPEN(UNIT=8,FILE='grfault.dxyz',STATUS='UNKNOWN')
        WRITE(8,2500)DOM,IOM,NOM,omi
        WRITE(8,2510)1
        DO 500 IOB=1,1
          DO 510 K=IOM,NOM
            S1(K)=DCMPLX(DREAL(S1(K)),-DIMAG(S1(K)))
            WRITE(8,2520)S1(K),S1(K),S1(K)
  510     CONTINUE
  500   CONTINUE
      CLOSE(8)
C
 2500 FORMAT(1X,F10.6,2I10,F10.6)
 2510 FORMAT(1X,I10,///)
 2520 FORMAT(1X,6E12.5)

      STOP
      end
