CC--PROGRAM TO INTEGRATE EQUATION OF MOTION SUBJECTED TO 
CC  ACCELERATION INPUT. EQUATION IS INTEGRATED FOR A GIVEN
CC  RANGE OF FREQUENCIES AND DAMPING PARAMETER.
CC  MAXIMUM RESPONSE WILL DEFINE RESPONSE SPECTRUM WHICH
CC  CAN BE USED IN MODAL DYNAMICS OPTION TO ESTIMATE MAXIMUM
CC  RESPONSE OF THE CONSIDERED STRUCTURE 
CC    PROGRAM RESPON
      SUBROUTINE ABQMAIN
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(2,2),B(2,2),Q(2500),QP(2500),QPP(2500),
     1ACC(2500),TAB(8),FR(1000),QM(1000),QPM(1000),
     1QPPM(1000)
      parameter(one=1.d0,two=2.d0,zero=0.d0,dt=1.d-2,fact=-386.09d0)
      DATA NACC/1000/
C     DATA NACC/2500/
C     DATA FREQMIN/0.1d0/
C     DATA FREQMAX/30.d0/
      DATA INT/501/
C          STORAGE ALLOCATED FOR A MAXIMUM OF 2500 TIMEPOINTS
C          IN THE EARTHQUAKE HISTORY
CC READ ACCELERATION HISTORY FROM FILE QUAKE.AMP
CC (COPIED FROM CANTILEVER_QUAKEDATA.INP). 
CC DISPLACEMENT SPECTRUM (FOR ABAQUS OPTION \*RESPONSE SPECTRUM)
CC WILL BE WRITTEN TO FILE QUAKE.DIS, VELOCITY SPECTRUM TO FILE QUAKE.VEL.
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 INT DEFINES NUMBER OF POINTS IN FREQUENCY RANGE
C** THIS INPUT ASSUMES THAT CANTILEVER_QUAKEDATA.INP HAS BEEN COPIED
C** TO QUAKE.AMP 
      OPEN(UNIT=100,STATUS='OLD',FILE='QUAKE.AMP')
C** OPEN the output files required, to generate figures 10,11 and 12
      OPEN(UNIT=116,STATUS='NEW',FILE='QUAKE0.DIS')
      OPEN(UNIT=117,STATUS='NEW',FILE='QUAKE2.DIS')
      OPEN(UNIT=118,STATUS='NEW',FILE='QUAKE4.DIS')
      OPEN(UNIT=119,STATUS='NEW',FILE='QUAKER0.DIS')
      OPEN(UNIT=120,STATUS='NEW',FILE='QUAKER2.DIS')
      OPEN(UNIT=121,STATUS='NEW',FILE='QUAKER4.DIS')
      OPEN(UNIT=122,STATUS='NEW',FILE='QUAKE0.VEL')
      OPEN(UNIT=123,STATUS='NEW',FILE='QUAKE2.VEL')
      OPEN(UNIT=124,STATUS='NEW',FILE='QUAKE4.VEL')
      OPEN(UNIT=125,STATUS='NEW',FILE='QUAKER0.VEL')
      OPEN(UNIT=126,STATUS='NEW',FILE='QUAKER2.VEL')
      OPEN(UNIT=127,STATUS='NEW',FILE='QUAKER4.VEL')
      FRAC=ONE/DBLE(INT-1)
C***   INITIATE AMAT,BMAT BEFORE TIME INTEGRATION
      DO 10 I=1,2
      DO 10 J=1,2
      A(I,J)=ZERO
      B(I,J)=ZERO
   10 CONTINUE
C** READ AMPLITUDE DATA AND STORE ON ACC(2500)
      DO 30 I=1,625
      READ(100,*) (TAB(LL),LL=1,8)
      DO 40 K=1,4
      ACC(K+4*(I-1))=TAB(2*K)*FACT
   40 CONTINUE
   30 CONTINUE
C** DO 990 is loop to generate data for frequency ranges of 0.1-30Hz
C**        and 0.01-5Hz, for data shown in figures 10,11 and 12
      DO 990 IFRQ=1,2
      IF (IFRQ.EQ.1) THEN
        FREQMIN=0.1
        FREQMAX=30
        IA=14
        IB=20
      ELSE IF (IFRQ.EQ.2) THEN
        FREQMIN=0.01
        FREQMAX=5
        IA=16
        IB=22
      END IF 
      DAMP=ZERO
      DFREQ=FREQMAX-FREQMIN
C** CHOOSE DAMPING. 
C** DAMPING MUST BE LESS THAN CRITICAL (BETWEEN 0.0 AND 1.0).
      DO 300 IKSI=1,3
      IF(DAMP.GT.ONE) WRITE(6,11)
   11 FORMAT(/,3X,50HTHIS PROGRAM IS WRITTEN FOR UNDERDAMPED CASES ONLY)
C** CHOOSE FREQUENCY FROM THE RANGE (FREQMIN,FREQMAX)
      DO 200 IFREQ=1,INT
      FREQIN=FREQMIN+FRAC*DFREQ*DBLE(IFREQ-1)
C** FREQMIN MUST BE GREATER THAN ZERO
      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
      DO 100 IT=1,NACC
         IF(IT.EQ.1)THEN
C** INITIAL CONDITIONS
            T=0.d0
            Q(1)=0.d0
            QP(1)=0.d0
            QPP(1)=0.d0
         ELSE
            T=T+DT
            Q(IT)=A(1,1)*Q(IT-1)+A(1,2)*QP(IT-1)+B(1,1)*ACC(IT-1)+B(1,2)
     1           *ACC(IT)
            QP(IT)=A(2,1)*Q(IT-1)+A(2,2)*QP(IT-1)+B(2,1)*ACC(IT-1)+
     1           B(2,2)*ACC(IT)
C           absolute acceleration record
            QPP(IT)=-QP(IT)*TWO*VXSI*FREQ-Q(IT)*FREQ*FREQ
         ENDIF
  100 CONTINUE
      QMAX=0.d0
      QPMAX=0.d0
      QPPMAX=0.d0
      DO 110 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)
  110 CONTINUE
      QM(IFREQ)=QMAX
      QPM(IFREQ)=QPMAX
      QPPM(IFREQ)=QPPMAX
      FR(IFREQ)=FREQIN
  200 CONTINUE
      DO 210 LI=1,INT
      WRITE(100+IFRQ+IKSI+IA,211) QM(LI),FR(LI),DAMP
      WRITE(100+IFRQ+IKSI+IB,211) QPM(LI),FR(LI),DAMP
      ID=1
  210 CONTINUE
  211 FORMAT(1X,E12.5,1H,,E12.5,1H,,E12.5)
      DAMP=DAMP+0.02d0
  300 CONTINUE
  990 CONTINUE
C
      CLOSE(100)
      CLOSE(116)
      CLOSE(117)
      CLOSE(118)
      CLOSE(119)
      CLOSE(120)
      CLOSE(121)
      CLOSE(122)
      CLOSE(123)
      CLOSE(124)
      CLOSE(125)
      CLOSE(126)
      CLOSE(127)
C
      STOP
      END