! ***********************************
! *	STATIC ANALYSIS FOR PLANE FRAME *
! *	 (METHOD OF LATTER TREATMENT)	*
! ***********************************
! 張琛 土木0605 0604050830
! Main program reads the control parameters and organizes
! the whole calculation by calling subroutines
  PROGRAM PFL
  DIMENSION X(50),Y(50),IJ(50,2),A(50),ZI(50),JR(20,4),PJ(50,3),PF(50,4),TK(150,150),P(150)
  CHARACTER*12 INDAT,OUTDAT
  WRITE(*,*)'please input primary data file name!'
  READ(*,'(A12)')INDAT
  WRITE(*,*)'please input calculation result file name!'
  READ(*,'(A12)')OUTDAT
  OPEN(1,FILE=INDAT,STATUS='OLD')
  OPEN(2,FILE=OUTDAT,STATUS='new')
  READ(1,*)NE,NJ,NR,NP,NF,E
  WRITE(2,10)NE,NJ,NR,NP,NF,E
  10 FORMAT (3X,'PLEASE FRAME STRUCTURE ANALYSIS'/5X,'NE=',I2,8X,'NJ=',I2,8X,'NR=',I2,/5X,'NP=',I2,8X,'NF=',I2,8X,'E=',E12.4)
  N=NJ*3
  CALL INPUT (NE,NJ,NR,NP,NF,X,Y,IJ,A,ZI,JR,PJ,PF)
  CAll TSM(NE,NJ,E,X,Y,IJ,A,ZI,TK,N)
  CALL JLP(NE,NJ,NP,NF,X,Y,IJ,PJ,PF,P,N)
  CALL ISC(NR,JR,TK,P,N)
   CALL GAUSS (TK,P,N)
  CALL MVN (NE,NJ,NF,E,X,Y,IJ,A,ZI,PF,P,N)
  CLOSE(1)
  CLOSE(2)
  END
! read and print primary data.
  SUBROUTINE INPUT(NE,NJ,NR,NP,NF,X,Y,IJ,A,ZI,JR,PJ,PF)
  DIMENSION X(NJ),Y(NJ),IJ(NE,2),A(NE),ZI(NE),JR(NR,4),PJ(NP,3),PF(NF,4)
 READ(1,*)(X(I),Y(I),I=1,NJ)
 READ(1,*)(IJ(I,1),IJ(I,2),A(I),ZI(I),I=1,NE)
 READ(1,*)((JR(I,J),J=1,4),I=1,NR)
 IF(NP.GT.0)READ(1,*)((PJ(I,J),J=1,3),I=1,NP)
 IF(NF.GT.0)READ(1,*)((PF(I,J),J=1,4),I=1,NF)
 WRITE(2,10)(I,X(I),Y(I),I=1,NJ)
 WRITE(2,20)(I,IJ(I,1),IJ(I,2),A(I),ZI(I),I=1,NE)
 WRITE(2,30)((JR(I,J),J=1,4),I=1,NR)
 IF(NP.GT.0)WRITE(2,40)((PJ(I,J),J=1,3),I=1,NP)
 IF(NF.GT.0)WRITE(2,50)((PF(I,J),J=1,4),I=1,NF)
 10 FORMAT (/2X,'COORDINATES OF JOINT'/6X,'JOINT',12X,'X',12X,'Y'/(6X,I4,5X,2F12.4))
 20 FORMAT (/2X,'INFORMATION OF ELEMENTS'/6X,'ELEMENT',4X,'JOINT-I',4X,'JOINT-J',10X,'A',12X,'ZI'/(2X,3I10,6X,2F12.6))
 30 FORMAT (/2X,'INFORMATION OF RESTRICTION'/6X,'RES.-JOINT',7X,'XR',8X,'YR',7X,'CETA'/(4X,4I10))
 40 FORMAT (/2X,'JOINT LOAD'/6X,'JOINT'8X,'XYM',12X,'LOAD'/(6X,F5.0,6X,F5.0,6X,F12.4))
 50 FORMAT (/2X,'NON-JOINT LOAD'/6X,'ELEMENT',8X,'TYPE',8X,'LOAD',12X,'C'/(6X,F6.0,6X,F6.0,4X,2F12.4))
 END
!Calculate length,sine and cosine of member.
SUBROUTINE LSC(M,NE,NJ,X,Y,IJ,BL,SI,CO)
DIMENSION X (NJ),Y(NJ),IJ(NE,2)
I=IJ(M,1)
J=IJ(M,2)
DX=X(J)-X(I)
DY=Y(J)-Y(I)
BL=SQRT(DX*DX+DY*DY)
SI=DY/BL
CO=DX/BL
END
!Calculate element stifness matrix referred to element
!coordinate system
SUBROUTINE ESM (M,NE,E,A,ZI,BL,EK)
DIMENSION A(NE),ZI(NE),EK(6,6)
DO 10 I=1,6
DO 10 J=1,6
10 EK(I,J)=0.0
EA=E*A(M)
EI=E*ZI(M)
Ek(1,1)=EA/BL
Ek(1,4)=-EK(1,1)
Ek(2,2)=12.0*EI/BL/BL/BL
Ek(2,3)=6.0*EI/BL/BL
Ek(2,5)=-EK(2,2)
Ek(2,6)=EK(2,3)
Ek(3,3)=4.0*EI/BL
Ek(3,5)=-EK(2,3)
Ek(3,6)=2.0*EI/BL
Ek(4,4)=EK(1,1)
Ek(5,5)=EK(2,2)
Ek(5,6)=-EK(2,3)
Ek(6,6)=EK(3,3)
DO 20 J=2,6
DO 20 I=1,J-1
20 EK(J,I)=EK(I,J)
END
!Form coordinate transformation matrix
SUBROUTINE CTM(SI,CO,T)
DIMENSION T(6,6)
DO 10 I=1,6
DO 10 J=1,6
10 T(I,J)=0.0
   T(1,1)=CO
   T(1,2)=SI
   T(2,1)=-SI
   T(2,2)=CO
   T(3,3)=1.0
DO 20 I=1,3
DO 20 J=1,3
20 T(I+3,J+3)=T(I,J)
END
! Caculate element stiffness matrix referred
!to global coordinate system
SUBROUTINE TTKT(EK,T)
DIMENSION EK(6,6),T(6,6),TE(6,6)
DO 10 I=1,6
DO 10 J=1,6
TE(I,J)=0.0
DO 10 K=1,6
TE(I,J)=TE(I,J)+T(K,I)*EK(K,J)
10 CONTINUE
DO 20 I=1,6
DO 20 J=1,6
EK(I,J)=0.0
DO 20 K=1,6
EK(I,J)=EK(I,J)+TE(I,K)*T(K,J)
20 CONTINUE
END
!Assemble total stiffness matrix
SUBROUTINE TSM(NE,NJ,E,X,Y,IJ,A,ZI,TK,N)
DIMENSION X(NJ),Y(NJ),IJ(NE,2),A(NE),ZI(NE),TK(N,N),EK(6,6),T(6,6),LV(6)
DO 10 I=1,N
DO 10 J=1,N
10 TK(I,J)=0.0
DO 40 M=1,NE
CALL LSC(M,NE,NJ,X,Y,IJ,BL,SI,CO)
CALL ESM(M,NE,E,A,ZI,BL,EK)
CALL CTM(SI,CO,T)
CALL TTKT(EK,T)
DO 20 K=1,3
LV(K)=3*(IJ (M,1) -1)+K
LV(3+K)=3*(IJ (M,2) -1)+K
20 CONTINUE
DO 30 L=1,6
I=LV(L)
DO 30 K=1,6
J=LV(K)
TK(I,J)=TK(I,J)+EK(L,K)
30 CONTINUE
40 CONTINUE
   END
! Calculate element fixed-end forces.
SUBROUTINE EFF(L,PF,NF,BL,FO)
DIMENSION PF(NF,4),FO(6)
NO=INT(PF(L,2))
Q=PF(L,3)
C=PF(L,4)
B=BL-C
C1=C/BL
C2=C1*C1
C3=C1*C2
DO 5 I=1,6
5 FO(I)=0.0
GO TO(10,20,30,40,50,60),NO
10 FO(2)=-Q*C*(1.0-C2+C3/2.0)
   FO(3)=-Q*C*C*(0.5-2.0*C1/3.0+0.25*C2)
   FO(5)=-Q*C*C2*(1.0-0.5*C1)
   FO(6)=Q*C*C*C1*(1.0/3.0-0.25*C1)
   RETURN
20 FO(2)=-Q*B*B*(1.0+2.0*C1)/BL/BL
   FO(3)=-Q*C*B*B/BL/BL
   FO(5)=-Q*C2*(1.0+2.0*B/BL)
   FO(6)=Q*C2*B
   RETURN
30 FO(2)=6.0*Q*C1*B/BL/BL
   FO(3)=Q*B*(2.0-3.0*B/BL)/BL
   FO(5)=-6.0*Q*C1*B/BL/BL
   FO(6)=Q*C1*(2.0-3.0*C1)
   RETURN
40 FO(2)=-Q*C*(0.5-0.75*C2+0.4*C3)
   FO(3)=-Q*C*C*(1.0/3.0-0.5*C1+0.2*C2)
   FO(5)=-Q*C*C2*(0.75-0.4*C1)
   FO(6)=Q*C*C*C1*(0.25-0.2*C1)
   RETURN
50 FO(1)=-Q*C*(1.0-0.5*C1)
   FO(4)=-0.5*Q*C*C1
   RETURN
60   FO(1)=-Q*B/BL
   FO(4)=-Q*C1
   END
