c
c *** Simulated Acceleration based on Kokuji No.1461 (2000) ***
c
      PROGRAM kokuji
c
c  Parameters:
c   mm: Number of Delta Time for Main Dimension (> nn) 
c   nn: Actual Number of Delta Time
c   dt: Delta Time (dt*nn = duration)
c    h: Dampling Factor
c   sa: Maximum Response Acceleration (gal)
c   sv: Maximum Response Velocity     (cm/sec)
c   sd: Maximum Response Displacement (cm)
c
c  Dimension Data:
c   ddy: Ground Acceleration (gal)
c   Fds: Fourier Spectra of Ground Acceleration (gal)
c   acc: Response Acceleration (gal)
c   vel: Response Velocity     (cm/sec)
c   dis: Response Displacement (cm)
c   resa: Response Acceleration (gal)
C
      PARAMETER(mm=66000)
      dimension time(mm),ddy(mm),eamp(mm)
      dimension acc(mm),vel(mm),dis(mm),resa(mm,0:10)
      COMPLEX Fds(mm),Fds1(mm)
      CHARACTER INFILE*20
C
C  Prompt for input RECON'ed file
C
      WRITE(*,1010) 
 1010 FORMAT(1X,'Enter Acceleration Wave Data')
      READ(*,1020) INFILE
 1020 FORMAT(A20)
      write(*,*)
   10 write(*,*)'Enter Number of Iterations (Nit: up to 10)'
      read(*,*)nit
      if(nit.gt.10)goto 10
c
c      write(*,*)
c      write(*,*)'Level 1 or 2? (Enter 1 or 2)'
c      read(*,*)L12
c
      tb=4.0
      tc=35.0
      td=60.0
      te=80.0
      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
      ampfac=1.0/float(nit)
      write(*,*)'Amplitue Factor for the Head and Tail of Waveform:'
      write(*,*)'  Initial Value (ampfac=1/Nit) =',ampfac
      write(*,*)'Change Those the Value? (Enter No=0 or Yes=1)'
      read(*,*)Nampfac
      if(Nampfac.eq.1)then
        write(*,*)'Enter ampfac (If ampfac=1, No Amplification):'
        read(*,*)ampfac
      end if
c
c      write(*,*)
c      write(*,*)'Skip Envelope Function? (Enter No=0 or Yes=1)'
c      read(*,*)Nenv
C
C *** Read input data file  ***
C
      OPEN(11,FILE=INFILE)
        read(11,*)
        i=1
   20   read(11,*,end=21)time(i),ddy(i)
c        if(L12.eq.1)ddy(i)=ddy(i)/5.0
        i=i+1
        goto 20
   21 CLOSE(11)
c
      NT=i-1
      dt=time(2)-time(1)
      write(*,*)'NT and DT=',NT,DT
c
C *** Find the Closest Number of 2**N ***
C
      NF=2
      do I=2,15
        NF=2*NF
        IF(NF.GE.NT)GOTO 30
      end do
      WRITE(*,*)' NT>2**15 !!!!!'
      STOP
   30 IF1=I
      WRITE(*,*)
      WRITE(*,*)' IF1,NF=',IF1,NF
c
C *** FOURIER TRANSFORM of the Original Wave ***
C
      do I=1,NT
        Fds(I)=CMPLX(ddy(I),0.0)
      end do
      do I=NT+1,NF
        Fds(I)=(0.0,0.0)
      end do
C
      CALL WVFAST(NF,Fds,mm,-1)
      DUR=FLOAT(NF)*DT
      DOM=6.2831853/DUR
      NOM=NF/2+1
      EOM=DOM*NOM
      WRITE(*,*)
      WRITE(*,*)' DUR,DOM,NOM,EOM=',DUR,DOM,NOM,EOM
      WRITE(*,*)' Highest Frequency (hz)=',EOM/6.2831853
      do I=1,NF
        Fds(I)=Fds(I)*DT
      end do
