C
C     QSRS.F   VER 1.1  JANUARY 3, 2012
C
C     BY TOM IRVINE
C     Email:  tomirvine@aol.com
C
C     Shock response spectrum of an arbitrary acceleration
C     base input time history supplied by the user.
C
C
C     Variables:
C
C        T = input time
C        Y = input base acceleration
C
C23456789
C
      INTEGER*4 MAX,NUM
      INTEGER MAX_NF
      INTEGER IRD
      DOUBLE PRECISION PI,TPI
      PARAMETER(PI=3.141592653589793D0,TPI=2D0*PI,MAX=4000000)
      PARAMETER(MAX_NF=10000)
      DOUBLE PRECISION AMP,TT
      DOUBLE PRECISION DTMIN,DTMAX,DT,SRMIN,SR,SRMAX
      DOUBLE PRECISION DAM
      DOUBLE PRECISION OCT
      DOUBLE PRECISION TEMP
      DOUBLE PRECISION FSTART,FEND
      DOUBLE PRECISION FFMAX,XXMAX
      CHARACTER :: A_filename*20,B_filename*20,C_filename*20
      INTEGER IOstatus,NSR,IDAMP,ALGO
      INTEGER stat_alloc
      DOUBLE PRECISION, ALLOCATABLE :: T(:),Y(:)
      DOUBLE PRECISION, ALLOCATABLE :: FN(:)
      DOUBLE PRECISION, ALLOCATABLE :: a1(:),a2(:)
      DOUBLE PRECISION, ALLOCATABLE :: b1(:),b2(:),b3(:)
      DOUBLE PRECISION, ALLOCATABLE :: rd_a1(:),rd_a2(:)
      DOUBLE PRECISION, ALLOCATABLE :: rd_b1(:),rd_b2(:),rd_b3(:)
      DOUBLE PRECISION, ALLOCATABLE :: XMAX(:),XMIN(:)
      DOUBLE PRECISION, ALLOCATABLE :: RMAX(:),RMIN(:)
      CHARACTER*1 tab
      tab=CHAR(9)
C
      PRINT '("")'
      PRINT '("The base input must have two columns: ")'
      PRINT '("  time(sec) & accel(G) ")'
      PRINT '("")'
      PRINT '( "Enter the base input acceleration filename.")'
      READ(*,*)A_filename
C
      OPEN(UNIT=10,FILE=A_filename,STATUS='OLD',ACTION='READ')
C
      DO I=1,MAX
            READ(10,*,IOSTAT=IOstatus) TT,AMP
            IF (IOstatus.EQ.0) THEN
                NUM=I
            ELSE
                EXIT
            ENDIF
      ENDDO
C
      ALLOCATE(T(NUM))
      ALLOCATE(Y(NUM))
C
      REWIND(10)
C
      DO I=1,NUM
            READ(10,*,IOSTAT=IOstatus) T(I),Y(I)
      ENDDO
C
      DTMIN=1.0D+90
      DTMAX=0.
C
      DO I=1,NUM-1
         DT=T(I+1)-T(I)
         IF(DT.GT.DTMAX)THEN
            DTMAX=DT
         ENDIF
         IF(DT<DTMIN)THEN
            DTMIN=DT
         ENDIF
      ENDDO
C
      SRMAX=1/DTMIN
      SRMIN=1/DTMAX
      DT=(T(NUM)-T(1))/(NUM-1)
      SR=1/DT
C
      PRINT '("")'
      PRINT '(" ",I10," lines read ")',NUM
C
      PRINT '("")'
      PRINT '(" DTMIN=",G10.4," sec ")',DTMIN
      PRINT '(" DT   =",G10.4," sec ")',DT
      PRINT '(" DTMAX=",G10.4," sec ")',DTMAX
      PRINT '("")'
      PRINT '(" SRMAX=",G10.4)',SRMAX
      PRINT '(" SR   =",G10.4)',SR
      PRINT '(" SRMIN=",G10.4)',SRMIN