! Form total joint load vector
SUBROUTINE JLP (NE,NJ,NP,NF,X,Y,IJ,PJ,PF,P,N)
DIMENSION X(NJ),Y(NJ),IJ(NE,2),PJ(NP,3),PF(NF,4),P(N),FO(6),PE(6)
DO 10 I=1,N
P(I)=0.0
10 CONTINUE
IF(NP.GT.0)THEN
DO 20 I=1,NP
J=INT(PJ(I,1))
L=3*(J-1)+INT(PJ(I,2))
P(L)=PJ(I,3)
20 CONTINUE
ENDIF
IF(NF.GT.0)THEN
DO 30 L=1,NF
M=INT(PF(L,1))
CALL LSC (M,NE,NJ,X,Y,IJ,BL,SI,CO)
CALL EFF(L,PF,NF,BL,FO)
PE(1)=-FO(1)*CO+FO(2)*SI
PE(2)=-FO(1)*SI-FO(2)*CO
PE(3)=-FO(3)
PE(4)=-FO(4)*CO+FO(5)*SI
PE(5)=-FO(4)*SI-FO(5)*CO
PE(6)=-FO(6)
I=IJ(M,1)
J=IJ(M,2)
P(3*I-2)=P(3*I-2)+PE(1)
P(3*I-1)=P(3*I-1)+PE(2)
P(3*I)=P(3*I)+PE(3)
P(3*J-2)=P(3*J-2)+PE(4)
P(3*J-1)=P(3*J-1)+PE(5)
P(3*J)=P(3*J)+PE(6)
30 CONTINUE
ENDIF
END
! Introduce support condition
SUBROUTINE ISC(NR,JR,TK,P,N)
DIMENSION TK(N,N),P(N),JR(NR,4)
DO 30 I=1,NR
J=JR(I,1)
DO 20 K=1,3
IF(JR(I,K+1).NE.0)THEN
L=3*(J-1)+K
DO 10 JJ=1,N
TK(L,JJ)=0.0
TK(JJ,L)=0.0
10 CONTINUE
TK(L,L)=1.0
P(L)=0.0
END IF
20 CONTINUE
30 CONTINUE
END
! Solution of simultaneous equations by the GAUSS
! elimination method.
SUBROUTINE GAUSS(A,B,N)
DIMENSION A(N,N),B(N)
DO 20 K=1,N-1
DO 20 I=K+1,N
A1=A(K,I)/A(K,K)
DO 10 J=K+1,N
A(I,J)=A(I,J)-A1*A(K,J)
10 CONTINUE
B(I)=B(I)-A1*B(K)
20 CONTINUE
B(N)=B(N)/A(N,N)
DO 40 I=N-1,1,-1
DO 30 J=I+1,N
B(I)=B(I)-A(I,J)*B(J)
30 CONTINUE
B(I)=B(I)/A(I,I)
40 CONTINUE
END
! Print joint displacements.Calculate and print
!member-end forces of elements
SUBROUTINE MVN(NE,NJ,NF,E,X,Y,IJ,A,ZI,PF,P,N)
DIMENSION X(NJ),Y(NJ),IJ(NE,2),A(NE),ZI(NE),P(N),PF(NF,4),FO(6),F(6),D(6),TD(6),T(6,6),EK(6,6)
WRITE(2,10)
10 FORMAT(//2X,'JOINT DISPLACEMENTS'/5X,'JOINT',12X,'u',14X,'v',11X,'ceta')
DO 20 I=1,NJ
WRITE(2,15)I,P(3*I-2),P(3*I-1),P(3*I)
15 FORMAT  (2x,I6,4X,3E15.6)
20 CONTINUE
WRITE(2,25)
25 FORMAT(/2X,'MEMBER-END FORCES OF ELEMENTS'/4X,'ELEMENT',13X,'N',17X,'V',17X,'M')
DO 90 M=1,NE
CALL LSC(M,NE,NJ,X,Y,IJ,BL,SI,CO)
CALL ESM(M,NE,E,A,ZI,BL,EK)
CALL CTM(SI,CO,T)
I=IJ(M,1)
J=IJ(M,2)
DO 30 K=1,3
D(K)= P(3*(I-1)+K)
D(K+3)=P(3*(J-1)+K)
30 CONTINUE
DO 40 I=1,6
TD(I)=0.0
DO 40 J=1,6
TD(I)=TD(I)+T(I,J)*D(J)
40 CONTINUE
DO 50 I=1,6
F(I)=0.0
DO 50 J=1,6
F(I)=F(I)+EK(I,J)*TD(J)
50 CONTINUE
IF (NF.GT.0) THEN
DO 70 L=1,NF
I=INT(PF(L,1))
IF(M.EQ.I)THEN
CALL EFF (L,PF,NF,BL,FO)
DO 60 J=1,6
F(J)=F(J)+FO(J)
60 CONTINUE
END IF
70 CONTINUE
END IF
WRITE(2,80) M, (F(I),I=1,6)
80 FORMAT (2X,I8,4X,'N1=',F12.4,3X,'V1=',F12.4,3X,'M1=',F12.4/14X,'N2=',F12.4,3X,'V2=',F12.4,3X,'M2=',F12.4)
90 CONTINUE
END