*     橋台落下のMicroAVS描画データ作成
*     MicroAVSで読み込むfall.inpを出力
*                                          2005/7 伊津野

      parameter (nd=100)                   ! 最大節点数
      real x(nd),y(nd),z(nd)               ! 節点のX座標，Y座標，Z座標
      character moji*11                    ! ステップ数を表示する文字変数

*---- ここを変えるといろいろ変化可能 ------
      h=0.25                               ! 減衰定数
      open(1,file='fall.inp')              ! AVS用出力ファイル
      d=0.5                                ! 橋台の幅
      s=3.0                                ! スパン
      g=0.2                                ! 桁高さ
      p=3.0                                ! 橋台高さ
      a=0.2                                ! パラペットの幅
      nodes=68                             ! 節点数
      members=9                            ! 部材数
      ido=50                               ! 移動に要するステップ数
      dt=0.01                              ! 時間刻み
      nstep=301                            ! 総ステップ数
      am=10000.                            ! 質量
      ak=100000.                           ! ケーブル剛性
      ag=9.8                               ! 重力加速度
      dm=0.1                               ! ケーブル太さのdに対する割合
*------------------------------------------

      b=am*ag/2./ak/s
      omega=sqrt(3.*ak/am)
      hs=sqrt(1-h*h)

* 桁
      x( 1)= -d ; y( 1)=g  ; z( 1)=d
      x( 2)=s+d ; y( 2)=g  ; z( 2)=d
      x( 3)=s+d ; y( 3)=g  ; z( 3)=0.
      x( 4)= -d ; y( 4)=g  ; z( 4)=0.

      x( 5)= -d ; y( 5)=0. ; z( 5)=d
      x( 6)=s+d ; y( 6)=0. ; z( 6)=d
      x( 7)=s+d ; y( 7)=0. ; z( 7)=0.
      x( 8)= -d ; y( 8)=0. ; z( 8)=0.

* 左橋台1（固定）
      x( 9)= -d ; y( 9)=0. ; z( 9)=d
      x(10)= 0. ; y(10)=0. ; z(10)=d
      x(11)= 0. ; y(11)=0. ; z(11)=0.
      x(12)= -d ; y(12)=0. ; z(12)=0.

      x(13)= -d ; y(13)=-p ; z(13)=d
      x(14)= 0. ; y(14)=-p ; z(14)=d
      x(15)= 0. ; y(15)=-p ; z(15)=0.
      x(16)= -d ; y(16)=-p ; z(16)=0.

* 左橋台2（パラペット）
      x(17)= -a-d ; y(17)=g    ; z(17)=d
      x(18)=   -d ; y(18)=g    ; z(18)=d
      x(19)=   -d ; y(19)=g    ; z(19)=0.
      x(20)= -a-d ; y(20)=g    ; z(20)=0.

      x(21)= -a-d ; y(21)=-p   ; z(21)=d
      x(22)= -a-d ; y(22)=-p   ; z(22)=0

* 右橋台1（移動）
      x(23)=s   ; y(23)=0. ; z(23)=d
      x(24)=s+d ; y(24)=0. ; z(24)=d
      x(25)=s+d ; y(25)=0. ; z(25)=0.
      x(26)=s   ; y(26)=0. ; z(26)=0.

      x(27)=s   ; y(27)=-p ; z(27)=d
      x(28)=s+d ; y(28)=-p ; z(28)=d
      x(29)=s+d ; y(29)=-p ; z(29)=0.
      x(30)=s   ; y(30)=-p ; z(30)=0.

* 右橋台2（パラペット）
      x(31)= s+d   ; y(31)=g    ; z(31)=d
      x(32)= s+d+a ; y(32)=g    ; z(32)=d
      x(33)= s+d+a ; y(33)=g    ; z(33)=0.
      x(34)= s+d   ; y(34)=g    ; z(34)=0.

      x(35)= s+d+a ; y(35)=-p   ; z(35)=d
      x(36)= s+d+a ; y(36)=-p   ; z(36)=0.

* 落橋防止ケーブル2本
      x(37)= s+d ; y(37)=g*0.5      ; z(37)=d*(1.+dm)
      x(38)= s+d ; y(38)=g*(0.5+dm) ; z(38)=d*(1.+dm)
      x(39)= s+d ; y(39)=g*(0.5+dm) ; z(39)=d
      x(40)= s+d ; y(40)=g*0.50     ; z(40)=d

      x(41)= s+d ; y(41)=g*0.5      ; z(41)=-d*dm
      x(42)= s+d ; y(42)=g*(0.5+dm) ; z(42)=-d*dm
      x(43)= s+d ; y(43)=g*(0.5+dm) ; z(43)=0.
      x(44)= s+d ; y(44)=g*0.5      ; z(44)=0.

      x(45)= s+d ; y(45)=g*0.5      ; z(45)=d*(1.+dm)
      x(46)= s+d ; y(46)=g*(0.5+dm) ; z(46)=d*(1.+dm)
      x(47)= s+d ; y(47)=g*(0.5+dm) ; z(47)=d
      x(48)= s+d ; y(48)=g*0.5      ; z(48)=d

      x(49)= s+d ; y(49)=g*0.5      ; z(49)=-d*dm
      x(50)= s+d ; y(50)=g*(0.5+dm) ; z(50)=-d*dm
      x(51)= s+d ; y(51)=g*(0.5+dm) ; z(51)=0.
      x(52)= s+d ; y(52)=g*0.5      ; z(52)=0.

* 地盤1
      x(53)=-2.*(d+a)  ; y(53)=-p   ; z(53)=2.*d
      x(54)=s/2.       ; y(54)=-p   ; z(54)=2.*d
      x(55)=s/2.+d+a   ; y(55)=-p   ; z(55)=-d
      x(56)=-2.*(d+a)  ; y(56)=-p   ; z(56)=-d
      x(57)=-2.*(d+a)  ; y(57)=-p-d ; z(57)=2.*d
      x(58)=s/2.-d     ; y(58)=-p-d ; z(58)=2.*d
      x(59)=s/2.+a     ; y(59)=-p-d ; z(59)=-d
      x(60)=-2.*(d+a)  ; y(60)=-p-d ; z(60)=-d

* 地盤2
      x(61)=s/2.       ; y(61)=-p   ; z(61)=2.*d
      x(62)=s+3.*(d+a) ; y(62)=-p   ; z(62)=2.*d
      x(63)=s+3.*(d+a) ; y(63)=-p   ; z(63)=-d
      x(64)=s/2.+d+a   ; y(64)=-p   ; z(64)=-d
      x(65)=s/2.-d     ; y(65)=-p-d ; z(65)=2.*d
      x(66)=s+3.*(d+a) ; y(66)=-p-d ; z(66)=2.*d
      x(67)=s+3.*(d+a) ; y(67)=-p-d ; z(67)=-d
      x(68)=s/2.+a     ; y(68)=-p-d ; z(68)=-d

* 初期状態
      time=0.                              ! 表示用の時間
      write(1,'(a9)') '# Falling'          ! ここからAVSファイル出力
      write(1,*) nstep
      write(1,'(a4)') 'geom'
      write(1,'(a5,f6.2,a1)') 'step1',time,'s'   ! コメントとして時間を使う
      write(1,*) nodes, members
* 節点番号と座標値
      do i=1,nodes
        write(1,*) i,x(i),y(i),z(i)
      end do
* 部材データ
      write(1,'(a31)') '1 0 hex  1  2  3  4  5  6  7  8'   ! 桁
      write(1,'(a31)') '2 0 hex  9 10 11 12 13 14 15 16'   ! 左橋台
      write(1,'(a31)') '3 0 hex 17 18 19 20 21 13 16 22'
      write(1,'(a31)') '4 0 hex 23 24 25 26 27 28 29 30'   ! 右橋台
      write(1,'(a31)') '5 0 hex 31 32 33 34 28 35 36 29'
      write(1,'(a31)') '6 0 hex 37 38 39 40 45 46 47 48'   ! ケーブル1
      write(1,'(a31)') '7 0 hex 41 42 43 44 49 50 51 52'   ! ケーブル2
      write(1,'(a31)') '8 0 hex 53 54 55 56 57 58 59 60'   ! 地盤1
      write(1,'(a31)') '9 0 hex 61 62 63 64 65 66 67 68'   ! 地盤2
* 節点データ
      write(1,*) '1 0'
      write(1,*) '1 1'
      write(1,*) 'material, No.' ! 部材によって色を変える
      do i=1,8
        write(1,*) i, 2.         ! 桁
      end do
      do i=9,36
        write(1,*) i, 3.         ! 橋台
      end do
      do i=37,52
        write(1,*) i, 4.         ! ケーブル
      end do
      do i=53,56
        write(1,*) i, 1          ! 地盤
      end do
      do i=57,60
        write(1,*) i, 0          ! 地下
      end do
      do i=61,64
        write(1,*) i, 1          ! 地盤
      end do
      do i=65,68
        write(1,*) i, 0          ! 地下
      end do

