C     
c *** Green's functions using Kamae's method
C *** FOR BUGS OR COMMENTS, PLEASE E-MAIL to "hisada@cc.kogakuin.ac.jp",
C     OR MAIL TO "Yoshiaki Hisada, Dept. of Architecture,
C     Kogakuin Univeristy, Nishi-Shinjuku 1-24-2, Tokyo 163-91, Japan"
C    (Originally made in 1992, and modified in 1998)
c
c * Parameters : 
c     MOM : Number of Omega
c     MOB : NUmber of Observation Points
c
      PROGRAM kamae
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MOM=32800,MOB=13,MFLT=227)
      COMPLEX*16 GX(MOM,MOB),FX(MOM)
      DIMENSION FD(MFLT,MFLT),TD(MFLT,MFLT),XR(MOM)
     *         ,XF(MFLT,MFLT),YF(MFLT,MFLT),ZF(MFLT,MFLT)
C
      WRITE(*,*)
      WRITE(*,*)'"Green functions using the Kamae et.al. method'
      WRITE(*,*)' Slip Duration by Somerville et.al(1999)'
      WRITE(*,*)' Function F by Irikura et.al.(1997)'
      WRITE(*,*)'(Originally made in 1992, and modified in 2000)'
      write(*,*)
      WRITE(*,*)'"kamae2.in" is the input data'
      write(*,*)
      WRITE(*,*)
     *' For questions, e-mail to "hisada@cc.kogakuin.ac.jp".'
      WRITE(*,*)' This code is open for the academic use.'
C
      CALL INDAT(MOM,MOB,XF,YF,ZF,TD,FD,MFLT)
      CALL BRSAH(MFLT,MOM,MOB,XF,YF,ZF,FD,TD,GX,FX,XR)
      CALL OTDAT(MOM,MOB,GX)
C
      STOP
      END
C
c
C *** DATA INPUT ***
C
      SUBROUTINE INDAT(MOM,MOB,XF,YF,ZF,TD,FD,MFLT)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION FD(MFLT,MFLT),TD(MFLT,MFLT)
     *         ,XF(MFLT,MFLT),YF(MFLT,MFLT),ZF(MFLT,MFLT)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(13),YO(13),ZO(13),NOB
     *    /STRC/DNS,VP,VS,Qp,Qs,Qpn,Qsn
     *    /FLT/STR,DIP,FL,FW,Vr,Fmax,Rtime
     *    /IFLT/NSTR,NDIP,idum,niter
     *    /Hypo/xhypo,yhypo,zhypo 
     *    /SIGN/ISIGN
c
C *** PARAMETERS IN COMMOM STATEMENT ***
C
      NMOB=13
      IF(MOB.LT.NMOB)NMOB=MOB
C
C *** OMEGA DATA ***
C
      write(*,*)
      OPEN (UNIT=7,FILE='kamae2.in',STATUS='OLD')
C
      read(7,1100)dt,nt
      read(7,1110)period
 1100 format(/f8.0,i8)
 1110 format(f8.0)
      duration=dt*nt
      dom=datan(1.d0)*8.d0/duration
      iom=1
      nom=duration/period

      WRITE(10,*)'FROM kamae.in:'
      WRITE(10,*)
      WRITE(10,1130)nt,dt,period
      WRITE(10,1140)NOM,DOM
 1130 format(/' Num. of Time, Delta Time and Min. Period (s)=',i5,2f8.4)
 1140 format( ' Num. of Omega, and Delta Omega (rad/sec) =',i5,f10.7)
C
      IF(NOM.GT.MOM)THEN
        WRITE(10,*)' NOM>MOM !!!!!!!!!!!!!'
        STOP
      END IF
C
C *** MEDIUM DATA ***
C
      READ(7,1200)DNS,VP,QP,VS,QS
c      WRITE(*,2210)DNS,VP,Qp,VS,Qs
      READ(7,1210)Qpn,Qsn