C
      IF( DTMIN.LT.0.99 *DTMAX)THEN
        PRINT '("")'
        PRINT '(" *** Time Step Warning *** ")'
        PRINT '(" The time step must be constant.")'
        PRINT '(" Enter new sample rate?  1=yes  2=no ")'
        READ(*,*)NSR
        IF(NSR.EQ.1)THEN
            PRINT '("")'
            PRINT '("Enter new sample rate ")'
            READ(*,*)SR
            DT=1/SR
        ELSE
           GOTO 100
        ENDIF
      ENDIF
C
      PRINT '("")'
      PRINT '(" Enter Damping format.")'
      PRINT '("  1= Damping ratio   2= Q  ")'
      READ(*,*) IDAMP
C
      IF(IDAMP.EQ.1)THEN
         PRINT '("")'
         PRINT '(" Enter Damping ratio (typically 0.05)")'
         READ(*,*)DAM
      ELSE
         PRINT '("")'
         PRINT '(" Enter amplification factor Q (typically 10)")'
         READ(*,*)DAM
	     DAM = 1./(2.*DAM)
      ENDIF
C
      PRINT '("")'
      PRINT '(" Enter starting frequency (Hz) ")'
      READ(*,*)FSTART
C
      FEND=SR/8.
      NFREQ=12*(INT(LOG(FEND/FSTART)/LOG(2.)))
      PRINT '("")'
      PRINT '(" NFREQ=",I4)',NFREQ
C
      ALLOCATE(FN(NFREQ))
C
      OCT=2.**(1./12.)
      FN(1)=FSTART
      DO J=2,NFREQ
         FN(J)=FN(J-1)*OCT
      ENDDO
C
      PRINT '("")'
      PRINT '(" Include residual?  1=yes 2=no ")'
      READ(*,*)IRE
C
      ALLOCATE(a1(NFREQ))
      ALLOCATE(a2(NFREQ))
      ALLOCATE(b1(NFREQ))
      ALLOCATE(b2(NFREQ))
      ALLOCATE(b3(NFREQ))
      ALLOCATE(rd_a1(NFREQ))
      ALLOCATE(rd_a2(NFREQ))
      ALLOCATE(rd_b1(NFREQ))
      ALLOCATE(rd_b2(NFREQ))
      ALLOCATE(rd_b3(NFREQ))
      ALLOCATE(XMAX(NFREQ))
      ALLOCATE(XMIN(NFREQ))
      ALLOCATE(RMAX(NFREQ))
      ALLOCATE(RMIN(NFREQ))
C
      CALL ALGORITHM(DT,a1,a2,b1,b2,b3,rd_a1,rd_a2,rd_b1,rd_b2,rd_b3,
     $NFREQ,DAM,TPI,FN,IRD)
C
      CALL SRS_ENGINE(a1,a2,b1,b2,b3,rd_a1,rd_a2,rd_b1,rd_b2,rd_b3,FN,
     $Y,XMAX,XMIN,RMAX,RMIN,NUM,NFREQ,IRE,IRD)
C
      PRINT '( "Enter the acceleration output filename. ")'
      READ(*,*)B_filename
      PRINT '("")'
      PRINT '(" Output Format:  fn(Hz)  Max(G)  MIN(G)")'
C
      OPEN(UNIT=11,FILE=B_filename,ACTION='WRITE')
C
      XXMAX=0
      FFMAX=0
C
      DO I=1,NFREQ
        XMIN(I)=ABS(XMIN(I))
        WRITE(11,20)FN(I),tab,XMAX(I),tab,XMIN(I)
C
        IF(XMAX(I).GT.XXMAX)THEN
           XXMAX=XMAX(I)
           FFMAX=FN(I)
        ENDIF
        IF(XMIN(I).GT.XXMAX)THEN
           XXMAX=XMIN(I)
           FFMAX=FN(I)
        ENDIF
C
      ENDDO
20    FORMAT(3X,F8.2,A1,E10.4,A1,E10.4)
C
      PRINT '(" ")'
      PRINT '(" Absolute Maximum Acceleration Response = ")'
      WRITE(*,66) XXMAX,FFMAX
