C***********************************************************************
C     CABLE-STAYED BRIDGE FRAM ANALYSIS PROGRAM
C                                                 BY ATSUSHI ITOH, 2002                
C***********************************************************************
C
	IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(IE=100,IS=100)
C
      DIMENSION X(IS),Y(IS),IES(2,IE)
      DIMENSION ITPCNP(IS*3),ITBCNP(IS*3),IADR(IS*3)
      DIMENSION TDK(IS*3,IS*3),DU(IS*3),DF(IS*3),FF(IS*3)
      DIMENSION EA(IE),EI(IE),DNQM(6,IE),DISLO(IE)
C
      OPEN(1,FILE='input.dat')          !　データファイル
	OPEN(2,FILE='output.dat')         !　出力ファイル
C
C----- 荷重：（tf），　　　　変位：（m）
C----- IE：要素数，          IS：節点数 DIMENSIONを決定するパラメータ
C----- X,Y：座標データ，     IES(2,IE)：要素を構成する節点番号
C----- TDK：剛性マトリクス ，DU：節点変位増分ベクトル
C----- DF：反力増分ベクトル，EA,EI：断面剛性
C----- FF：外力荷重ベクトル，DNQM：断面力，DISLO：当分布荷重ベクトル
C----- ITPCNP,ITBCNP,IADR：剛性マトリクス等の並び替えに必要な情報
C
C----- １．データ入力 -----
      CALL DATAIN(X,Y,ITPCNP,ITBCNP,IADR,EA,EI,IES,IEL,INODE,IPN,IBN)
C
C----- ２．外力荷重ベクトルの設定 {FF} -----
      CALL LOAD(INODE,IEL,FF,IES,X,Y,DISLO)
C
C----- ３．剛性マトリクスの作成 [TDK] -----
      CALL MKMATRIX(TDK,X,Y,IADR,IES,IEL,EA,EI,INODE)
C
C----- ４．剛性方程式を解く {DF}=[TDK]{DU} -----
      CALL CALMAT(TDK,FF,DU,ITBCNP,ITPCNP,IBN,IPN,INODE,DF)
C
C----- ５．断面力の算出 DNQM{6,IE} = N1,Q1,M1,N2,Q2,M2 -----
      CALL SECTION(DU,IEL,X,Y,IES,EA,EI,DNQM,DISLO)
C
C----- ６．計算結果の出力 -----
      CALL OUTPUT(DU,DF,DNQM,IEL,INODE)
C
      WRITE(6,*) 'Press any key to continue'
      READ(5,*)
      STOP
      END
C***********************************************************************
      SUBROUTINE DATAIN(X,Y,ITPCNP,ITBCNP,IADR,EA,EI,IES,IEL,INODE,
     -                  IPN,IBN)
C***********************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(IE=100,IS=100)
C
      DIMENSION X(IS),Y(IS),ITPCNP(IS*3),ITBCNP(IS*3),IADR(IS*3)
      DIMENSION EA(IE),EI(IE),IES(2,IE)
      DIMENSION IBCNP(3,IS)
C
C----- INPUT INODE = 総節点数
      READ(1,100) INODE
      WRITE(*,100) INODE
  100 FORMAT(I10)
C
C----- INPUT X,Y = 節点座標，IBCNP = 幾何学的境界条件 1:固定，0:自由
      DO 10 I1=1,INODE
      READ(1,110) X(I1),Y(I1),(IBCNP(J,I1),J=1,3)
      WRITE(*,110) X(I1),Y(I1),(IBCNP(J,I1),J=1,3)
   10 CONTINUE
  110 FORMAT(2F10.3,4X,3I2)
C
C----- INPUT IEL = 要素数
      READ(1,120) IEL
      WRITE(*,120) IEL
  120 FORMAT(I10)
C
C----- INPUT IES = 要素を構成する節点番号，EA,EI = 断面剛性
      DO 20 I1=1,IEL
      READ(1,130)  IES(1,I1),IES(2,I1),DA,DI,E
      WRITE(*,130) IES(1,I1),IES(2,I1),DA,DI,E
      EA(I1)=E*DA
      EI(I1)=E*DI
   20 CONTINUE
  130 FORMAT(2I5,3D10.3)
C
C----- ITPCNP：並び替えた自由なその自由度番号をもとに並べた時の位置
C----- ITBCNP：並び替えた拘束なその自由度番号をもとに並べた時の位置
C----- IPN：自由自由度の総数，IBN：拘束自由度の総数
      MN=0                                                              
      NN=0                                                              
      KK=0                                                              
      DO 30 I=1,INODE
      DO 30 J=1,3                                                       
      KK=KK+1                                                           
      IF(IBCNP(J,I).EQ.1) THEN                                          
      MN=MN+1                                                           
      ITBCNP(MN)=KK                                                     
      ENDIF                                                             
      IF(IBCNP(J,I).EQ.0)THEN                                           
      NN=NN+1                                                           
      ITPCNP(NN)=KK                                                     
      ENDIF                                                             
   30 CONTINUE                                                          
      IBN=MN                                                            
      IPN=NN
C----- 順番に並んだ自由度番号が，剛性マトリクスの何処に格納されるか
      M=1
      DO 40 I=1,INODE*3
      IF(I.EQ.ITPCNP(M)) THEN
      IADR(I)=M
      M=M+1
      ENDIF
   40 CONTINUE
      DO 50 I=1,INODE*3
      IF(I.EQ.ITBCNP(M-IPN)) THEN
      IADR(I)=M
      M=M+1
      ENDIF
   50 CONTINUE
C
      RETURN
      END
C***********************************************************************
      SUBROUTINE LOAD(INODE,IEL,FF,IES,X,Y,DISLO)
C***********************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(IE=100,IS=100)
      CHARACTER*30 AAAA,BBBB
C
      DIMENSION FF(IS*3),IES(2,IE),X(IS),Y(IS),DISLO(IE)
C
      DO 10 I=1,INODE*3
   10  FF(I)=0.0D0
C
C----- 集中荷重の入力 -----
C----- INPUT IDIR = 荷重の方向，1 : X 方向，2 : Y 方向，3 : モーメント
C----- INPUT INODE1 = 節点番号，DPPP = その節点に作用させる荷重ベクトル
C----- IDIR = 0 とすると，集中荷重の入力終了
      READ(1,*) AAAA                !　使用しない
      WRITE(*,*) AAAA
  100 READ(1,555) INODE1,IDIR,DP
      WRITE(*,555)  INODE1,IDIR,DP
      IF(IDIR.EQ.0) GOTO 200
  555 FORMAT(I5,I5,D10.3)
      IF(IDIR.EQ.1) INO=INODE1*3-2
      IF(IDIR.EQ.2) INO=INODE1*3-1
      IF(IDIR.EQ.3) INO=INODE1*3
      FF(INO)=DP
      GOTO 100
  200 CONTINUE
C
C----- 分布荷重の入力 -----
C----- INPUT IELEM = 要素番号，DPP = 等分布荷重（tonf/m)
C----- IELEM = 0 とすると，分布荷重の入力終了
      DO 250 I1=1,IEL
  250  DISLO(I1)=0.0D0
      READ(1,*) BBBB                !　使用しない
      WRITE(*,*) BBBB
  300 READ(1,999) IELEM,DPP
      WRITE(*,999) IELEM,DPP
      IF(IELEM.EQ.0) GOTO 400
       DISLO(IELEM)=DPP
  999 FORMAT(I5,D15.3)
C
      ELL = ( (X(IES(2,IELEM))- X(IES(1,IELEM)))**2 +
     -        (Y(IES(2,IELEM))- Y(IES(1,IELEM)))**2 )**0.5D0
      SIN = (Y(IES(2,IELEM))- Y(IES(1,IELEM)))/ELL
      COS = (X(IES(2,IELEM))- X(IES(1,IELEM)))/ELL
      WH = DPP * SIN
      WV = DPP * COS