c      WRITE(*,2220)Qpn,Qsn
 1200 FORMAT(/6f8.0)
 1210 FORMAT(2f8.0)
C
      WRITE(10,2200)
      WRITE(10,2210)DNS,VP,Qp,VS,Qs
      WRITE(10,2220)Qpn,Qsn
 2200 FORMAT(/' Density , Vp(m/s), Qp0     '
     *        ', Vs(m/s) , Qs0    ,Thickness(m)')
 2210 FORMAT(5F10.3,F12.3)
 2220 format(' n of Q*f^n for Qp and Qs =',2f8.2)
C
C *** FAULT DATA ***
C
      READ(7,1300)SLENGTH,DWIDTH
      write(*,*)'SLENGTH,DWIDTH ',SLENGTH,DWIDTH
      READ(7,1310)NSTR,NDIP
      aND=dsqrt(dfloat(ndip*nstr))
      write(*,*)'NSTR,NDIP,aND ',NSTR,NDIP,aND
      IF(NSTR.GT.MFLT)THEN
        WRITE(10,*)' NSTR.GT.MFLT !!!!!!'
        WRITE(*,*)' NSTR.GT.MFLT !!!!!!'
        STOP
      END IF
      IF(NDIP.GT.MFLT)THEN
        WRITE(10,*)' NDIP.GT.MFLT !!!!!!'
        WRITE(*,*)' NDIP.GT.MFLT !!!!!!'
        STOP
      END IF
      READ(7,1320)DELAY,niter
      write(*,*)'DELAY,Num. of Iterations ',DELAY,Niter
      READ(7,1330)SHYPO,DHYPO
      write(*,*)'SHYPO,DHYPO ',SHYPO,DHYPO
      READ(7,1330)xflt,yflt,zflt
      write(*,*)'xflt,yflt,zflt ',xflt,yflt,zflt
      READ(7,1330)STR,DIP
      write(*,*)'STR,DIP',STR,DIP
      str=str*0.01745329252d0
      dip=dip*0.01745329252d0
      READ(7,1340)VR,idum
      write(*,*)'Vr,idum ',Vr,idum
 1300 FORMAT(/2f10.0)
 1310 FORMAT(3I5)
 1320 FORMAT(F10.0,i10)
 1330 format(3f10.0)
 1340 format(f10.0,i10)
c
      


      read(7,*)
      aslip=0.0
      do i=1,ndip
      do j=1,nstr
        READ(7,*)FD(i,j)
c        write(*,*)'i,j,FD(i,j,k)',i,j,FD(i,j)
         aslip=aslip+FD(i,j)
      end do
      end do
      aslip=aslip/(ndip*nstr)
      

      SMo=Dns*Vs*Vs*aslip*SLENGTH*DWIDTH*10000000000.0d0
      rtime=2.03*(SMo**(1.0/3.0))/1000000000.0
      aMag=(dlog10(SMo)-17.0)/1.33
      Tw=10.0**(0.31*aMag-0.77)
      fc=10.0**((23.8-dlog10(SMo))/3.0)
      fmax=7310.0*SMo**(-0.12)
      write(*,*)
      write(*,*)'Avg Slip=',aslip
      write(*,*)'M,Mo    = ',aMag,SMo
      write(*,*)'fc,fmax = ',fc,fmax
      write(*,*)'rtime,Tw= ',rtime,Tw
      write(*,*)'Enter rtime (sec)'
      read(*,*)rtime
      write(10,*)
      write(10,*)'Avg Slip=',aslip
      write(10,*)'M,Mo    = ',aMag,SMo
      write(10,*)'fc,fmax = ',fc,fmax
      write(10,*)'rtime,Tw= ',rtime,Tw
c
c   * Locations of Sub_Fault Centers *  
C
      fl=slength/nstr
      fw=dwidth/ndip

      cdp=dcos(dip)
      sdp=dsin(dip)
      cst=dcos(str)
      sst=dsin(str)
      xhypo=xflt+shypo*cst+dhypo*cdp*sst
      yhypo=yflt+shypo*sst-dhypo*cdp*cst
      zhypo=zflt-dhypo*sdp
