c
      PROGRAM RandACC1
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MOM=8200,MOB=13,MFLT=7)
c
      COMPLEX*16 FX(MOM),FX0,EDT
      DIMENSION XR(MOM)
c
      write(*,*)'Enter rise time'
      read(*,*)rtime
      write(*,*)'Enter n-dash'
      read(*,*)ndash
      SA=1.0

      dt=0.005
      nt=4096
      period=0.01
      duration=dt*nt
      dom=datan(1.d0)*8.d0/duration
      nom=duration/period
      write(*,*)
      write(*,*)'duration,dom,nom',duration,dom,nom
c
      num=nom
      nom1=nom+1
      pi=datan(1.d0)*4.d0
      DUR = pi*2.d0/DOM
      NFAST = NUM*2
      write(*,*)
      write(*,*)'dur,nfast,dt',dur,nfast,dt
c
      do io=1,nfast
        fx(io)=(0.0,0.0)
      end do

c      dash=dfloat(ndash)
c      FX(1)=dcmplx(dash,0.d0)
      FX(1)=dcmplx(1.d0,0.d0)
      DO 100 IO=2,nom1      
        FX(IO)=(0.d0,0.d0)
        OM=DOM*DFLOAT(IO-1)
        f=om/(2.d0*pi)
        ran=ran0(idum)*pi*2.d0
        dr=SA*dcos(ran)
        di=SA*dsin(ran)
        FX0=dcmplx(dr,di)
c	FX(io)=FX0
	do k=1,ndash
          dtime=rtime*dfloat((k-1))/dfloat(ndash)
          EDT=CDEXP(DCMPLX(0.D0,OM*dtime))
	  FX(io)=FX(io)+FX0*EDT
	end do
	FX(io)=FX(io)/dfloat(ndash)
  100 CONTINUE
c
      do io=1,nom1
         OM=DOM*DFLOAT(IO-1)
         f=om/(2.d0*pi)
         write(3,9999)om,f,cdabs(fx(io))
      end do
 9999 format(2f12.6,f12.5)
c
c   * Waveform using the inverse FFT ***
C
      DO io = NOM1+1,NUM+1
        FX(io) = DCMPLX(0.0,0.0)
      end do
C
      DO I=2,NFAST/2
        FX(NFAST+2-I)=dCMPLX(dREAL(FX(I)),-DIMAG(FX(I)))
      end do
C
      CALL dWVFAST(NFAST,FX,MOM,1)
c
      DO I = 1,NFAST
        XR(I)=dREAL(FX(I))/DUR
      end do
c
      do i=1,nfast
        write(4,9998)dt*(i-1),xr(i),xr(i),xr(i)
      end do
 9998 format(f8.4,3f12.8)
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