C
      FF(IES(1,IELEM)*3-2)=FF(IES(1,IELEM)*3-2)+WH*ELL/2.0D0
      FF(IES(1,IELEM)*3-1)=FF(IES(1,IELEM)*3-1)+WV*ELL/2.0D0
      FF(IES(1,IELEM)*3  )=FF(IES(1,IELEM)*3  )+WV*ELL**2/12.0D0
      FF(IES(2,IELEM)*3-2)=FF(IES(2,IELEM)*3-2)+WH*ELL/2.0D0
      FF(IES(2,IELEM)*3-1)=FF(IES(2,IELEM)*3-1)+WV*ELL/2.0D0
      FF(IES(2,IELEM)*3  )=FF(IES(2,IELEM)*3  )-WV*ELL**2/12.0D0
      GOTO 300
  400 CONTINUE
C
      RETURN
      END
C***********************************************************************
      SUBROUTINE MKMATRIX(TDK,X,Y,IADR,IES,IEL,EA,EI,INODE)
C***********************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(IE=100,IS=100)
C
      DIMENSION TDK(IS*3,IS*3),IADR(IS*3),IES(2,IE),EA(IE),EI(IE)
      DIMENSION TKL(6,6),TT(6,6),DTEK(6,6),TEK(6,6),IA(6),X(IS),Y(IS)
C
      DO 5 I1=1,INODE*3
      DO 5 I2=1,INODE*3
    5 TDK(I1,I2)=0.0D0
C
      DO 999 I=1,IEL
C
      CALL KLOCAL(I,X,Y,IES,EA,EI,TKL)
      CALL  TRANS(I,X,Y,TT,IES)
C
C----- MAKE DTEK[6,6] MATRIX [TT]t[K]local -----
      DO 50 I1=1,6
      DO 50 J1=1,6
      DTEK(I1,J1)=0.0D0
      DO 50  K=1,6
   50 DTEK(I1,J1)=DTEK(I1,J1)+TT(K,I1)*TKL(K,J1)
C----- MAKE TEK[6,6] MATRIX [DTEK][TT] -----
      DO 60 I1=1,6
      DO 60 J1=1,6
      TEK(I1,J1)=0.0D0
      DO 60  K=1,6
   60 TEK(I1,J1)=TEK(I1,J1)+DTEK(I1,K)*TT(K,J1)
C----- ADD TO STIFFNESS MATRIX TDK[INODE*3,INODE*3] -----------------
      IA(1)=3*IES(1,I)-2
      IA(2)=3*IES(1,I)-1
      IA(3)=3*IES(1,I)
      IA(4)=3*IES(2,I)-2
      IA(5)=3*IES(2,I)-1
      IA(6)=3*IES(2,I)
C
      DO 100 I1=1,6
      DO 100 I2=1,6
      TDK(IADR(IA(I1)),IADR(IA(I2)))=TDK(IADR(IA(I1)),IADR(IA(I2)))
     *                              +TEK(I1,I2)
  100 CONTINUE
C
  999 CONTINUE
C
      RETURN
      END
C***********************************************************************
      SUBROUTINE KLOCAL(I,X,Y,IES,EA,EI,SM)
C***********************************************************************
      IMPLICIT REAL*8(A-H,O-Z)                                          
      PARAMETER(IE=100,IS=100)
C
      DIMENSION X(IS),Y(IS),IES(2,IE),EA(IE),EI(IE),SM(6,6)
C
C----- 局所座標系における要素剛性を求める
      DO 10 I1=1,6
      DO 10 J1=1,6
   10 SM(I1,J1)=0.0D0
C
      XL2=(X(IES(2,I))-X(IES(1,I)))**2
      YL2=(Y(IES(2,I))-Y(IES(1,I)))**2
      ELL=DSQRT(XL2+YL2)