* 移動                           ! 最初の50ステップで右橋台を移動
      dy=d/float(ido)
      do istep=2,ido+1
        do i=23,44
          x(i)=x(i)+dy
        end do
        x(58)=x(58)+dy
        x(59)=x(59)+dy
        x(61)=x(61)+dy
        x(64)=x(64)+dy
        x(65)=x(58)
        x(68)=x(59)

        write(moji,*) istep
        moji='step'//moji(2:)                  ! ステップ番号の数字を文字に変換
        write(1,'(a11,f6.2,a1)') moji,time,'s'
        write(1,*) nodes, members
        do i=1,nodes
          write(1,*) i,x(i),y(i),z(i)
        end do
        write(1,'(a31)') '1 0 hex  1  2  3  4  5  6  7  8'   ! 桁
        write(1,'(a31)') '2 0 hex  9 10 11 12 13 14 15 16'   ! 左橋台
        write(1,'(a31)') '3 0 hex 17 18 19 20 21 13 16 22'
        write(1,'(a31)') '4 0 hex 23 24 25 26 27 28 29 30'   ! 右橋台
        write(1,'(a31)') '5 0 hex 31 32 33 34 28 35 36 29'
        write(1,'(a31)') '6 0 hex 37 38 39 40 45 46 47 48'   ! ケーブル1
        write(1,'(a31)') '7 0 hex 41 42 43 44 49 50 51 52'   ! ケーブル2
        write(1,'(a31)') '8 0 hex 53 54 55 56 57 58 59 60'   ! 地盤1
        write(1,'(a31)') '9 0 hex 61 62 63 64 65 66 67 68'   ! 地盤2

      end do

* 桁落下
      do istep=ido+2,nstep
        if (istep.le.60) then      ! 落下開始後10ステップは右橋台をさらに移動
          do i=23,44
            x(i)=x(i)+dy
          end do
          x(58)=x(58)+dy
          x(59)=x(59)+dy
          x(61)=x(61)+dy
          x(64)=x(64)+dy
          x(65)=x(58)
          x(68)=x(59)
        end if

        time=time+dt
        theta=b*(1.-exp(-h*omega*time)*cos(hs*omega*time))   ! 落下角度の計算

        x(6)=(s+d)*cos(theta)      ; y(6)=-(s+d)*sin(theta)
        x(2)=x(6)+g*sin(theta)     ; y(2)=y(6)+g*cos(theta)
        x(5)=-d*cos(theta)         ; y(5)=d*sin(theta)
        x(1)=x(5)+g*sin(theta)     ; y(1)=y(5)+g*cos(theta)
        x(4)=x(1)                  ; y(4)=y(1)
        x(3)=x(2)                  ; y(3)=y(2)
        x(7)=x(6)                  ; y(7)=y(6)
        x(8)=x(5)                  ; y(8)=y(5)
        x(45)=x(6)+g*sin(theta)*0.5 ; y(45)=y(6)+g*cos(theta)*0.5
        x(46)=x(6)+g*sin(theta)*(0.5+dm)
        y(46)=y(6)+g*cos(theta)*(0.5+dm)
        x(47)=x(46) ; y(47)=y(46) ; x(48)=x(45) ; y(48)=y(45)
        x(50)=x(46) ; y(50)=y(46) ; x(49)=x(45) ; y(49)=y(45)
        x(51)=x(46) ; y(51)=y(46) ; x(52)=x(45) ; y(52)=y(45)

        write(moji,*) istep
        moji='step'//moji(2:)                  ! ステップ番号の数字を文字に変換
        write(1,'(a11,f6.2,a1)') moji,time,'s'
        write(1,*) nodes, members
        do i=1,nodes
          write(1,*) i,x(i),y(i),z(i)
        end do
        write(1,'(a31)') '1 0 hex  1  2  3  4  5  6  7  8'   ! 桁
        write(1,'(a31)') '2 0 hex  9 10 11 12 13 14 15 16'   ! 左橋台
        write(1,'(a31)') '3 0 hex 17 18 19 20 21 13 16 22'
        write(1,'(a31)') '4 0 hex 23 24 25 26 27 28 29 30'   ! 右橋台
        write(1,'(a31)') '5 0 hex 31 32 33 34 28 35 36 29'
        write(1,'(a31)') '6 0 hex 37 38 39 40 45 46 47 48'   ! ケーブル1
        write(1,'(a31)') '7 0 hex 41 42 43 44 49 50 51 52'   ! ケーブル2
        write(1,'(a31)') '8 0 hex 53 54 55 56 57 58 59 60'   ! 地盤1
        write(1,'(a31)') '9 0 hex 61 62 63 64 65 66 67 68'   ! 地盤2

      end do

      close(1)
      end