66    FORMAT(5X,G9.4,'G  at ',G9.4,'Hz')
C
      IF(IRD.EQ.1)THEN
         PRINT '(" ")'
         PRINT '(" ")'
         PRINT '( "Enter the relative displacement output filename. ")'
         READ(*,*)C_filename
         PRINT '("")'
         PRINT '(" Output Format:  fn(Hz)  Max(inch)  MIN(inch)")'
C
         OPEN(UNIT=12,FILE=C_filename,ACTION='WRITE')
C
         XXMAX=0
         FFMAX=0
C
         DO I=1,NFREQ
C
           RMIN(I)=ABS(RMIN(I))
           WRITE(12,20)FN(I),tab,RMAX(I),tab,RMIN(I)
C
           IF(RMAX(I).GT.XXMAX)THEN
              XXMAX=RMAX(I)
              FFMAX=FN(I)
           ENDIF
           IF(RMIN(I).GT.XXMAX)THEN
              XXMAX=RMIN(I)
              FFMAX=FN(I)
           ENDIF
C
         ENDDO
C
         PRINT '(" ")'
         PRINT '(" Absolute Maximum Relative Displacement Response =")'
         WRITE(*,67) XXMAX,FFMAX
67       FORMAT(5X,G9.3,' inch  at ',G9.4,'Hz')
      ENDIF
C
      DEALLOCATE(T,stat=stat_alloc)
      DEALLOCATE(Y,stat=stat_alloc)
      DEALLOCATE(FN,stat=stat_alloc)
      DEALLOCATE(a1,stat=stat_alloc)
      DEALLOCATE(a2,stat=stat_alloc)
      DEALLOCATE(b1,stat=stat_alloc)
      DEALLOCATE(b2,stat=stat_alloc)
      DEALLOCATE(b3,stat=stat_alloc)
C
      DEALLOCATE(rd_a1,stat=stat_alloc)
      DEALLOCATE(rd_a2,stat=stat_alloc)
      DEALLOCATE(rd_b1,stat=stat_alloc)
      DEALLOCATE(rd_b2,stat=stat_alloc)
      DEALLOCATE(rd_b3,stat=stat_alloc)
      DEALLOCATE(XMAX,stat=stat_alloc)
      DEALLOCATE(XMIN,stat=stat_alloc)
      DEALLOCATE(RMAX,stat=stat_alloc)
      DEALLOCATE(RMIN,stat=stat_alloc)
C
100   STOP
      END
C
      SUBROUTINE SRS_ENGINE(a1,a2,b1,b2,b3,rd_a1,rd_a2,rd_b1,rd_b2,
     $rd_b3,FN,Y,XMAX,XMIN,RMAX,RMIN,NUM,NFREQ,IRE,IRD)
      INTEGER NFREQ
      DOUBLE PRECISION Y(NUM)
      DOUBLE PRECISION X(NFREQ),XB(NFREQ),XBB(NFREQ)
	  DOUBLE PRECISION XMAX(NFREQ),XMIN(NFREQ)
      DOUBLE PRECISION R(NFREQ),RB(NFREQ),RBB(NFREQ)
	  DOUBLE PRECISION RMAX(NFREQ),RMIN(NFREQ)
C
 	  DOUBLE PRECISION YB(NFREQ),YBB(NFREQ)
C
      DOUBLE PRECISION a1(NFREQ),a2(NFREQ)
      DOUBLE PRECISION b1(NFREQ),b2(NFREQ),b3(NFREQ)
      DOUBLE PRECISION rd_a1(NFREQ),rd_a2(NFREQ)
      DOUBLE PRECISION rd_b1(NFREQ),rd_b2(NFREQ),rd_b3(NFREQ)
C
      DOUBLE PRECISION YY
      DOUBLE PRECISION THREE_QUARTERS_PERIOD
      DOUBLE PRECISION FN(NFREQ)
      INTEGER*4 NUM
      INTEGER IRE
