c
c *** Simulated Acceleration based on Jennings's method ***
c
      PROGRAM jennings
c
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(mm=8200)
      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.01
      nt=8192
      DUR=dt*nt
      DOM=pi*2.d0/DUR
      EOM=pi/dt
      NOM=EOM/DOM
      WRITE(*,*)
      WRITE(*,*)' DUR,DOM,NOM,EOM=',DUR,DOM,NOM,EOM
      WRITE(*,*)' Highest Frequency (hz)=',EOM/(2.0*pi)
c
      per1=0.64
c      per2=0.16
      fx(1)=(0.d0,0.d0)
      DO 100 IO=1,NOM
        OM=DOM*DFLOAT(IO)
        period=2.0*pi/OM
        ran=ran0(idum)*pi*2.d0
        dr=dcos(ran)
        di=dsin(ran)
        AA=81.5
        IF(period.lt.per1)AA=127.3*period
c        IF(period.lt.per2)AA=(50.9+477.*period)*period
        FX(io+1)=dcmplx(dr,di)*AA
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
      write(*,*)
      write(*,*)'Enter a Seed Acceleration'
      read(*,*)INFILE
      open(9,file=INFILE,status='unknown')
      write(9,*)'time(sec),seed-Acc(gal)'
      t1=4.0
      t2=35.0
c      t3=80.0
      DO I = 1,NFAST
        time=dt*(i-1)
        acc(I)=dREAL(FX(I))/DUR
        if(time.lt.t1)then
          AA=time*time/16.0
        else
          if(time.lt.t2)then
            AA=1.0
          else
c            if(time.lt.t3)then
c              AA=dexp(-0.0357*(time-35.0))
              AA=dexp(-0.07*(time-35.0))
c            else
c              AA=0.05+0.0000938*(120.0-time)*(120.0-time)
c            end if
          end if
        end if
        AA=AA*5.0
c        write(4,9999)time,AA,acc(i),acc(i)*AA
        write(9,9999)time,acc(i)*AA
 9999 format(f8.4(',',e12.4))
      end do
      close(9)

      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
