c
c *** Simulated Acceleration based on Jennings's method ***
c
      PROGRAM jennings
c
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(mm=131080)
      dimension acc(mm)
      COMPLEX*16 FX(mm)
      character INFILE*50
C

C  Prompt for input RECON'ed file
C
      write(*,*)
      write(*,*)'Enter Random Number (1-999)'
      read(*,*)idum
c
      pi=datan(1.d0)*4.d0
      dt=0.01d0
      nt=2**16
      DUR=dt*nt
      DOM=pi*2.d0/DUR
      EOM=pi/dt
      NOM=nt/2
      WRITE(*,*)
      WRITE(*,*)' dt,nt,duration(=dt*nt)=',dt,nt,dur
      WRITE(*,*)' DUR,DOM,NOM,EOM=',DUR,DOM,NOM,EOM
      WRITE(*,*)' Highest Frequency (hz)=',EOM/(2.0d0*pi)
c
c *** 2000 Kokuji Spectrum ***
c
      per1=0.16d0
      per2=0.64d0
      fx(1)=(0.d0,0.d0)
      DO 100 IO=1,NOM
        OM=DOM*DFLOAT(IO)
        period=2.0d0*pi/OM
        ran=ran0(idum)*pi*2.d0
        dr=dcos(ran)
        di=dsin(ran)
        Sa=800.d0
        IF(period.lt.per1)Sa=(320.d0+3000.d0*period)
        IF(period.ge.per2)Sa=512.d0/period
        pSv=Sa*period/(pi*pi)
        FX(io+1)=dcmplx(dr,di)*pSv
c        write(3,9999)period,fx(io+1),cdabs(fx(io+1))
c     *   ,cdabs(fx(io+1))*OM
  100 CONTINUE
c 9999 format(10e12.4)
c
      NFAST=nt
      DO I=2,NFAST/2
        FX(NFAST+2-I)=dCMPLX(dREAL(FX(I)),-DIMAG(FX(I)))
      end do
      CALL dWVFAST(NFAST,FX,MOM,1)
c
      tb=4.0d0
      tc=35.0d0
      td=80.0d0
      te=dur
      Afac=0.1d0
c
      write(*,*)
      write(*,*)'Corner Times for Waveform:'
      write(*,*)'  tb,tc,td,te(s)=',tb,tc,td,te
      write(*,*)'Change Those Data? (Enter No=0 or Yes=1)'
      read(*,*)Ncorner
      if(Ncorner.eq.1)then
        write(*,*)'Enter tb,tc,td,te(s):'
        read(*,*)tb,tc,td,te
      end if
c
      write(*,*)
      write(*,*)'Enter a Seed Acceleration for OUTPUT'
      read(*,*)INFILE
      open(9,file=INFILE,status='unknown')
      write(9,*)'time(sec),seed-Acc(gal)'
c
      DO I = 1,NFAST
        time=dt*(i-1)
        acc(I)=dREAL(FX(I))/DUR
        if(time.lt.tb)then
          AA=time*time/(tb*tb)
        else
          if(time.lt.tc)then
            AA=1.0d0
          else
            if(time.lt.td)then
              AA=dexp(dlog(Afac)*(time-tc)/(td-tc))
            else
              if(time.lt.te)then
                AA=Afac*((te-time)/(te-td))**2
              else
                AA=0.d0
              end if
            end if
          end if
        end if
        acc(I)=acc(I)*AA
        write(9,9999)time,acc(I)
 9999 format(f8.2,',',e12.4)
      end do
      close(9)
c
      stop
      end
c
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