C
      DO I=1,NFREQ
            X(I)=0
           XB(I)=0
          XBB(I)=0
           YB(I)=0
	      YBB(I)=0
C
         XMAX(I)=0.
         XMIN(I)=0.
C
      ENDDO
C
C  Digital Recursive Filtering Relationship
C
      DO I=1,NUM
         DO J=1,NFREQ
C
            X(J)= a1(J)*XB(J)+a2(J)*XBB(J)
     $           +b1(J)*Y(I)+b2(J)*YB(J)+b3(J)*YBB(J)
C
	        IF(X(J) .GT. XMAX(J))THEN
		        XMAX(J)=X(J)
	        ENDIF
	        IF(X(J) .LT. XMIN(J))THEN
		        XMIN(J)=X(J)
            ENDIF
C
	        XBB(J)=XB(J)
	         XB(J)=X(J)
	        YBB(J)=YB(J)
	         YB(J)=Y(I)
C
         ENDDO
      ENDDO
C
      IF(IRE.EQ.1)THEN
C
         DO J=1,NFREQ
            THREE_QUARTERS_PERIOD=0.75/FN(J)
            YY=0
            DO I=1,NUM
               TR=(I-1)*DT
               IF(TR.GT.THREE_QUARTERS_PERIOD)THEN
                  EXIT
               ENDIF
C
            X(J)= a1(J)*XB(J)+a2(J)*XBB(J)
     $           +b1(J)*YY+b2(J)*YB(J)+b3(J)*YBB(J)
C
               IF(X(J) .GT. XMAX(J))THEN
		          XMAX(J)=X(J)
	           ENDIF
	           IF(X(J) .LT. XMIN(J))THEN
		          XMIN(J)=X(J)
               ENDIF
C
	           XBB(J)=XB(J)
                XB(J)=X(J)
	           YBB(J)=YB(J)
                YB(J)=YY
C
            ENDDO
         ENDDO
      ENDIF
C
      IF(IRD.EQ.1)THEN
         DO I=1,NFREQ
              R(I)=0
             RB(I)=0
            RBB(I)=0
             YB(I)=0
	        YBB(I)=0
C
           RMAX(I)=0.
           RMIN(I)=0.
C
        ENDDO
C
        DO I=1,NUM
           DO J=1,NFREQ
C
            R(J)= rd_a1(J)*RB(J)+rd_a2(J)*RBB(J)
     $           +rd_b1(J)*Y(J)+rd_b2(J)*YB(J)+rd_b3(J)*YBB(J)
C
               IF(R(J) .GT. RMAX(J))THEN
		          RMAX(J)=R(J)
	           ENDIF
	           IF(R(J) .LT. RMIN(J))THEN
		          RMIN(J)=R(J)
               ENDIF
C
	          RBB(J)=RB(J)
	           RB(J)=R(J)
	          YBB(J)=YB(J)
	           YB(J)=Y(I)
C
           ENDDO
        ENDDO
C
        IF(IRE.EQ.1)THEN
C
           DO J=1,NFREQ
              THREE_QUARTERS_PERIOD=0.75/FN(J)
              YY=0
              DO I=1,NUM
                 TR=(I-1)*DT
                 IF(TR.GT.THREE_QUARTERS_PERIOD)THEN
                    EXIT
                 ENDIF
C
            R(J)= rd_a1(J)*RB(J)+rd_a2(J)*RBB(J)
     $           +rd_b1(J)*YY+rd_b2(J)*YB(J)+rd_b3(J)*YBB(J)
C
                 IF(R(J) .GT. RMAX(J))THEN
		            RMAX(J)=R(J)
	             ENDIF
	             IF(R(J) .LT. RMIN(J))THEN
		            RMIN(J)=R(J)
                 ENDIF
C
	             RBB(J)=RB(J)
                  RB(J)=R(J)
	             YBB(J)=YB(J)
                  YB(J)=YY
C
              ENDDO
           ENDDO
        ENDIF