C
      S0=     EA(I)/ELL
      S1=  2.*EI(I)/ELL
      S2=  6.*EI(I)/ELL**2
      S3= 12.*EI(I)/ELL**3
C
      SM(1,1)=  S0
      SM(1,4)= -S0
      SM(2,2)=  S3
      SM(2,3)=  S2
      SM(2,5)= -S3
      SM(2,6)=  S2
      SM(3,3)=  S1*2
      SM(3,5)= -S2
      SM(3,6)=  S1
      SM(4,4)=  S0
      SM(5,5)=  S3
      SM(5,6)= -S2
      SM(6,6)=  S1*2
C
      DO 20 I1=2,6
      DO 20 J1=1, I1-1
   20  SM(I1,J1)=SM(J1,I1)
C
      RETURN
      END
C***********************************************************************
      SUBROUTINE TRANS(I,X,Y,TT,IES)
C***********************************************************************
      IMPLICIT REAL*8(A-H,O-Z)                                          
      PARAMETER(IE=100,IS=100)
C
      DIMENSION IES(2,IE),X(IS),Y(IS),TT(6,6)
C
C----- MAKE TC[6,6] MATRIX　座標変換マトリクスを作る
      XL2=(X(IES(2,I))-X(IES(1,I)))**2
      YL2=(Y(IES(2,I))-Y(IES(1,I)))**2
      ELL=DSQRT(XL2+YL2)
      S=(Y(IES(2,I))-Y(IES(1,I)))/ELL
      C=(X(IES(2,I))-X(IES(1,I)))/ELL
C
      DO 30 I1=1,6
      DO 30 J1=1,6
   30 TT(I1,J1)=0.0D0
C
      TT(1,1)=C
      TT(1,2)=S
      TT(2,1)=-S
      TT(2,2)=C
      TT(3,3)=1.0D0
      DO 35 I1=1,3
      DO 35 J1=1,3
   35 TT(I1+3,J1+3)=TT(I1,J1)
C
      RETURN
      END
C***********************************************************************
      SUBROUTINE CALMAT(TDK,FF,DU,ITBCNP,ITPCNP,IBN,IPN,INODE,DF)
C***********************************************************************                                      
      IMPLICIT REAL*8(A-H,O-Z) 
      LOGICAL IERR                           
      PARAMETER(IS=100)
      DIMENSION TDK(IS*3,IS*3),FF(IS*3),DU(IS*3),DF(IS*3)
      DIMENSION ITPCNP(IS*3),ITBCNP(IS*3),F1(IS*3),DFu2(IS*3)
      DIMENSION DUk(IS*3),DUu(IS*3),DFk(IS*3),DFu(IS*3),DFu1(IS*3)
      DIMENSION TK11(IS*3,IS*3),DQ(IS*3)
C
      DO 10 I1=1,INODE*3
       DU(I1)=0.0D0                 ! 節点変位
       DQ(I1)=0.0D0                 ! 等価節点力
   10  DF(I1)=0.0D0                 ! 反力
C----- 既知変位成分は全て　ゼロ（ただし，変位を制御する場合もある）
      DO 15 I1=1,IBN
   15 DUk(I1)=0.0D0
C----- {F1}=[K12]{Xk} ----- ゼロの場合は必要ない
      DO 20 I1=1,IPN
       F1(I1)=0.0D0
      DO 20 J1=1,IBN
   20  F1(I1)=F1(I1)+TDK(I1,IPN+J1)*DUk(J1)
C----- 剛性マトリクスの保存（逆マトリクスを計算するため）
      DO 30 I1=1,IPN
      DO 30 J1=1,IPN
       TK11(I1,J1)=TDK(I1,J1)
   30 CONTINUE
C----- 逆マトリクスを求める
      CALL GJINV(TK11,IS*3,IERR)
C----- {Xu}=[K11]^(-1){{Fk}-{F1}} 未知変位ベクトルの計算
C----- 剛性マトリクスは並び替えているが，{FF}は並び替えていないので
C----- {ITOCNP}を使用して{FF}を並び替えている
      DO 40 I1=1,IPN
      DUu(I1)=0.0D0
      DO 40 J1=1,IPN
      DUu(I1)=DUu(I1)+TK11(I1,J1)*(FF(ITPCNP(J1))-F1(J1))
   40 CONTINUE
C----------------------------------------------------------------------
C----- 既知外力{Fk}を保存
      DO 55 I1=1,IPN
   55  DFk(I1)=FF(ITPCNP(I1))
C----- {DFu1}=[K21]{Xu} 
      DO 60 I1=1,IBN
       DFu1(I1)=0.0D0
      DO 60 J1=1,IPN
   60  DFu1(I1)=DFu1(I1)+TDK(IPN+I1,J1)*DUu(J1)
C----- {DFu2}=[K22]{Xk} （{Xk}={0}であれば必要ない） 
      DO 70 I1=1,IBN
       DFu2(I1)=0.0D0
      DO 70 J1=1,IBN
   70  DFu2(I1)=DFu2(I1)+TDK(IPN+I1,IPN+J1)*DUk(J1)
C----- {未知外力}={DFu1}+{DFu2}
      DO 80 I1=1,IBN
   80  DFu(I1)=DFu1(I1)+DFu2(I1)
C----- FOR EVALUATE EQUIVALENT NODAL FORCE and DISPLACEMENT -----
C----- 以上で等価節点力が算出，全てを元に並び替える
      DO 90 I1=1,IPN
       DQ(ITPCNP(I1))=DFk(I1)
   90  DU(ITPCNP(I1))=DUu(I1)
C
      DO 100 I1=1,IBN
       DQ(ITBCNP(I1))=DFu(I1)
  100  DU(ITBCNP(I1))=DUk(I1)
C----- FOR EVALUATE REACTION FORCE -----
C----- 反力の計算 -----
      DO 110 I1=1,INODE*3
  110 DF(I1)=DQ(I1)-FF(I1)
C
      RETURN
      END
C***********************************************************************
      SUBROUTINE GJINV(A,N,IERR)
C***********************************************************************
      REAL*8 A(N,N)
      INTEGER N
      LOGICAL IERR
C
C----- [A(N,N)]の逆行列を求めるサブルーチン
      IERR= .FALSE.
      DO 990 K=1, N
         S= A(K, K)
         IF(ABS(S).LT.1E-6) GO TO 999
         DO 991 J = 1, N
            A(K,J) = A(K,J)/S
  991    CONTINUE
         A(K,K) = 1.0/S
         DO 992 I=1, N
            IF(I.EQ.K) GO TO 992
            S= A(I,K)
            DO 993 J=1, N
               A(I,J) = A(I,J) -S*A(K,J)
  993       CONTINUE
            A(I,K) = -S*A(K,K)
  992    CONTINUE
  990 CONTINUE
      RETURN
  999 CONTINUE
         IERR= .TRUE.
      RETURN
      END
C***********************************************************************
      SUBROUTINE SECTION(DU,IEL,X,Y,IES,EA,EI,DNQM,DISLO)
C***********************************************************************
      IMPLICIT REAL*8(A-H,O-Z)                                          
      PARAMETER(IE=100,IS=100)
C
      DIMENSION DU(IS*3),X(IS),Y(IS),IES(2,IE),EA(IE),EI(IE),F(6)
      DIMENSION TT(6,6),DNQM(6,IE),DUL(6),DUL1(6),TKL(6,6),DISLO(IE)
C----- 断面力を算出する
      DO 999 I=1,IEL
C
      CALL KLOCAL(I,X,Y,IES,EA,EI,TKL)
      CALL  TRANS(I,X,Y,TT,IES)
