📄 inecal.for
字号:
$DEBUG
PROGRAM INECAL
DIMENSION ACEL(30000),Q(15),UMAX(200,15),HE(200,15)
DIMENSION R(200,15),XMAX(200,15),DUC(10),C(4)
WRITE(*,*)
WRITE(*,*)
WRITE(*,*)
WRITE(*,*) ' SE ESTAN CALCULANDO LOS ESPECTROS INELASTICOS...'
WRITE(*,*)
WRITE(*,*) ' CUANDO TERMINE EL CALCULO'
WRITE(*,*) ' CIERRE LA VENTANA DOS Window'
WRITE(*,*)
WRITE(*,*)
WRITE(*,*)
OPEN(5,FILE='ACEL.DAT',STATUS='OLD')
OPEN(6,FILE='SDAT.DAT',STATUS='OLD')
OPEN(7,FILE='UMAX.SAL',STATUS='UNKNOWN')
OPEN(8,FILE='HE.SAL',STATUS='UNKNOWN')
OPEN(9,FILE='R.SAL',STATUS='UNKNOWN')
OPEN(10,FILE='XMAX.SAL',STATUS='UNKNOWN')
READ(6,'(I8,F8.3,3I8,3F8.3)')NACEL,DT,NQ,ND,NSAIN,DTSAIN,SAINI
1 ,XI
READ(6,'(20F8.3)') (Q(I),I=1,NQ)
READ(6,'(20F8.3)') (DUC(I),I=1,ND)
DO I=1,NACEL
READ(5,'(E12.4)') ACEL(I)
END DO
DO I=1,NSAIN
W=2.0*3.1416/(SAINI+(I-1)*DTSAIN)
XT0=0
VT0=0
SD=0
DO K=1,NACEL
A0=ACEL(K)
PA=(ACEL(K+1)-ACEL(K))/DT
D=-PA/W/W
CC=(-A0+2*XI*PA/W)/W/W
G1=XT0-CC
G2=(VT0+XI*W*G1+PA/W/W)/W
XT0=EXP(-XI*W*DT)*(G1*COS(W*DT)+G2*SIN(W*DT))+CC+D*DT
Vt0=EXP(-XI*W*DT)*(-G1*W*SIN(W*DT)+G2*W*COS(W*DT))
1 -XI*W*EXP(-XI*W*DT)*(G1*COS(W*DT)+G2*SIN(W*DT))+D
IF (SD .LT. ABS(XT0)) THEN
SD=ABS(XT0)
END IF
END DO
R(I,1)=W*W*SD
XMAX(I,1)=SD
UMAX(I,1)=1.0
DO J=1,NQ
HE(I,J+1)=0.0
UMAX(I,J+1)=1.0
R(I,J+1)=R(I,1)/Q(J)
X0=0
V0=0
A0=0
R0=0
P0=0
IND=1
DO K=1,NACEL
IF (IND .EQ. 1) THEN
DP=(ACEL(K)-P0)-A0*DT/2.*(4*XI*W+DT*W*W)-W*W*DT*V0
DK=1.+XI*W*DT+W*W*DT*DT/4.
DA=DP/DK
XU=X0+DT*DT/4.*(DA+2.*A0)+DT*V0
RU=(Xu-X0)*W*W+R0
IF (ABS(RU) .GE. R(I,J+1)) THEN
RU=RU/ABS(RU)*R(I,J+1)
DX=(RU-R0)/W/W
C(1)=4.*DX
C(2)=4.*(XI*W*DX-V0)
C(3)=W*W*DX-4.*XI*W*V0-2.*A0
C(4)=(P0-ACEL(K))/DT
CALL RAICES(3,C,DT,DTT)
DP=(ACEL(K)-P0)/DT*DTT-A0*DTT/2.
1 *(4.*XI*W+DTT*W*W)-W*W*DTT*V0
DK=1.+XI*W*DTT+W*W*DTT*DTT/4.
DA=DP/DK
XU=X0+DTT*DTT/4.*(DA+2.*A0)+DTT*V0
HE(I,J+1)=HE(I,J+1)+(RU+R0)*(XU-X0)/2.
R0=RU
X0=XU
V0=V0+DTT/2.*(DA+2.*A0)
A0=A0+DA
P0=P0+(ACEL(K)-P0)/DT*DTT
DTT=DT-DTT
DA=ACEL(K)-P0
XU=X0+DTT*DTT/4.*(DA+2.*A0)+DTT*V0
VU=V0+DTT/2.*(DA+2.*A0)
IF (VU*V0 .LE. 0) THEN
DV=-V0
C(1)=2.*DV
C(2)=-2.*A0
C(3)=(P0-ACEL(K))/DTT
CALL RAICES (2,C,DTT,DTTT)
DA=(ACEL(K)-P0)/DTT*DTTT
XU=X0+DTTT*DTTT/4.*(DA+2.*A0)+DTTT*V0
IF (ABS(XU)/R(I,J)*W*W .GT. UMAX(I,J)) THEN
UMAX(I,J+1)=ABS(XU)/R(I,J+1)*W*W
XMAX(I,J+1)=ABS(XU)
END IF
HE(I,J+1)=HE(I,J+1)+R0*(XU-X0)
X0=XU
V0=0.
A0=A0+DA
P0=P0+(ACEL(K)-P0)/DTT*DTTT
DTTT=DTT-DTTT
DP=(ACEL(K)-P0)-A0*DTTT/2.*(4.*XI*W+DTTT*W*W)
1 -W*W*DTTT*V0
DK=1.+XI*W*DTTT+W*W*DTTT*DTTT/4.
DA=DP/DK
XU=X0+DTTT*DTTT/4.*(DA+2.*A0)+DTTT*V0
RU=(XU-X0)*W*W+R0
C AQUI PUEDE HABER PROBLEMAS
HE(I,J+1)=HE(I,J+1)+(RU+R0)*(XU-X0)/2.
R0=RU
X0=XU
V0=V0+DTTT/2.*(DA+2.*A0)
A0=A0+DA
P0=ACEL(K)
IND=1
ELSE
HE(I,J+1)=HE(I,J+1)+R0*(XU-X0)
X0=XU
V0=V0+DTT/2.*(DA+2.*A0)
A0=A0+DA
P0=ACEL(K)
IND=2
END IF
ELSE
HE(I,J+1)=HE(I,J+1)+(RU+R0)*(XU-X0)/2.
R0=RU
X0=XU
V0=V0+DT/2.*(DA+2.*A0)
A0=A0+DA
P0=ACEL(K)
END IF
ELSE
DA=ACEL(K)-P0
XU=X0+DT*DT/4.*(DA+2.*A0)+DT*V0
VU=V0+DT/2.*(DA+2.*A0)
IF (VU*V0 .LE. 0) THEN
DV=-V0
C(1)=2.*DV
C(2)=-2.*A0
C(3)=(P0-ACEL(K))/DT
CALL RAICES (2,C,DT,DTT)
DA=(ACEL(K)-P0)/DT*DTT
XU=X0+DTT*DTT/4.*(DA+2.*A0)+DTT*V0
IF (ABS(XU)/R(I,J+1)*W*W .GT. UMAX(I,J+1)) THEN
UMAX(I,J+1)=ABS(XU)/R(I,J+1)*W*W
XMAX(I,J+1)=ABS(XU)
END IF
HE(I,J+1)=HE(I,J+1)+R0*(XU-X0)
X0=XU
V0=0
A0=A0+DA
P0=P0+(ACEL(K)-P0)/DT*DTT
DTT=DT-DTT
DP=(ACEL(K)-P0)-A0*DTT/2.*(4*XI*W+DTT*W*W)
1 -W*W*DTT*V0
DK=1.+XI*W*DTT+W*W*DTT*DTT/4.
DA=DP/DK
XU=X0+DTT*DTT/4.*(DA+2.*A0)+DTT*V0
RU=(XU-X0)*W*W+R0
C AQUI PUEDE HABER PROBLEMAS
HE(I,J+1)=HE(I,J+1)+(RU+R0)*(XU-X0)/2.
R0=RU
X0=XU
V0=V0+DTT/2.*(DA+2.*A0)
A0=A0+DA
P0=ACEL(K)
IND=1
ELSE
HE(I,J+1)=HE(I,J+1)+R0*(XU-X0)
X0=XU
V0=V0+DT/2.*(DA+2.*A0)
A0=A0+DA
P0=ACEL(K)
END IF
END IF
END DO
IF ((UMAX(I,J+1) .GT. DUC(ND)) .AND.
1 (Q(J) .GE. DUC(ND))) THEN
DO L=J+1,NQ
UMAX(I,L+1)=UMAX(I,J+1)+L+1
HE(I,L+1)=HE(I,J+1)
R(I,L+1)=R(I,1)/Q(L)
XMAX(I,L+1)=XMAX(I,J+1)
END DO
GOTO 500
END IF
END DO
500 CONTINUE
WRITE(7,'(21E12.5)') (UMAX(I,J),J=1,NQ+1)
WRITE(8,'(21E12.5)') (HE(I,J),J=1,NQ+1)
WRITE(9,'(21E12.5)') (R(I,J),J=1,NQ+1)
WRITE(10,'(21E12.5)') (XMAX(I,J),J=1,NQ+1)
END DO
CLOSE (5)
CLOSE (6)
CLOSE (7)
CLOSE (8)
CLOSE (9)
CLOSE (10)
STOP
END
SUBROUTINE RAICES(NG,C,DT,DTT)
DIMENSION C(NG+1)
REAL*8 Q,P,D,DD,R,FI
IF (NG .EQ. 2) THEN
IF (C(3) .NE. 0) THEN
PP=SQRT(C(2)*C(2)-4.*C(3)*C(1))/2./C(3)
QQ=-C(2)/2./C(3)
R1=QQ+PP
R2=QQ-PP
IF ((R1 .GE. 0.) .AND. (R1 .LE. DT)) THEN
DTT=R1
ELSE
DTT=R2
END IF
ELSE
DTT=-C(1)/C(2)
END IF
RETURN
ELSE
IF (C(4) .EQ. 0) THEN
IF (C(3) .NE. 0) THEN
PP=SQRT(C(2)*C(2)-4.*C(3)*C(1))/2./C(3)
QQ=-C(2)/2./C(3)
R1=QQ+PP
R2=QQ-PP
IF ((R1 .GE. 0.) .AND. (R1 .LE. DT)) THEN
DTT=R1
ELSE
DTT=R2
END IF
ELSE
DTT=-C(1)/C(2)
END IF
RETURN
ELSE
A=C(3)/3./C(4)
Q=(C(3)/C(4))**3/27.-C(3)*C(2)/6./C(4)/C(4)+C(1)/C(4)/2.
P=(3.*C(4)*C(2)-C(3)*C(3))/9./C(4)/C(4)
D= Q*Q+P*P*P
IF (D .GT. 0.) THEN
DD= SQRT(D)
DTT=(-Q+DD)/ABS(-Q+DD)*(ABS(-Q+DD))**(1./3.)+
1 (-Q-DD)/ABS(-Q-DD)*(ABS(-Q-DD))**(1./3.)-A
ELSE
R=Q/ABS(Q)*SQRT(ABS(P))
FI=ATAN(SQRT(R**6.-Q*Q)/ABS(Q))
R1=-2.*R*COS(FI/3.)-A
R2=2.*R*COS(1.0472-FI/3.)-A
R3=2.*R*COS(1.0472+FI/3.)-A
IF ((R1 .GE. 0) .AND. (R1 .LE. DT)) THEN
DTT=R1
ELSE
IF ((R2 .GE. 0) .AND. (R2 .LE. DT)) THEN
DTT=R2
ELSE
DTT=R3
END IF
END IF
END IF
RETURN
END IF
END IF
STOP
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -