C
c *** Calculation of Shindo (Coded by Y. Hisada, 1997) ***
c *** If you have any questions, please ask to hisada@cc.kogakuin.ac.jp ***
c
      PROGRAM shindo1
C
      PARAMETER(MWV=33000)
      REAL TIME(MWV),XR1(MWV),XR2(MWV),XR3(MWV),xx(mwv)
      COMPLEX*8 FX1(MWV),FX2(MWV),FX3(MWV)
      character infile*20
C
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) 'Calculation of Shindo (Coded by Y. Hisada, 1997)'
      WRITE(*,*) '------------------------------------------------'
      WRITE(*,*)
C      
C *** Read input data file  ***
C
      write(*,*)'Enter acceleration file name:'
      read(*,*)infile
      open(11,file=infile,status='old')
      i=0
   10 i=i+1
      read(11,*,end=20)time(i),xr1(i),xr2(i),xr3(i)
      goto 10
   20 close(11)
      nt=i-1
      dt=time(2)-time(1)
      write(*,*)' nt, dt=',nt,dt
C
C *** FOURIER TRANSFORM ***
C
      NF=2
      DO 200 I=2,15
        NF=2*NF
        IF(NF.GE.NT)GOTO 30
  200 CONTINUE
      WRITE(*,*)' NT>2**15 !!!!!'
      STOP
   30 IF1=I
      WRITE(*,*)
      WRITE(*,*)' IF1,NF=',IF1,NF
C
      DO 300 I=1,NT
        FX1(I)=CMPLX(XR1(I),0.0)
        FX2(I)=CMPLX(XR2(I),0.0)
        FX3(I)=CMPLX(XR3(I),0.0)
  300 CONTINUE

      DO 310 I=NT+1,NF
        FX1(I)=(0.0,0.0)
        FX2(I)=(0.0,0.0)
        FX3(I)=(0.0,0.0)
  310 CONTINUE
C
      CALL WVFAST(NF,FX1,MWV,-1)
      CALL WVFAST(NF,FX2,MWV,-1)
      CALL WVFAST(NF,FX3,MWV,-1)
      DUR=FLOAT(NF)*DT

      DO 400 I=1,NF
        FX1(I)=FX1(I)*DT
        FX2(I)=FX2(I)*DT
        FX3(I)=FX3(I)*DT
  400 CONTINUE
C
c      FOM=0.0
      DOM=6.2831853070e0/DUR
      NOM=NF/2+1
      EOM=DOM*NOM
      WRITE(*,*)
      WRITE(*,*)' DUR,DOM,NOM,EOM=',DUR,DOM,NOM,EOM
C     
C *** BAND PASS filter ***
C
      do i=1,nom

	freq=float(i)/dur
        f1=sqrt(1.0/freq)

        y=freq/10.0
        y2=y*y
        y4=y2*y2
        y6=y2*y4
        y8=y4*y4
        y10=y4*y6
        y12=y4*y8
        f2=1.0/sqrt(1+0.694*y2+0.241*y4+0.0557*y6+0.009664*y8+
     *               0.00134*y10+0.000155*y12)  

	y=2.0*freq
        y3=y*y*y
        f3=sqrt(1.0-exp(-y3))

        aa=f1*f2*f3
	FX1(I) = AA*FX1(I)
	FX2(I) = AA*FX2(I)
	FX3(I) = AA*FX3(I)

      end do
C
C *** INVERSE FOURIER TRANSFORM  ***
C
      DO 600 I=2,NF/2
           FX1(NF+2-I)=CMPLX(REAL(FX1(I)),-AIMAG(FX1(I)))
           FX2(NF+2-I)=CMPLX(REAL(FX2(I)),-AIMAG(FX2(I)))
           FX3(NF+2-I)=CMPLX(REAL(FX3(I)),-AIMAG(FX3(I)))
  600 CONTINUE
C
      CALL WVFAST(NF,FX1,MWV,1)
      CALL WVFAST(NF,FX2,MWV,1)
      CALL WVFAST(NF,FX3,MWV,1)
C
C *** filtered WAVE DATA ***
C
      write(*,*)
      WRITE(*,*)'wave1_dat (filtered wave):' 
      OPEN(11,FILE='wave1_dat',STATUS='unknown')

      amax=0.0
      DO 700 I = 1,Nt
         TIME(I) = FLOAT(I-1)*DT
         X1=REAL(FX1(I))/DUR
         X2=REAL(FX2(I))/DUR
         X3=REAL(FX3(I))/DUR
         xx(i)=sqrt(x1*x1+x2*x2+x3*x3)
         if(xx(i).gt.amax)amax=xx(i)
         write(11,2000)dt*(i-1),XX(I)
  700 CONTINUE
 2000 format(4e12.4)
      CLOSE(11)

      write(*,*)
      write(*,*)' Amax of Filtered wave (gal) =',amax
c
c *** Calculation of the A0-Tau	relation ***
c
      write(*,*)
      WRITE(*,*)'a0tau_dat (A0-Tau Relation):' 
      OPEN(11,FILE='a0tau_dat',STATUS='unknown')

      na=1000
      da=amax/na
      jf=0
      do 800 j=1,na
	aa=amax-da*j
        n=0
	do 820 i=1,nt
	  if(xx(i).gt.aa)n=n+1
  820	continue
        write(11,2000)aa,n*dt
        if((n*dt.gt.0.299).and.(jf.eq.0))then
          a0=aa
          jf=1
        end if
  800 continue
      close(11)

      write(*,*)
      write(*,*)' A0 (gal) =',a0
      write(10,*)' A0 (gal) =',a0
c
c *** Keisoku Shindo ***
c
      shindo=2.0*log10(a0)+0.94
      write(*,*)' Shindo   =',shindo
      write(10,*)' Shindo   =',shindo
      write(*,*)
      write(*,*)'Data are stored in fort.10 and wave1_dat'
      read(*,*)dummy
 
      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 MWV > 2**13(=8192), CAN NOT USE !!!
C *** IND : = -1 : TRANS.,   = 1 : INVERSE TRANS.
C
      COMPLEX*8 A(MOM),TEMP,THETA
C
C *** CALUCURATION OF FFT
C
      PAI=3.1415926540e0
      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
