      program resp40

c     This program calculates the response spectrum 

      real f(100), t(62), sa(3), sv(3), sd(3), psv(3), psa(3)
      real acc(8200), maxpga(3)
      real w, pi, beta, dt
      character*80 fileout, filein     
      integer nzero, npts, nskip, ev, stanum, tb, te, npole
      integer mpad
      real mag, delta, depth, azimuth
      complex cx(8200)
      common cx
      common /tmp1/acc
      data t/0.040, 0.046, 0.051, 0.056, 0.062, 0.068, 0.075,
     1       0.083, 0.091, 0.100, 0.110, 0.121, 0.133, 0.147,
     2       0.162, 0.178, 0.196, 0.215, 0.237, 0.261, 0.287,
     3       0.316, 0.348, 0.383, 0.422, 0.464, 0.511, 0.562,
     4       0.619, 0.681, 0.750, 0.825, 0.909, 1.000, 1.10,
     5       1.21, 1.33, 1.47, 1.62, 1.78, 1.96, 2.15, 2.37,
     6       2.61, 2.87, 3.16, 3.48, 3.83, 4.22, 4.64, 5.11,
     7       5.62, 6.19, 6.81, 7.50, 8.25, 9.09, 10.00, 11.00,
     8       12.12, 13.3, 15.0/

c     SET PARAMETERS
      fc = 0.1
      npole = 4
      nmin = 0

c tb = taper beginning percent.
      tb = 3
c te = taper ending percent.
      te = 3

      pi = 3.1415926

c calculate at 62 diferent periods.
      nper = 62

      write (*,*) 'Enter the value of damping Beta.'
      read (*,*) beta

c      write (*,*) 'Enter the number of stations.'
c      read (*,*) nsta     
      nsta = 1

      write (*,*) 'Enter the event,magnitude,delta,depth and azimuth.'
      read (*,*) ev,mag,delta,depth,azimuth

      write (*,*) 'Enter the input data file name.'
      read (*,'(a80)') filein

      Write (*,*) 'Enter the time steps and npts for the data.'
      read (*,*) dt,npts

      write (*,*) 'Enter the plotxy output file name.'  
      read (*,'( a80)') fileout

      open (20,file=fileout,status='unknown')
      open (4,file=filein,status='old')

      write (20,'( 2x,f10.4,''  damping'')') beta
      write (20,'(i5,4f10.4)') ev,mag,delta,depth,azimuth

c calculate the 62 different freq.
      nfreq = nper
      do 2 i=1,nfreq
        f(i) = 1./T(i)
  2   continue

c Number of components.
      ncomp = 1

      do 600 jj=1,nsta

        do 300 icomp=1,ncomp
          nskip = 0

c         READ THE TIME HISTORY FILE
c          call rdat (acc,dt,npts)
          do 17 k=1,npts,1
             read (4,18) time,acc(k)
c            write (*,18) time,acc(k)
 17       continue
 18       format (f10.3,f10.3)

c         REMOVE THE DC
          call rdc (acc,npts,0)

c         TAPER
          call taper (acc,npts,tb,te)

c         PAD TO POWER OF 2
          call pad (acc,npts,nmin,mpad)
          df = 1./(npts*dt)
          write (*,'( 2x,''mpad ='',i4)') mpad

c         LOAD INTO COMPLEX ARRAY
          do 10 i=1,npts
            cx(i) = cmplx(acc(i),0.)
  10      continue

c         FFT
          write (*,'( 2x,''FFT'')')
          call cool (-1.,mpad)

c         FILTER
          write (*,'( 2x,''FILTER'')')
          call hpass(npole,npts,fc,df)

c         INVERSE FFT
          write (*,'( 2x,''INVERSE FFT'')')
          call cool (1.,mpad)

c         COPY TO REAL ARRAY
          do 15 i=1,npts
            acc(i) = real(cx(i))/npts
  15      continue