c
c *** Response Spectrum of Kokuji No.1461 (2000) ***
c
      h=0.05
      OM1=6.283185/0.64
      OM2=6.2831853/0.16
      do j=2,nom
        OMEGA=FLOAT(j-1)*DOM
        period=6.2831853/omega
        resa(j,0)=800.0
        IF(OMEGA.LT.OM1)resa(j,0)=512.0/period
        IF(OMEGA.GT.OM2)resa(j,0)=period*3000.0+320.0
c        if(L12.eq.1)resa(j,0)=resa(j,0)/5.0
      end do
C
C *** Modification of Original Acceleration ***
c
      do I = 1,NF
         t=dt*(i-1)
         if(t.le.tb)then
           eamp(I)=t*t/(tb*tb)
         else
           if(t.gt.te)then
             eamp(i)=0.0
           else
             if(t.gt.td)then
               eamp(I)=0.1*((te-t)/(te-td))**2
             else
	       if(t.gt.tc)then
                 eamp(I)=exp(alog(0.1)*(t-tc)/(td-tc))
               else
                 eamp(i)=1.0
               end if
             end if
           end if
        end if
      end do
c
      write(*,*)
c
      do 100 iit=1,nit
      
        cr1=9999.0
        cr2=0.0
        cr2a=0.0
        cr2b=0.0
        cr3=0.0
        cr4=0.0
        
        write(*,*)' Iteration Number =',iit
        Fds(1)=(0.0,0.0)
        Fds1(1)=(0.0,0.0)
        do j=2,nom
          OMEGA=FLOAT(j-1)*DOM
          period=6.2831853/omega
          call resp(h,omega,dt,nt,ddy,acc,vel,dis,mm,sa,sv,sd)
          Fds(j)=Fds(j)*resa(j,0)/sa
          Fds1(j)=Fds(j)
        end do
C
        do I=2,NF/2
           Fds1(NF+2-I)=CMPLX(REAL(Fds1(I)),-IMAG(Fds1(I)))
        end do
c
        CALL WVFAST(NF,Fds1,mm,1)
        do I = 1,NF
           ddy(I)=REAL(Fds1(I))/DUR
           ddy(I)=ddy(I)*(eamp(i)**ampfac)
        end do
        
        ncr=0
        do j=2,nom
          OMEGA=FLOAT(j-1)*DOM
          period=6.2831853/omega
          call resp(h,omega,dt,nt,ddy,acc,vel,dis,mm,sa,sv,sd)
          resa(j,iit)=sa
          if((period.ge.0.02).and.(period.le.10.0))then
            ncr=ncr+1
            eps=sa/resa(j,0)
            if(eps.lt.cr1)cr1=eps
            cr3=cr3+(eps-1.0)**2
            cr4=cr4+eps
          end if
          if((period.ge.1.0).and.(period.le.5.0))then
            cr2a=cr2a+sa
            cr2b=cr2b+resa(j,0)
          end if
        end do
        cr2=cr2a/cr2b
        cr3=sqrt(cr3/float(ncr))
        cr4=cr4/float(ncr)
        write(*,*)'   Error Criteria 1 (>=0.85):',cr1
        write(*,*)'   Error Criteria 2 (>=1.0) :',cr2
        write(*,*)'   Error Criteria 3 (<=0.05):',cr3
        write(*,*)'   Error Criteria 4 (>=0.95):',cr4
c
  100 CONTINUE
C
c *** output ***
c
      write(*,*)
      write(*,*)'Enter Response Spectra Data'
      READ(*,*)INFILE
      OPEN(11,FILE=INFILE,STATUS='unknown')
      write(11,4567)(I,I=1,nit)
      period0=9999.9
      do 200 j=nom,2,-1
        OMEGA=FLOAT(j-1)*DOM
        period1=6.2831853/omega
        if(abs(period1-period0).lt.0.01)goto 200
        write(11,2000)period1,(resa(j,i),i=0,nit)
        period0=period1
  200 continue
      CLOSE(11)
 4567 format('period(sec),kokuji1461',20(',Iteration No.',I3))
 2000 format(f8.4,20(',',e12.4))
