程序11.9 RESP

1、主程序为：

DIMENSION DDY(24256),ACC(24256),VEL(24256),DIS(24256)
DIMENSION TIM(24256)
INTEGER::NN
REAL::DT
CHARACTER*5::SNN,SDT
DATA H/0.05/,T/0.3/
!
open(5,file='ua0259.DAT')
open(7,file='OUTPUT.DAT')
!------------INPUT DATA CONVERSION------------
DO J=1,11
READ(5,*)
ENDDO
!
READ(5,501) SNN,SDT
READ(SNN,502) NN
READ(SDT(1:5),503) DT
WRITE(*,*)SNN,SDT
WRITE(*,*)NN,DT
!------------INPUT CALCULATED DATA------------
DO J=1,4
READ(5,*)
ENDDO
!
READ(5,601)(DDY(M),M=1,NN)
W=6.283185/T
CALL RESP(H,W,DT,NN,DDY,ACC,VEL,DIS,24256,SA,SV,SD)
!
DO J=1,NN
TIM(J)=REAL(J-1)*DT 
ENDDO
!
WRITE (7,701) SA,SV,SD,(MK,TIM(MK),ACC(MK),VEL(MK),DIS(MK),MK=1,NN)
STOP
!
501 FORMAT(T17,A5,TR32,A5)
502 FORMAT(I5)
503 FORMAT(F5.3)
!
601 FORMAT(8ES15.6)
!
701 FORMAT(T35,'EXAMPLE WAVE'//T15,'--ACCELERATION, VELOCITY AND DISPLACEMENT SPECTRUM--'&
          ///T8,'SA',TR5,'SV',TR5,'SD'//(T5,F10.3,TR5,F10.3,TR5,F10.3)&
          ///T8,'MK',TR5,'TIME',TR5,'ACCELERATION SPECTRUM',TR2,'VELOCITY SPECTRUM',&
          TR2,'DISPLACEMENT SPECTRUM'//(T5,I5,TR3,F10.3,TR7,F10.3,TR5,F10.3,TR5,F10.3))
END

2、子程序为：

!******************************************************************
!          SUBROUTINE FOR RESPONSE OF SINGLE-DOF SYSTEM
!******************************************************************
!
!
!
SUBROUTINE RESP(H,W,DT,NN,DDY,ACC,VEL,DIS,ND,SA,SV,SD)
!
DIMENSION DDY(ND),ACC(ND),VEL(ND),DIS(ND)
!
W2=W*W
HW=H*W
WD=W*SQRT(1.-H*H)
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
ACC(1)=2.*H*W*DDY(1)*DT
VEL(1)=-DDY(1)*DT
DIS(1)=0.
DX=VEL(1)
X=0.
SA=0.
SV=0.
SD=0.
DO M=2,NN
  DXF=DX
  XF=X
  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
  ACC(M)=DDX
  VEL(M)=DX
  DIS(M)=X
  SA=AMAX1(SA,ABS(DDX))
  SV=AMAX1(SV,ABS(DX))
  SD=AMAX1(SD,ABS(X))
ENDDO
RETURN
END

3、输出结果（见附件 OUTPUT.TXT），结果竟然是无穷大，无比郁闷，不知道哪里出了问题.

OUTPUT.DAT

转载本文请联系原作者获取授权，同时请注明本文来自黄艳科学网博客。
链接地址：https://blog.sciencenet.cn/blog-456941-457269.html

上一篇：阅读手记之六——地震动的谱分析入门
下一篇：阅读手记之八——地震动的谱分析入门