c         FIND MAX PGA
          maxpga(icomp) = 0.0
          do 20 i=1,npts
            maxpga(icomp) = amax1(maxpga(icomp),abs(acc(i)))
  20      continue

c         LOOP FOR EACH FREQUENCY
          do 200 l=1,nfreq
            w = 2.*pi*f(l)
            call coeff (w,beta,dt)

c           CALCULATE THE ACCELERATION RESPONSE 
            call brs (acc,w,beta,npts,sa1,sv1,sd1)

            SA(icomp) = sa1/980.
            SV(icomp) = sv1
            SD(icomp) = sd1
            PSV(icomp) = w*sd1
            PSA(icomp) = w**2*sd1
            write (*,5500) l,PSV(icomp),PSA(icomp),SA(icomp),SV(icomp),
     1                     SD(icomp)
 5500       format( 2x,i2,3x,5f9.3)

            write (20,'( 2x,7f9.3)') f(l),psv(icomp),psa(icomp),
     1         sa(icomp),sa(icomp)/maxpga(icomp),SV(icomp),SD(icomp)
 200      continue
 300    continue
 600  continue

      close (20)
      close (4)

      stop
      end

c************************************************************
c      subroutine rdat(acc,dt,npts)
c This subroutine will read the data file from the aux01 disk
c of the Loma Prieta Mainshock data. These records have been 
c processed by CSMIP.

c      real acc(8200)

c     Read the data

c     do 100 i=npts+1,1
c       read (4,1001) time,acc(i)
c  100 continue

c      return
c 1001 format(f10.3,f10.3)
c      end


c*********************************
      subroutine coeff (w,beta,dt)
      common /coef/a11,a12,a21,a22,b11,b12,b21,b22
c
c     Set up reapeated terms
      t1 = sqrt(1.-beta**2)
      t2 = sin (w*t1*dt)
      t3 = cos (w*t1*dt)
      t4 = exp (-beta*w*dt)
      s1 = (2.*beta**2-1.) / (w**2*dt)
      s2 = 2.*beta / (w**3*dt)
c
c     calculate the a's
      a11 = t4*(beta*t2/t1+t3)
      a12 = t4*t2 / (w*t1)
      a21 = -t4*w*t2 / t1
      a22 = t4*(t3-beta*t2/t1)
c
c     calculate the b's
      b11 = t4*((s1+beta/w)*t2 / (w*t1) + (s2+1./w**2)*t3) - s2
      b12 = -t4*(s1*t2/(w*t1)+s2*t3) - 1./w**2 + s2
      b21 = (s1+beta/w) * (t3-beta*t2/t1)
      b21 = t4*(b21 - (s2+1./w**2)*(w*t1*t2+beta*w*t3)) + 1./(w**2*dt)
      b22 = s1*(t3-beta*t2/t1)
      b22 = -t4*(b22 - s2*(w*t1*t2+beta*w*t3)) - 1./(w**2*dt)
      return
      end

c********************************************
      subroutine brs (x,w,beta,npts,sa,sv,sd)
      real x(1), sa, sv, sd, w, beta
      integer npts
      common /coef/ a11,a12,a21,a22,b11,b12,b21,b22
c
c     initialize
      t1 = 2.*beta*w
      t2 = w**2
      d = 0.
      v = 0.
      a = 0.
      Sa = 0.0
      Sv = 0.0
      Sd = 0.0
c
c     calculate the response
      do 10 i=1,npts
        ap1 = x(i)
        dp1 = a11*d + a12*v + b11*a + b12*ap1
        vp1 = a21*d + a22*v + b21*a + b22*ap1
        z = -(t1*vp1 + t2*dp1)
        Sa = amax1(Sa,abs(z))
        Sv = amax1(Sv,abs(vp1))
        Sd = amax1(Sd,abs(dp1))
        a = ap1
        v = vp1
        d = dp1
  10  continue
      return
      end 

