C     Program to integrate equation of motion subjected to 
C     acceleration input. Equation is integrated for a given range of 
c     frequencies and damping parameter.  Maximum response will define
C     response spectrum which can be used in modal dynamics option to
C     estimate maximum response of the considered structure. 
C
      subroutine abqmain
C    
C     program respon
      implicit real*8(a-h,o-z)
      parameter(nacc=1001, freqmin=0.1d0, freqmax=40.d0, intxxx=501,
     *     zero=0.d0, one=1.d0, two=2.d0)
      dimension a(2,2),b(2,2),q(nacc),qp(nacc),qpp(nacc),acc(2500),
     *     tab(8),fr(intxxx),qm(intxxx),qpm(intxxx),qppm(intxxx)
C     
C     Storage allocated for a maximum of 2500 timepoints in the 
C     earthquake acceleration history read from file QUAKE.AMP
C     (copied from CANTILEVER_QUAKEDATA.INP). Acceleration spectrum
C     will be written to file SPECTRUM.ACC. 
C     
C     Time integration for linear acceleration (exact solution)
C     Data damp denotes damping as percentage of critical damping
C     Data freqmin and freqmax define frequency range
C     Data intxxx defines number of points in frequency range
C     
      open(unit=101,status='old',file='QUAKE.AMP')
      open(unit=117,status='new',file='SPECTRUM.ACC')
C
      dt=0.01d0
      fact=-386.09d0
      frac=one/dble(intxxx-1)
      damp=0.02d0
C***  initiate amat,bmat before time integration
      do i=1,2
         do j=1,2
            a(i,j)=zero
            b(i,j)=zero
         end do
      end do
C***  read amplitude data and store on acc(2500)
      do i=1,625
         read(101,*) (tab(ll),ll=1,8)
         do k=1,4
            acc((i-1)*4+k)=tab(2*k)*fact
         end do
      end do
      dfreq=freqmax-freqmin
C**   choose damping
      if(damp.gt.one) write(6,11)
 11   format(/,3x,
     *     'This program is written for underdamped cases only')
C**   choose frequency from the range (freqmin,freqmax)
      do ifreq=1,intxxx
         freqin=freqmin+frac*dfreq*(ifreq-1)
C     ***   freq greater than zero but less than one
         pi=two*asin(one)
         vxsi=damp
         freq=freqin*two*pi
         den=sqrt(one-vxsi*vxsi)
         vratio=one/den
         freqe=den*freq
         xsiwt=vxsi*freq*dt
         etau=exp(-xsiwt)
         sinwt=sin(freqe*dt)
         coswt=cos(freqe*dt)
         a(1,1)=etau*(vxsi*vratio*sinwt+coswt)
         a(1,2)=etau*sinwt/freqe
         a(2,1)=-etau*freq*vratio*sinwt
         a(2,2)=etau*(coswt-vxsi*vratio*sinwt)
C     
         txsi=(two*vxsi*vxsi-one)/freq/freq/freqe/dt
         xsif=vxsi/freq/freqe
         freqi=one/freq/freq
C     
         b(1,1)=etau*(-(xsif+txsi)*sinwt-
     1        (freqi+two*vxsi*freqi/freq/dt)*coswt)+
     2        two*vxsi*freqi/freq/dt
         b(1,2)=etau*(txsi*sinwt+two*vxsi*freqi/freq/dt*coswt)+
     1        freqi-two*vxsi*freqi/freq/dt
         b(2,1)=etau*(-(freqe*coswt-vxsi*freq*sinwt)*(txsi+xsif)+
     1        (freqe*sinwt+vxsi*freq*coswt)*(freqi+two*vxsi*freqi/
     1        freq/dt))-freqi/dt
         b(2,2)=etau*((freqe*coswt-vxsi*freq*sinwt)*txsi-
     1        (freqe*sinwt+vxsi*freq*coswt)*two*vxsi*freqi/
     2        freq/dt)+freqi/dt
C**   initial conditions
         t=zero
         q(1)=zero
         qp(1)=zero
         qpp(1)=zero
         do it=2,nacc
            t=t+dt
            aq=q(it-1)
            aqp=qp(it-1)
            aacc=acc(it-1)
            q(it)=a(1,1)*aq+a(1,2)*aqp+b(1,1)*aacc+b(1,2)*acc(it)
            qp(it)=a(2,1)*aq+a(2,2)*aqp+b(2,1)*aacc+b(2,2)*acc(it)
C           absolute acceleration record
            qpp(it)=-qp(it)*two*vxsi*freq-q(it)*freq*freq
         end do
         qmax=zero
         qpmax=zero
         qppmax=zero
         do ii=1,nacc
            qabs=abs(q(ii))
            qpabs=abs(qp(ii))
            qppabs=abs(qpp(ii))
            qmax=max(qmax,qabs)
            qpmax=max(qpmax,qpabs)
            qppmax=max(qppmax,qppabs)
         end do
         qm(ifreq)=qmax
         qpm(ifreq)=qpmax
         qppm(ifreq)=qppmax
         fr(ifreq)=freqin
      end do
      do li=1,intxxx
         write(117,211) qppm(li),fr(li),damp
      end do
 211  format(1x,e12.5,',',e12.5,',',e12.5)
      damp=damp+0.02d0
C
      close(101)
      close(117)
C
      stop
      end