c							
c   * Dip direction *     
c
      do 200 i=1,ndip

        da=fw*(i-1)+fw/2
	da1=da-dhypo

c     * Strike direction *
        do 210 j=1,nstr   
   
          sa=fl*(j-1)+fl/2
	  sa1=sa-shypo
          dis=dsqrt(da1*da1+sa1*sa1)
          ta=dis/Vr
          xa=xflt+sa*cst+da*cdp*sst
          ya=yflt+sa*sst-da*cdp*cst
          za=zflt-da*sdp
          if(za.lt.0.d0)then
             write(*,*)'za < 0 !!'
             stop
          end if
          xf(i,j)=xa
          yf(i,j)=ya
          ZF(i,j)=za
          td(i,j)=delay+ta

  210   continue
  200 continue
C
      WRITE(10,2300)nstr,ndip
      WRITE(10,2305)STR,DIP
      WRITE(10,2315)slength,dwidth
      WRITE(10,2320)FL,FW    
      write(10,2322)shypo,dhypo
      write(10,2324)xhypo,yhypo,zhypo
      write(10,2325)xflt,yflt,zflt
      WRITE(10,2330)VR,niter,Rtime
 2300 FORMAT(//' Seismic Source Data:' 
     *      //'   Numbers of Sub_Fault 0 (NSTR,NDIP)=',2I5)
 2305 FORMAT( '   Strike Angle(rd)=',F12.8,',  Dip Angle   (rd)=',F12.8)
 2315 FORMAT( '   Fault Length (m)=',F12.3,',  Width(m)        =',F12.3)
 2320 FORMAT( '   Sub-Length   (m)=',F12.3,',  Sub-Width    (m)=',F12.3)
 2322 FORMAT( '   Str_Hypo.    (m)=',F12.3,',  Dip_Hypo     (m)=',F12.3)
 2324 FORMAT( '   X_Hypocenter (m)=',F12.3,',  Y_Hypocenter (m)=',F12.3
     *       /'   Z_Hypocenter (m)=',F12.3)
 2325 FORMAT( '   X Fault Edge (m)=',F12.3,',  Y Fault Edge (m)=',F12.3
     *       /'   Z Fault Edge (m)=',F12.3)
 2330 FORMAT( '   Vr         (m/s)=',F12.3,
     *       /'   Num of Iteration=',I12  ,',  Rise Time  (sec)=',F12.3)
c      
      DO 300 I=1,NDIP
      DO 300 J=1,NSTR
      WRITE(10,2350)I,J
      write(10,2355)FD(i,j),TD(i,j)
      write(10,2360)XF(i,j),YF(i,j),ZF(i,j)
  300 continue
 2350 FORMAT(/'   Sub_Fault 0 Number   = (',I2,',',I2,')')
 2355 format( '   Dislocation (m) and Delay Time (s) =',2F10.2) 
 2360 format( '   Locations of Center: X, Y, Z   (m)=',3f10.2)
C
C *** CHANGE OF SIGNS FOR IMAGINARY PARTS ***
C
      READ(7,1700)ISIGN
      IF(ISIGN.EQ.1)WRITE(10,2700)
      IF(ISIGN.EQ.0)WRITE(10,2710)
 1700 FORMAT(/I10)
 2700 FORMAT(/' SIGNS OF IMAGINARY PARTS FOR RESULTS ARE CHANGED')
 2710 FORMAT(/' SIGNS OF IMAGINARY PARTS FOR RESULTS ARE NOT CHANGED')
C
C *** OP DATA ***
C
      READ(7,1500)NOB
      WRITE(10,2500)NOB
 1500 FORMAT(/I10)
 2500 FORMAT(//' TOTAL NUMBER OF OBSERVATION POINTS =',I6/)
C
      IF(NOB.GT.NMOB)THEN
        WRITE(10,*)' NOB>NMOB !!!!!!!!!!!!!'
        STOP
      END IF
C
      DO 400 I=1,NOB
        READ(7,1520)XO(I),YO(I),ZO(I)
        WRITE(10,2520)I,XO(I),YO(I),ZO(I)
  400 CONTINUE
 1520 FORMAT(3F10.0)
 2520 FORMAT('   OB NUM.=',I5,',     (X,Y,Z) =',3F12.2)
C
      CLOSE(UNIT=7)
C
c *** output of fault and observation data ***
c
      open(9,file='fault5_dat',status='unknown')
        write(9,9500)XHYPO,YHYPO,ZHYPO
        write(9,9500)xflt,yflt,zflt
        write(9,9500)STR,DIP,0.0
        write(9,9500)slength,dwidth,VR
        write(9,9510)NSTR,NDIP,0
        write(9,9510)1,1,0
        do 700 i=1,ndip
        do 700 j=1,nstr
        write(9,9500)xf(i,j),yf(i,j),ZF(i,j)
  700   continue
 9500   format(3f12.3)
 9510   format(3i12)
      close(9) 
c
      open(9,file='obs_dat',status='unknown')
        do i=1,nob
        write(9,9500)xo(i),yo(i),zo(i)
        end do
      close(9) 
c
      open(9,file='time5_dat',status='unknown')
        write(9,9600)SLENGTH,DWIDTH,0
        write(9,9610)NSTR,NDIP,0
        write(9,9610)1,1,0
        do 800 i=1,ndip
	do 800 j=1,nstr
          da=fw*(i-1)+fw/2
          sa=fl*(j-1)+fl/2
          write(9,9620)sa,da,TD(i,j)
  800   continue
      close(9) 
 9600 format(2f10.1,i5)
 9610 format(3i5)
 9620 format(2f10.1,f10.4)

      RETURN
      END
C
C
C *** DATA OUTPUT ***
C
      SUBROUTINE OTDAT(MOM,MOB,GX)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 GX(MOM,MOB)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(13),YO(13),ZO(13),NOB
     *    /STRC/DNS,VP,VS,Qp,Qs,Qpn,Qsn
     *    /FLT/STR,DIP,FL,FW,Vr,Fmax,Rtime
     *    /IFLT/NSTR,NDIP,idum,niter
     *    /Hypo/xhypo,yhypo,zhypo 
     *    /SIGN/ISIGN
C
C *** CHECK ***
C
      NC=13
      IF(NOB.GT.NC)THEN
         WRITE(10,*)'FROM OTDAT:NOB>NC',NOB,NC
         STOP
      END IF
C
C *** OUTPUT OF DATA (GX, GY AND GZ)***
C
      IF(ISIGN.EQ.1)THEN
        DO 100 IOB=1,NOB
        DO 110 K=IOM,NOM
          GX(K,IOB)=DCMPLX(DREAL(GX(K,IOB)),-DIMAG(GX(K,IOB)))
  110   CONTINUE
  100   CONTINUE
      END IF
C      
      OPEN(UNIT=8,FILE='grfault.d',STATUS='UNKNOWN')
        WRITE(8,2500)DOM,IOM,NOM
        WRITE(8,2510)NOB
        DO 500 IOB=1,NOB
          DO 510 K=IOM,NOM
            WRITE(8,2520)GX(K,IOB)
  510     CONTINUE
  500   CONTINUE
      CLOSE(8)
C
 2500 FORMAT(1X,F10.6,2I10)
 2510 FORMAT(1X,I10,///)
 2520 FORMAT(1X,6E12.5)
C
      RETURN
      END
C
C
C *** COMPUTATION OF GREEN FUNCTIONS of the Homogeneous Full-Space ***
C
      SUBROUTINE BRSAH(MFLT,MOM,MOB,XF,YF,ZF,FD,TD,GX,FX,XR)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 GX(MOM,MOB),FX(MOM),FX1(MOM),EDT,cmp,F11
      DIMENSION FD(MFLT,MFLT),TD(MFLT,MFLT),XR(MOM)
     *         ,XF(MFLT,MFLT),YF(MFLT,MFLT),ZF(MFLT,MFLT)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(13),YO(13),ZO(13),NOB
     *    /STRC/DNS,VP,VS,Qp,Qs,Qpn,Qsn
     *    /FLT/STR,DIP,FL,FW,Vr,Fmax,Rtime
     *    /IFLT/NSTR,NDIP,idum,niter
     *    /SIGN/ISIGN
C
C *** CHECK OF DIMENSION ***
C
      NC=13
      IF(NOB.GE.NC)THEN
        WRITE(10,*)' FROM BRQVW:NOB>NC !!!!!!!!!!'
        WRITE(10,*)' NOB,NC=',NOB,NC
        STOP
      END IF
C
      NC=MOB
      IF(NOB.GE.NC)THEN
        WRITE(10,*)' FROM BRQVW:NOB>MOB !!!!!!!!!!'
        WRITE(10,*)' NOB,NC=',NOB,MOB
        STOP
      END IF
C
C   * Find closest (but larger) 2**J value to the input data dimension *
C
c      nom1=nom+1
      nom1=nom
      DO J = 1,999
        N = 2**J
        IF(NOM1.LE.N)THEN
             NUM = N
             GO TO 10
        END IF
      END DO
C
   10 IF(NUM*2.GT.MOM)THEN
        WRITE(*,*)
        WRITE(*,*)' NUM*2.GT.MOM !!!!!'
        WRITE(*,*)' Increase MOM in the parameter statement!'
        STOP
      ENDIF
c
      pi=datan(1.d0)*4.d0
      DUR = pi*2.d0/DOM
      NFAST = NUM*2
      DT = DUR/dFLOAT(NFAST)
      write(*,*)
      write(*,*)'dur,nfast,dt',dur,nfast,dt
      write(*,*)'New NOM for FFT, and Final F',num,dom*num/(pi*2)
      write(*,*)
c
c      Tdu=dsqrt(FL*FL+FW*FW)/Vr
c      Tw=2.d0*Tdu
c      Afac=dur/Tw
c      write(*,*)'Tw,Afac ',Tw,Afac
C
c   * Source Spectra *
c
      aND=dsqrt(dfloat(ndip*nstr))
      Rpatn=0.63d0
      PRTITN=1.d0/dsqrt(2.d0)
      write(*,*)'Rpatn,PRTITN =',Rpatn,PRTITN
c
      cmp=DCMPLX(0.D0,1.d0)
      DO 200 IOB=1,NOB

	do i=1,mom
          GX(i,iob)=(0.d0,0.d0)
	end do
c
        do 200 is=1,ndip
        do 200 js=1,nstr

          RX=XF(is,js)-XO(IOB)
          RY=YF(is,js)-YO(IOB)
          RZ=ZF(is,js)-ZO(IOB)
          RR=DSQRT(RX*RX+RY*RY+RZ*RZ)
c          write(*,*)'XF,YF,ZF',XF(is,js),YF(is,js),ZF(is,js)
c          write(*,*)'XO,YO,ZO',XO(iob),YO(iob),ZO(iob)
          slip=FD(is,js)/aND
c
          SMo=Dns*Vs*Vs*slip*FL*FW*10000000000.0d0
          aMag=(dlog10(SMo)-17.0)/1.33
          Tw=10.0**(0.31*aMag-0.77)
          fc=10.0**((23.8-dlog10(SMo))/3.0)
          fmax=7310.0*SMo**(-0.12)
          write(*,*)
          write(*,*)'iob,is,js = ',iob,is,js
          write(*,*)'Mo,M   = ',SMo,aMag
          write(*,*)'fc,fmax= ',fc,fmax
          write(*,*)'Tw,RR  = ',Tw,RR
          write(10,*)
          write(10,*)'iob,is,js = ',iob,is,js
          write(10,*)'Mo,M   = ',SMo,aMag
          write(10,*)'fc,fmax= ',fc,fmax
          write(10,*)'Tw,RR  = ',Tw,RR
c
          FX(1)=(0.d0,0.d0)
	  CM=Rpatn*PRTITN/(4.d0*pi*dns*Vs*Vs*Vs*10000000000.0d0)
          DO 100 IO=IOM,Num      
            OM=DOM*DFLOAT(IO)
            f=om/(2.d0*pi)
	    FFc=f/Fc
	    S=om*om/(1.d0+FFc*FFc)
            FFmax=f/Fmax
c   * by Boore (1983) *
c            P=1.d0/dsqrt(1.d0+FFmax*FFmax)
c   * by Kamae et al. (1991) *
            P=1.d0/(1.d0+FFmax)
              Q=Qs*f**Qsn
            E=dexp(-pi*f*RR/(Vs*Q))
            SA=CM*SMo*S*P*E/RR
            ran=ran0(idum)*pi*2.d0
            dr=SA*dcos(ran)
            di=SA*dsin(ran)
c   * by Irikura (1986) *
c          box=om*Rtime/2.d0
c          EDT=CDEXP(DCMPLX(0.D0,box))
c          F11=1.d0+(aND-1.d0)*dsin(box)*EDT/box
c   * by Irikura et al, (1996) *
          EDT=CDEXP(-1.d0-cmp*om*Rtime)
          aa=1.d0-dexp(-1.d0)
          F11=1.d0+(aND-1.d0)*(1.d0-EDT)/(1.d0+cmp*om*Rtime)/aa
c
          FX(io+1)=dcmplx(dr,di)*F11
          write(1,9999)om/(2*pi),cdabs(F11),cdabs(fx(io+1))*200
          FX1(io+1)=FX(io+1)
c
  100     CONTINUE
 9999 format(f12.4,4f12.5)
c
c   * Original waveform considering Boore's waveform functions *
c
          DO I=2,NFAST/2
            FX(NFAST+2-I)=dCMPLX(dREAL(FX(I)),-DIMAG(FX(I)))
          end do
          CALL dWVFAST(NFAST,FX,MOM,1)
          DO I = 1,NFAST
            XR(I)=dREAL(FX(I))/DUR
          end do
c
c   * by Boore (1983) *
c	  Tdu=1.d0/Fc
c	  Tw=2.d0*Tdu
c	  Afac=dur/Tw
c       write(*,*)'Tw,Afac = ',Tw,Afac
c
c   * by Kamae et al. (1991) *
c          aMag=(dlog10(SMo)-17.0)/1.33
c          Tw=10.0**(0.31*aMag-0.77)
c        write(*,*)'Mo,M,Tw= ',SMo,aMag,Tw
	  Afac=dur/Tw
c
	  b=1.2531d0
          c=6.266/Tw
          a=(2.718281829d0/(0.2*Tw))**b
          do i=1,nfast
	    time=dt*(i-1)
            XR(i)=XR(i)*Afac*a*(time**b)*dexp(-c*time)
c       write(2,9999)time,xr(i)
	  end do
c
c  * Iterations to obtain Boore's waveform functions *
c
       do 150 iter=1,Niter

	  do i=1,nfast
            FX(I)=dCMPLX(XR(I),0.0d0)
	  end do
          CALL dWVFAST(nfast,FX,MOM,-1)
          FX(1)=FX(1)*dt
          DO K=2,num+1
            FX(K)=FX(K)*dt
c       write(3,9999)dom*(k-1),fx(k),cdabs(fx(k))
            FX(K)=FX(K)*cdabs(FX1(K))/cdabs(FX(K))
c       write(4,9999)dom*(k-1),fx(k),cdabs(fx(k))
	  end do
c
          DO I=2,NFAST/2
            FX(NFAST+2-I)=dCMPLX(dREAL(FX(I)),-DIMAG(FX(I)))
          end do
          CALL dWVFAST(NFAST,FX,MOM,1)
          DO I = 1,NFAST
            XR(I)=dREAL(FX(I))/DUR
	    time=dt*(i-1)
c       write(7,9999)time,xr(i)
            if(time.lt.0.1)then
              XR(i)=XR(i)*time*10.0
            else
              if(time.gt.Tw)then
                 if(time.lt.2*tw)then
                   XR(i)=XR(i)*(2.0-time/Tw)
                 else
                   XR(i)=0.0
                 end if
              end if
            end if
c       write(8,9999)time,xr(i)*200.0
	  end do
c
  150   CONTINUE
c
c * Fourier Transform of Boore's waveform functions *
c
	  do i=1,nfast
            FX(I)=dCMPLX(XR(I),0.0d0)
	  end do
          CALL dWVFAST(nfast,FX,MOM,-1)
          FX(1)=(0.0,0.0)
          DO K=1,num+1
            om=dom*(k-1)
            FX(K)=FX(K)*dt
            FX1(K)=FX(K)*CDEXP(DCMPLX(0.D0,-OM*(TD(is,js)+RR/Vs)))
c       write(14,9999)om/(2*pi),cdabs(fx(k))*200
	  end do
c
c          DO I=2,NFAST/2
c            FX1(NFAST+2-I)=dCMPLX(dREAL(FX1(I)),-DIMAG(FX1(I)))
c          end do
c          CALL dWVFAST(NFAST,FX1,MOM,1)
c          DO I = 1,NFAST
c            XR(I)=dREAL(FX1(I))/DUR
c	    time=dt*(i-1)
c       write(15,9999)time,xr(i)*200
c	  end do
c
          DO I=1,num
	    om=dom*i
            EDT=CDEXP(DCMPLX(0.D0,OM*(TD(is,js)+RR/Vs)))
            if(ISIGN.eq.1)FX(I+1)=dCMPLX(dREAL(FX(I+1)),-DIMAG(FX(I+1)))
            GX(i,iob)=GX(i,iob)-FX(i+1)*EDT/(om*om)
          end do
C
  200 CONTINUE
c
      RETURN
      END
c
c *** Randam Function **
c
      function ran0(idum)
      IMPLICIT REAL*8(A-H,O-Z)
      parameter(IA=16807,IM=2147483647,AM=1.d0/IM,IQ=127773
     *         ,IR=2836,MASK=123459876)
c
      idum=ieor(idum,MASK)
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if(idum.lt.0)idum=idum+IM
      ran0=AM*idum
      idum=ieor(idum,MASK)
c
      return
      end
c
C
C *** FOURIER TRANSFORM OR INVERSE FOURIER TRANSFORM
C
      SUBROUTINE dWVFAST(N,A,MOM,IND)
c
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 A(MOM),TEMP,THETA
C
C * FFT and inverse FFT *
C
      PAI=3.1415926540d0
      J=1
      DO 140 I=1,N
        IF(I.GE.J) GO TO 110
        TEMP=A(J)
        A(J)=A(I)
        A(I)=TEMP
  110   M=N/2
  120   IF(J.LE.M) GO TO 130
        J=J-M
        M=M/2
        IF(M.GE.2) GO TO 120
  130   J=J+M
  140 CONTINUE
      KMAX=1
  150 IF(KMAX.GE.N) RETURN
      ISTEP=KMAX*2
      DO 170 K=1,KMAX
        THETA=dCMPLX(0.d0,PAI*dFLOAT(IND*(K-1))/dFLOAT(KMAX))
        DO 160 I=K,N,ISTEP
          J=I+KMAX
          TEMP=A(J)*cdEXP(THETA)
          A(J)=A(I)-TEMP
          A(I)=A(I)+TEMP
  160   CONTINUE
  170 CONTINUE
      KMAX=ISTEP
      GO TO 150
C
      END