c**********************************
      subroutine rdc (x,npts,iflag)

c     This subroutine removes a dc shift from the data

c     IFLAG = 0    remove the mean
c     IFLAG = 1    remove the mean value of the first 10 points
c     IFLAG = 2    manual set of DC value to be removed

      real x(1),sum,mean
      integer npts,iflag

      if (iflag .eq. 0) then
        sum = 0.0
        do 10 i=1,npts
          sum = x(i) + sum
  10    continue
        mean = sum/float(npts)

      elseif (iflag .eq. 1) then
        sum = 0.0
        do 20 i=1,10
          sum = x(i) + sum
  20    continue
        mean = sum / float(npts)

      else
        write (*,1000)
        read (*,1001) mean

      endif

      do 100 i=1,npts
        x(i) = x(i) - mean
  100 continue
      write (*,1010) mean
 1010 format( 2x,'Remove DC of ',f12.8)

      return
 1000 format( 2x,'Enter mean to be removed')
 1001 format( f12.8)
      end

c************************************
      subroutine taper (x,npts,tb,te)

c     This subroutine tapers the x array
    
      real x(1), arg
      integer npts,tb,te
      pio2 = 3.1415926/2.

      if (tb .ne. 0.) then
        n = (npts*tb)/100
        do 10 i=1,n
          arg = pio2*float(i-1)/float(n)
          x(i) = x(i)*sin(arg)
  10    continue
        write (*,1000) tb
 1000   format( 2x,'Taper beginning ',i2,' percent')

      endif

      if (te .ne. 0.) then
        n = (npts*te)/100
        do 20 i=1,n
          arg = pio2*float(i-1)/float(n)
          x(npts-i+1) = x(npts-i+1) * sin(arg)
  20    continue
        write (*,1001) te
 1001   format( 2x,'Taper end ',i2,' percent')

      endif
      return
      end

c***********************************
      subroutine pad (x,npts,nmin,m)

c     This subroutine pads the x array to a power of 2

      real x(1)
      integer test,npts,m,nmin

      do 10 i=1,20
        test = 2**i
        if (test .ge. npts .and. test .ge. nmin) goto 20
  10  continue
  20  m = i
      write (*,1000) test,m
 1000 format( 2x,'Pad to ',i4,' points,  m =',i4)

      do 30 i=npts,test
        x(i) = 0.0
  30  continue
      npts = test

      return
      end

c******************************
      subroutine cool (signi,n)
      complex cx(8200),carg,temp,cw
      common cx
      pi = 3.1415926*signi
      lx = 2**n
      j = 1
      do 30 i=1,lx
        if (i .gt. j) goto 10
        temp = cx(j)
        cx(j) = cx(i)
        cx(i) = temp
  10    m = lx/2
  20    if (j .le. m) goto 25
        j = j-m
        m = m/2
        if (m .ge. 1) goto 20
  25    j = j+m
  30  continue
      l = 1
  40  istep = l+l
      do 50 m=1,l
        carg = cmplx(0.,pi*float(m-1)/float(l))
        cw = cexp(carg)
        do 45 i=m,lx,istep
          temp = cw*cx(i+l) 
          cx(i+l) = cx(i)-temp
          cx(i) = cx(i)+temp
  45    continue
  50  continue
      l = istep
      if (l .lt. lx) goto 40
      return
      end

c***********************************
      subroutine hpass(n,npts,fc,df)
c     This subroutine is a "n" pole high-pass Butterworth filter
      complex cx(8200)
      real omega, df, fc
      integer npts
      common cx
      
      domega = df/fc
      do 100 i=2,npts/2
        omega = (i-1)*domega
        YH = sqrt(omega**(2*n)/(1.+omega**(2*n)))
        cx(i) = cx(i)*YH
        cx(npts-i+2) = cx(npts-i+2)*YH
 100  continue
      cx(1) = cmplx(0.,0.)
      return
      end