C----- その要素の節点の変位成分を取り出す
      DUL1(1)=DU(IES(1,I)*3-2)
      DUL1(2)=DU(IES(1,I)*3-1)
      DUL1(3)=DU(IES(1,I)*3  )
      DUL1(4)=DU(IES(2,I)*3-2)
      DUL1(5)=DU(IES(2,I)*3-1)
      DUL1(6)=DU(IES(2,I)*3  )
C----- 取り出した全体系の変位を局所系に座標変換する
      DO 100 I1=1,6
       DUL(I1)=0.0D0
      DO 100 J1=1,6
  100  DUL(I1)=DUL(I1)+TT(I1,J1)*DUL1(J1)
C----- 断面力の算出 DNQM{6,要素}=[K 局所系]{X 局所系}
      DO 200 I1=1,6
       DNQM(I1,I)=0.0D0
      DO 200 J1=1,6
  200  DNQM(I1,I)=DNQM(I1,I)+TKL(I1,J1)*DUL(J1)
C----- 要素の長さ，角度を求める
C----- {DISLO}：要素に作用する分布荷重（INPUT DATA値），Sub.LOAD 参照
      ELL = ( (X(IES(2,I))- X(IES(1,I)))**2 +
     -        (Y(IES(2,I))- Y(IES(1,I)))**2 )**0.5D0
      SIN = (Y(IES(2,I))- Y(IES(1,I)))/ELL
      COS = (X(IES(2,I))- X(IES(1,I)))/ELL
      WH = DISLO(I) * SIN
      WV = DISLO(I) * COS
C----- 等分布荷重による部材固定端力
      F(1)=+WH*ELL/2.0D0
      F(2)=+WV*ELL/2.0D0
      F(3)=+WV*ELL**2/12.0D0
      F(4)=+WH*ELL/2.0D0
      F(5)=+WV*ELL/2.0D0
      F(6)=-WV*ELL**2/12.0D0
C----- このプログラムでは鉛直上方向を正としている
C----- 変形による力と等分布荷重による部材固定端力の和
      DO 300 J1=1,6
  300  DNQM(J1,I)=DNQM(J1,I)+(-F(J1))
C
  999 CONTINUE
C
      RETURN
      END
C***********************************************************************
      SUBROUTINE OUTPUT(DU,DF,DNQM,IEL,INODE)
C***********************************************************************
      IMPLICIT REAL*8(A-H,O-Z)                                          
      PARAMETER(IE=100,IS=100)
C
      DIMENSION DU(IS*3),DF(IS*3),DNQM(6,IE)
C
C----- 節点力の出力 -----
C
      WRITE(2,*)
      WRITE(2,*) '節点変位 '
      WRITE(2,*)"     NO.   u(m)        v(m)        r(rad)   "
      DO 10 I1=1,INODE
   10 WRITE(2,100) I1,DU(3*I1-2),DU(3*I1-1),DU(3*I1)
  100 FORMAT(6X,I2,3F12.8)
C
C----- 支点反力の出力 -----
C
      WRITE(2,*)
      WRITE(2,*) '支点反力 '
      WRITE(2,*)"     NO.     Ru(tf)      Rv(tf)      M(tf-m)   "
      DO 20 I1=1,INODE
   20 WRITE(2,200) I1,DF(3*I1-2),DF(3*I1-1),DF(3*I1)
  200 FORMAT(6X,I2,3F12.4)
C
C----- 断面力の出力 -----
C
      WRITE(2,*)
      WRITE(2,*) '断面力 '
      WRITE(2,300)  "NO.","N1(tf)","N2(tf)","S1(tf)","S2(tf)",
     &               "M1(tfm)","M2(tfm)"
  300 FORMAT(1X,A8,1X,4A8,2A9)
C
      DO 30 I=1,IEL
      WRITE(2,400) I, DNQM(1,I), -DNQM(4,I), DNQM(2,I), -DNQM(5,I),
     &               -DNQM(3,I),  DNQM(6,I)
   30 CONTINUE
  400 FORMAT(1X,I8,1X,4F8.2,2F9.2)
C
      RETURN
      END

