C  * * * * * * * * * * * * * * * * * * * * * * * * *                    
C     SUBROUTINE FOR EARTHQUAKE RESPONSE SPECTRA                        
C  * * * * * * * * * * * * * * * * * * * * * * * * *                    
C                                                                       
C                              CODED BY Y.OHSAKI                        
C                                                                       
C     PURPOSE                                                           
C       TO COMPUTE ABSOLUTE ACCELERATION, RELATIVE VELOCITY OR          
C       RELATIVE DISPLACEMENT RESPONSE SPECTRA OF THE GIVEN ACCELERA-   
C       TION TIME HISTORY FOR SPECIFIED DAMPING FACTORS. THE MAXIMUM    
C       INPUT ACCELERATION, VELOCITY OR DISPLACEMENT ARE ALSO COMPUTED  
C                                                                       
C     USAGE                                                             
C       CALL ERES(NH,H,ND1,NT,T,ND2,DT,NN,DDY,ND3,IND,QMAX,RES)         
C                                                                       
C     DESCRIPTION OF ARGUMENTS                                          
C       NH           - TOTAL NUMBER OF DAMPING FACTORS                  
C       H(ND1)       - DAMPING FACTORS IN DECIMAL FRACTION              
C       ND1          - DIMENSION OF H,RES IN CALLING PROGRAM  ND1.GE.NH 
C       NT           - TOTAL NUMBER OF PERIODS FOR RESPONSE COMPUTATION 
C       T(ND2)       - PERIODS IN SEC FOR RESPONSE COMPUTATION          
C       ND2          - DIMENSION OF T,RES IN CALLING PROGRAM  ND2.GE.NT 
C       DT           - TIME INCREMENT IN THE ACCELERATION TIME HISTORY  
C                      IN SEC                                           
C       NN           - TOTAL NUMBER OF DATA IN THE ACCELERATION TIME    
C                      HISTORY                                          
C       DDY(ND3)     - THE GIVEN ACCELERATION TIME HISTORY IN GALS      
C       ND3          - DIMENSION OF DDY IN CALLING PROGRAM  ND3.GE.NN   
C       IND          - 1 : FOR ABSOLUTE ACCELERATION RESPONSE SPECTRA   
C                      2 : FOR RELATIVE VELOCITY RESPONSE SPECTRA       
C                      3 : FOR RELATIVE DISPLACEMENT RESPONSE SPECTRA   
C       QMAX         - IF IND=1 : MAX. INPUT ACCELERATION               
C                      IF IND=2 : MAX. INPUT VELOCITY                   
C                      IF IND=3 : MAX. INPUT DISPLACEMENT               
C       RES(ND2,ND1) - RESPONSE SPECTRA CORRESPONDING TO THE IND        
C                                                                       
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     
C       NONE                                                            
C                                                                       
      SUBROUTINE ERES(NH,H,ND1,NT,T,ND2,DT,NN,DDY,ND3,IND,QMAX,RES)     
C                                                                       
      DIMENSION H(ND1),T(ND2),DDY(ND3),RES(ND2,ND1)                     
      DIMENSION EMAX(3),RMAX(3)                                         
      PARAMETER (P2=6.283185)                                           
C                                                                       
C     MAXIMA OF INPUT MOTION                                            
C                                                                       
      EMAX(1)=ABS(DDY(1))                                               
      EMAX(2)=0.                                                        
      EMAX(3)=0.                                                        
      DDYF=DDY(1)                                                       
      DYF=0.                                                            
      YF=0.                                                             
      DO 110 M=2,NN                                                     
      DDYM=DDY(M)                                                       
      DY=DYF+(DDYF+DDYM)*DT/2.                                          
      Y=YF+DYF*DT+(DDYF/3.+DDYM/6.)*DT**2                               
      EMAX(1)=AMAX1(EMAX(1),ABS(DDYM))                                  
      EMAX(2)=AMAX1(EMAX(2),ABS(DY))                                    
      EMAX(3)=AMAX1(EMAX(3),ABS(Y))                                     
      DDYF=DDYM                                                         
      DYF=DY                                                            
      YF=Y                                                              
  110 CONTINUE                                                          
      QMAX=EMAX(IND)                                                    
C                                                                       
C     RESPONSE COMPUTATION                                              
C                                                                       
      DO 150 L=1,NH                                                     
      DO 140 K=1,NT                                                     
      IF(T(K).EQ.0.) GO TO 130                                          
      W=P2/T(K)                                                         
      W2=W*W                                                            
      HW=H(L)*W                                                         
      WD=W*SQRT(1.-H(L)**2)                                             
      WDT=WD*DT                                                         
      E=EXP(-HW*DT)                                                     
      CWDT=COS(WDT)                                                     
      SWDT=SIN(WDT)                                                     
      A11= E*(CWDT+HW*SWDT/WD)                                          
      A12= E*SWDT/WD                                                    
      A21=-E*W2*SWDT/WD                                                 
      A22= E*(CWDT-HW*SWDT/WD)                                          
      SS=-HW*SWDT-WD*CWDT                                               
      CC=-HW*CWDT+WD*SWDT                                               
      S1=(E*SS+WD)/W2                                                   
      C1=(E*CC+HW)/W2                                                   
      S2=(E*DT*SS+HW*S1+WD*C1)/W2                                       
      C2=(E*DT*CC+HW*C1-WD*S1)/W2                                       
      S3=DT*S1-S2                                                       
      C3=DT*C1-C2                                                       
      B11=-S2/WDT                                                       
      B12=-S3/WDT                                                       
      B21=(HW*S2-WD*C2)/WDT                                             
      B22=(HW*S3-WD*C3)/WDT                                             
      RMAX(1)=2.*HW*ABS(DDY(1))*DT                                      
      RMAX(2)=ABS(DDY(1))*DT                                            
      RMAX(3)=0.                                                        
      DXF=-DDY(1)*DT                                                    
      XF=0.                                                             
      DO 120 M=2,NN                                                     
      DDYM=DDY(M)                                                       
      DDYF=DDY(M-1)                                                     
      X= A12*DXF+A11*XF+B12*DDYM+B11*DDYF                               
      DX=A22*DXF+A21*XF+B22*DDYM+B21*DDYF                               
      DDX=-2.*HW*DX-W2*X                                                
      RMAX(1)=AMAX1(RMAX(1),ABS(DDX))                                   
      RMAX(2)=AMAX1(RMAX(2),ABS(DX))                                    
      RMAX(3)=AMAX1(RMAX(3),ABS(X))                                     
      DXF=DX                                                            
      XF=X                                                              
  120 CONTINUE                                                          
      RES(K,L)=RMAX(IND)                                                
      GO TO 140                                                         
  130 RES(K,L)=0.                                                       
      IF(IND.EQ.1) RES(K,L)=EMAX(1)                                     
  140 CONTINUE                                                          
  150 CONTINUE                                                          
      RETURN                                                            
      END                                                               