c
      WRITE(*,2010) 
 2010 FORMAT(/1X,'Enter output acceleration file name')
      READ(*,2020) INFILE
 2020 FORMAT(A20)
      OPEN(12,FILE=INFILE,STATUS='unknown')
      write(12,*)'time(sec),Acc(gal),eamp,eamp*ampfac'
      do i=1,nt
      write(12,2000)dt*(i-1),ddy(I),eamp(i),eamp(i)**ampfac
      end do
      CLOSE(12)
c
        STOP
	END
C
C *** FOURIER TRANSFORM OR INVERSE FOURIER TRANSFORM
C
      SUBROUTINE WVFAST(N,A,MOM,IND)
C
C *** N   : NUMBER OF FOURIER COEFF. OR WAVE ( = 2**I )
C *** A(I): FOURIER COEFFICIENTS ( N/2+1 ) OR WAVE ( REAL PART )
C    ( FOURIER SPECTRUM = FOURIER COEFFICIENTS * DURATION/2.0 )
C *** MOM : N IN MAIN DIMMENSION
C         : IF mm > 2**13(=8192), CAN NOT USE !!!
C *** IND : = -1 : TRANS.,   = 1 : INVERSE TRANS.
C
      COMPLEX A(MOM),TEMP,THETA
C
C *** CALUCURATION OF FFT
C
      PAI=3.1415927
      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=CMPLX(0.0,PAI*FLOAT(IND*(K-1))/FLOAT(KMAX))
        DO 160 I=K,N,ISTEP
          J=I+KMAX
          TEMP=A(J)*CEXP(THETA)
          A(J)=A(I)-TEMP
          A(I)=A(I)+TEMP
  160   CONTINUE
  170 CONTINUE
      KMAX=ISTEP
      GO TO 150
C
      END
c
c
c	Subroutine for Response of Single-Degree of System
c
	subroutine resp(h,w,dt,nn,ddy,acc,vel,dis,nd,sa,sv,sd)
c
	dimension ddy(nd),acc(nd),vel(nd),dis(nd)
c
	w2=w*w
	hw=h*w
	wd=w*sqrt(1.-h*h)
	wdt=wd*dt
	e=exp(-hw*dt)
	cwdt=cos(wdt)
	swdt=sin(wdt)
	a11= e*(cwdt+hW*swdt/wd)
	a12= e*swdt/wd
	a21=-e*w2*swdt/wd
	a22= e*(cwdt-hw*swdt/wd)
	ss=-hw*swdt-wd*cwdt
	cc=-hw*cwdt+wd*swdt
	s1=(e*ss+wd)/w2
	c1=(e*cc+hw)/w2
	s2=(e*dt*ss+hw*s1+wd*c1)/w2
	c2=(e*dt*cc+hw*c1-wd*s1)/w2
	s3=dt*s1-s2
	c3=dt*c1-c2
	b11=-s2/wdt
	b12=-s3/wdt
	b21=(hw*s2-wd*c2)/wdt
	b22=(hw*s3-wd*c3)/wdt
	acc(1)=2.*h*w*ddy(1)*dt
	vel(1)=-ddy(1)*dt
	dis(1)=0.
	dx=vel(1)
	x=0
	sa=0
	sv=0	
	sd=0
c
	do 110 m=2,nn
	dxf=dx
	xf=x
	ddym=ddy(m)
	ddyf=ddy(m-1)
	x = a12*dxf+a11*xf+b12*ddym+b11*ddyf
	dx= a22*dxf+a21*xf+b22*ddym+b21*ddyf
	ddx=-2.*hw*dx-w2*x
	acc(m)=ddx
	vel(m)=dx
	dis(m)=x
	sa=amax1(sa,abs(ddx))
	sv=amax1(sv,abs(dx))
	sd=amax1(sd,abs(x))
  110   continue
	return
	end