C
        RMAX=RMAX*386.
        RMIN=RMIN*386.
C
      ENDIF
C
      RETURN
      END
C
      SUBROUTINE ALGORITHM(DT,a1,a2,b1,b2,b3,rd_a1,rd_a2,rd_b1,rd_b2,
     $rd_b3,NFREQ,DAM,TPI,FN,IRD)
C
      INTEGER NFREQ,IRD
      DOUBLE PRECISION a1(NFREQ),a2(NFREQ)
      DOUBLE PRECISION b1(NFREQ),b2(NFREQ),b3(NFREQ)
      DOUBLE PRECISION rd_a1(NFREQ),rd_a2(NFREQ)
      DOUBLE PRECISION rd_b1(NFREQ),rd_b2(NFREQ),rd_b3(NFREQ)
      DOUBLE PRECISION FN(NFREQ)
      DOUBLE PRECISION DAM
      DOUBLE PRECISION DT,TPI
      DOUBLE PRECISION OMEGA,OMEGAD
      DOUBLE PRECISION COSD,SIND,DOMEGADT
      DOUBLE PRECISION C,E,K,S,Sp
      INTEGER ALGO
C
      PRINT '("")'
      PRINT '(" Select algorithm ")'
	  PRINT '("   1=Kelly-Richman  2=Smallwood ")'
      READ(*,*)ALGO
C
      IF(ALGO.EQ.1)THEN
         DO J=1,NFREQ
			OMEGA=TPI*FN(J)
			OMEGAD=OMEGA*SQRT(1.-(DAM**2))
			COSD=COS(OMEGAD*DT)
			SIND=SIN(OMEGAD*DT)
			DOMEGADT=DAM*OMEGA*DT
			a1(J)=2.*EXP(-DOMEGADT)*COSD
			a2(J)=-EXP(-2.*DOMEGADT)
			b1(J)=2.*DOMEGADT
			b2(J)=OMEGA*DT*EXP(-DOMEGADT)
			b2(J)=b2(J)*((OMEGA/OMEGAD)*(1.-2.*(DAM**2))*SIND -2.*DAM*COSD)
			b3(J)=0.
         ENDDO
      ELSE
         DO J=1,NFREQ
			OMEGA=TPI*FN(J)
			OMEGAD=OMEGA*SQRT(1.-(DAM**2))
			E=EXP(-DAM*OMEGA*DT)
			K=OMEGAD*DT
			C=E*COS(K)
			S=E*SIN(K)
			Sp=S/K
			a1(J)=2*C
			a2(J)=-(E**2)
			b1(J)=1.-Sp
			b2(J)=2.*(Sp-C)
			b3(J)=(E**2)-Sp
         ENDDO
      ENDIF
C
      PRINT '("")'
      PRINT '(" Calculate Relative Displacement?  1=yes 2=no ")'
      READ(*,*)IRD
C
      IF(IRD.EQ.1)THEN
         DO J=1,NFREQ
			OMEGA=TPI*FN(J)
			OMEGAD=OMEGA*SQRT(1.-(DAM**2))
			COSD=COS(OMEGAD*DT)
			SIND=SIN(OMEGAD*DT)
			DOMEGADT=DAM*OMEGA*DT
			a1(J)=2.*EXP(-DOMEGADT)*COSD
			a2(J)=-EXP(-2.*DOMEGADT)
			b1(J)=2.*DOMEGADT
			b2(J)=OMEGA*DT*EXP(-DOMEGADT)
			b2(J)=b2(J)*((OMEGA/OMEGAD)*(1.-2.*(DAM**2))*SIND -2.*DAM*COSD)
			b3(J)=0.
C
			rd_a1(J)=2.*EXP(-DOMEGADT)*COSD
            rd_a2(J)=-EXP(-2.*DOMEGADT)
            rd_b1(J)=0.
            rd_b2(J)=-(DT/OMEGAD)*EXP(-DOMEGADT)*SIND
            rd_b3(J)=0
C
         ENDDO
      ENDIF
C
      RETURN
      END