*     橋脚地震応答のMicroAVS描画データ作成
*     2構造物の比較バージョン
*     resp.inp というMicroAVSファイルを出力
*                             2005/7 伊津野
      parameter (nd=300)                   ! 最大節点数
      real x(nd),y(nd),z(nd)               ! 左構造物節点のX座標，Y座標，Z座標
      real u(nd),v(nd),w(nd)               ! 右構造物の座標
      character moji*11                    ! ステップ数を表示する文字変数
      character fname*80                   ! ファイル名記憶用

      d=0.25                               ! 橋脚の幅の半分
      s=3.0                                ! スパン
      h=0.1                                ! 支承高さ
      g=0.2                                ! 桁高さ
      p=3.0                                ! 橋脚高さ
      dp=0.2                               ! 橋脚分割高さ
      nn=15                                ! 分割数
      nodes=8*(nn+4)                       ! 節点数
      members=2*nn+3                       ! 部材数
      a=4.                                 ! 2つの構造物間距離
      cs=0.05 ; cg=0.02                    ! 支承と桁の色

* 各ファイルに橋脚天端と，支承（桁−橋脚相対変位）の応答変位が
* 1つ以上の空白で区切られたものを準備．

      write(*,*) 'Time Increment : '       ! 時間刻み
      read (*,*) dt
      write(*,*) 'Input file name for 1'
      read(*,'(a80)') fname
      open(2,file=fname,status='old')      ! 応答結果ファイル
      write(*,*) 'Input file name for 2'
      read(*,'(a80)') fname
      open(3,file=fname,status='old')      ! 応答結果ファイル

* 橋脚
      zz=p
      m=nn*4
      do i=1,m+1,4
        x(i    )=-d ; y(i    )= -d ; z(i    )=zz
        x(i+  1)= d ; y(i+  1)= -d ; z(i+  1)=zz
        x(i+  2)= d ; y(i+  2)=  d ; z(i+  2)=zz
        x(i+  3)=-d ; y(i+  3)=  d ; z(i+  3)=zz
        x(i+m+4)=-d ; y(i+m+4)=s-d ; z(i+m+4)=zz
        x(i+m+5)= d ; y(i+m+5)=s-d ; z(i+m+5)=zz
        x(i+m+6)= d ; y(i+m+6)=s+d ; z(i+m+6)=zz
        x(i+m+7)=-d ; y(i+m+7)=s+d ; z(i+m+7)=zz
        zz=zz-dp
      end do
* 桁
      k=(nn+1)*8
      x(k+1)=-d ; y(k+1)= -d ; z(k+1)=p+h+g
      x(k+2)= d ; y(k+2)= -d ; z(k+2)=p+h+g
      x(k+3)= d ; y(k+3)=s+d ; z(k+3)=p+h+g
      x(k+4)=-d ; y(k+4)=s+d ; z(k+4)=p+h+g
      x(k+5)=-d ; y(k+5)= -d ; z(k+5)=p+h
      x(k+6)= d ; y(k+6)= -d ; z(k+6)=p+h
      x(k+7)= d ; y(k+7)=s+d ; z(k+7)=p+h
      x(k+8)=-d ; y(k+8)=s+d ; z(k+8)=p+h
* 支承
      ks=k+9; ke=k+17
      yy=0.
      do j=ks,ke,8
        zz=p+h
        do i=j,j+4,4
          x(i  )=-d/2. ; y(i  )=-d/2.+yy ; z(i  )=zz
          x(i+1)= d/2. ; y(i+1)=-d/2.+yy ; z(i+1)=zz
          x(i+2)= d/2. ; y(i+2)= d/2.+yy ; z(i+2)=zz
          x(i+3)=-d/2. ; y(i+3)= d/2.+yy ; z(i+3)=zz
          zz=zz-h
        end do
        yy=yy+s
      end do
      do i=1,nodes
        u(i)=x(i)
        v(i)=y(i)
        w(i)=z(i)
      end do

      time=0.                              ! 表示用の時間
      nstep=501                            ! ステップ数

      open(1,file='resp.inp')              ! AVS用出力ファイル

      write(1,'(a21)') '# Earthquake Response'   ! ここからAVSファイル出力
      write(1,*) nstep
      write(1,'(a9)') 'data_geom'
      write(1,'(a5,f6.2,a1)') 'step1',time,'s'   ! コメントとして時間を使う
      write(1,*) 2*nodes, 2*members
* 節点番号と座標値
      do i=1,nodes
        write(1,*) i,x(i),y(i),z(i)
      end do
      do i=1,nodes
        write(1,*) nodes+i,u(i),v(i)+a,w(i)
      end do
* 部材データ   左橋脚
      j=1
      do i=1,m,4
        write(1,'(i4,a6,8i4)') j,' 0 hex',i,i+1,i+2,i+3,i+4,i+5,i+6,i+7
        j=j+1
      end do
* 右橋脚
      do i=m+5,k-7,4
        write(1,'(i4,a6,8i4)') j,' 0 hex',i,i+1,i+2,i+3,i+4,i+5,i+6,i+7
        j=j+1
      end do
* 桁
      write(1,'(i4,a6,8i4)') j,' 0 hex',
     :          k+1,k+2,k+3,k+4,k+5,k+6,k+7,k+8
* 左支承
      write(1,'(i4,a6,8i4)') j+1,' 0 hex',
     :          ks,ks+1,ks+2,ks+3,ks+4,ks+5,ks+6,ks+7
* 右支承
      write(1,'(i4,a6,8i4)') j+2,' 0 hex',
     :          ks+8,ks+9,ks+10,ks+11,ks+12,ks+13,ks+14,ks+15

* 構造物2
      j=j+3
      do i=nodes+1,nodes+m,4
        write(1,'(i4,a6,8i4)') j,' 0 hex',i,i+1,i+2,i+3,i+4,i+5,i+6,i+7
        j=j+1
      end do
      do i=nodes+m+5,nodes+k-7,4
        write(1,'(i4,a6,8i4)') j,' 0 hex',i,i+1,i+2,i+3,i+4,i+5,i+6,i+7
        j=j+1
      end do
      write(1,'(i4,a6,8i4)') j,' 0 hex',
     :          nodes+k+1,nodes+k+2,nodes+k+3,nodes+k+4,
     :          nodes+k+5,nodes+k+6,nodes+k+7,nodes+k+8
      write(1,'(i4,a6,8i4)') j+1,' 0 hex',
     :          nodes+ks,nodes+ks+1,nodes+ks+2,nodes+ks+3,
     :          nodes+ks+4,nodes+ks+5,nodes+ks+6,nodes+ks+7
      write(1,'(i4,a6,8i4)') j+2,' 0 hex',
     :          nodes+ks+8,nodes+ks+9,nodes+ks+10,nodes+ks+11,
     :          nodes+ks+12,nodes+ks+13,nodes+ks+14,nodes+ks+15

* 節点データ
      write(1,*) '1 0'
      write(1,*) '1 1'
      write(1,*) 'moment, kN*m'
      do i=1,k
        write(1,*) i, 0.
      end do
      do i=k+1,k+8
        write(1,*) i,cg
      end do
      do i=ks,nodes
        write(1,*) i,cs
      end do
      do i=nodes+1,nodes+k
        write(1,*) i, 0.
      end do
      do i=nodes+k+1,nodes+k+8
        write(1,*) i,cg
      end do
      do i=nodes+ks,2*nodes
        write(1,*) i,cs
      end do

* ここから各ステップの応答表示

      do istep=2,nstep
        time=time+dt                           ! 時間を進める
        read(2,*) pd,bd
        read(3,*) pd2,bd2
        bd=bd+pd                               ! 支承変位から桁変位に変換
        bd2=bd2+pd2

        write(moji,*) istep
        moji='step'//moji(2:)                  ! ステップ番号の数字を文字に変換
        write(1,'(a11,f6.2,a1)') moji,time,'s'
        write(1,*) 2*nodes, 2*members

* 橋脚応答
        zz=p
        do i=1,m+1,4
          zp=zz/p
          theta=3./2./p*zp*(2.-zp)*pd
          disp=zp*zp*(3.-zp)*pd/2.
          y(i    )= -cos(theta)*d+disp ; z(i    )=zz+sin(theta)*d
          y(i+  1)= -cos(theta)*d+disp ; z(i+  1)=zz+sin(theta)*d
          y(i+  2)=  cos(theta)*d+disp ; z(i+  2)=zz-sin(theta)*d
          y(i+  3)=  cos(theta)*d+disp ; z(i+  3)=zz-sin(theta)*d
          y(i+m+4)=s-cos(theta)*d+disp ; z(i+m+4)=zz+sin(theta)*d
          y(i+m+5)=s-cos(theta)*d+disp ; z(i+m+5)=zz+sin(theta)*d
          y(i+m+6)=s+cos(theta)*d+disp ; z(i+m+6)=zz-sin(theta)*d
          y(i+m+7)=s+cos(theta)*d+disp ; z(i+m+7)=zz-sin(theta)*d
          zz=zz-dp
        end do

        zz=p
        do i=1,m+1,4
          zp=zz/p
          theta=3./2./p*zp*(2.-zp)*pd2
          disp=zp*zp*(3.-zp)*pd2/2.
          v(i    )= -cos(theta)*d+disp ; w(i    )=zz+sin(theta)*d
          v(i+  1)= -cos(theta)*d+disp ; w(i+  1)=zz+sin(theta)*d
          v(i+  2)=  cos(theta)*d+disp ; w(i+  2)=zz-sin(theta)*d
          v(i+  3)=  cos(theta)*d+disp ; w(i+  3)=zz-sin(theta)*d
          v(i+m+4)=s-cos(theta)*d+disp ; w(i+m+4)=zz+sin(theta)*d
          v(i+m+5)=s-cos(theta)*d+disp ; w(i+m+5)=zz+sin(theta)*d
          v(i+m+6)=s+cos(theta)*d+disp ; w(i+m+6)=zz-sin(theta)*d
          v(i+m+7)=s+cos(theta)*d+disp ; w(i+m+7)=zz-sin(theta)*d
          zz=zz-dp
        end do

* 桁応答
        y(k+1)= -d+bd
        y(k+2)= -d+bd
        y(k+3)=s+d+bd
        y(k+4)=s+d+bd
        y(k+5)= -d+bd
        y(k+6)= -d+bd
        y(k+7)=s+d+bd
        y(k+8)=s+d+bd

        v(k+1)= -d+bd2
        v(k+2)= -d+bd2
        v(k+3)=s+d+bd2
        v(k+4)=s+d+bd2
        v(k+5)= -d+bd2
        v(k+6)= -d+bd2
        v(k+7)=s+d+bd2
        v(k+8)=s+d+bd2

* 支承応答
        theta=3.*pd/2./p
        yy=0.
        do i=ks,ke,8
          y(i  )=-d/2.+bd+yy
          y(i+1)=-d/2.+bd+yy
          y(i+2)= d/2.+bd+yy
          y(i+3)= d/2.+bd+yy
          y(i+4)= -cos(theta)*d/2.+pd+yy ; z(i+4)=p+sin(theta)*d/2.
          y(i+5)= -cos(theta)*d/2.+pd+yy ; z(i+5)=p+sin(theta)*d/2.
          y(i+6)=  cos(theta)*d/2.+pd+yy ; z(i+6)=p-sin(theta)*d/2.
          y(i+7)=  cos(theta)*d/2.+pd+yy ; z(i+7)=p-sin(theta)*d/2.
          yy=yy+s
        end do

        theta=3.*pd2/2./p
        yy=0.
        do i=ks,ke,8
          v(i  )=-d/2.+bd2+yy
          v(i+1)=-d/2.+bd2+yy
          v(i+2)= d/2.+bd2+yy
          v(i+3)= d/2.+bd2+yy
          v(i+4)= -cos(theta)*d/2.+pd2+yy ; w(i+4)=p+sin(theta)*d/2.
          v(i+5)= -cos(theta)*d/2.+pd2+yy ; w(i+5)=p+sin(theta)*d/2.
          v(i+6)=  cos(theta)*d/2.+pd2+yy ; w(i+6)=p-sin(theta)*d/2.
          v(i+7)=  cos(theta)*d/2.+pd2+yy ; w(i+7)=p-sin(theta)*d/2.
          yy=yy+s
        end do

* 出力
        do i=1,nodes
          write(1,*) i,x(i),y(i),z(i)
        end do
        do i=1,nodes
          write(1,*) nodes+i,u(i),v(i)+a,w(i)
        end do

        j=1
        do i=1,m,4
          write(1,'(i4,a6,8i4)') 
     :          j,' 0 hex',i,i+1,i+2,i+3,i+4,i+5,i+6,i+7
          j=j+1
        end do
        do i=m+5,k-7,4
          write(1,'(i4,a6,8i4)') 
     :          j,' 0 hex',i,i+1,i+2,i+3,i+4,i+5,i+6,i+7
          j=j+1
        end do
        write(1,'(i4,a6,8i4)') j,' 0 hex',
     :          k+1,k+2,k+3,k+4,k+5,k+6,k+7,k+8
        write(1,'(i4,a6,8i4)') j+1,' 0 hex',
     :          ks,ks+1,ks+2,ks+3,ks+4,ks+5,ks+6,ks+7
        write(1,'(i4,a6,8i4)') j+2,' 0 hex',
     :          ks+8,ks+9,ks+10,ks+11,ks+12,ks+13,ks+14,ks+15

        j=j+3
        do i=nodes+1,nodes+m,4
          write(1,'(i4,a6,8i4)') 
     :            j,' 0 hex',i,i+1,i+2,i+3,i+4,i+5,i+6,i+7
          j=j+1
        end do
        do i=nodes+m+5,nodes+k-7,4
          write(1,'(i4,a6,8i4)') 
     :            j,' 0 hex',i,i+1,i+2,i+3,i+4,i+5,i+6,i+7
          j=j+1
        end do
        write(1,'(i4,a6,8i4)') j,' 0 hex',
     :          nodes+k+1,nodes+k+2,nodes+k+3,nodes+k+4,
     :          nodes+k+5,nodes+k+6,nodes+k+7,nodes+k+8
        write(1,'(i4,a6,8i4)') j+1,' 0 hex',
     :          nodes+ks,nodes+ks+1,nodes+ks+2,nodes+ks+3,
     :          nodes+ks+4,nodes+ks+5,nodes+ks+6,nodes+ks+7
        write(1,'(i4,a6,8i4)') j+2,' 0 hex',
     :          nodes+ks+8,nodes+ks+9,nodes+ks+10,nodes+ks+11,
     :          nodes+ks+12,nodes+ks+13,nodes+ks+14,nodes+ks+15

* 力
        write(1,*) '1 0'
        write(1,*) '1 1'
        write(1,*) 'moment, kN*m'
* 橋脚
        zz=p
        do i=1,nn*4+1,4
          force=abs(pd*(p-zz)/p)
          write(1,*) i    , force
          write(1,*) i+  1, force
          write(1,*) i+  2, force
          write(1,*) i+  3, force
          zz=zz-dp
        end do
        zz=p
        do i=(nn+1)*4+1,k-3,4
          force=abs(pd*(p-zz)/p)
          write(1,*) i    , force
          write(1,*) i+  1, force
          write(1,*) i+  2, force
          write(1,*) i+  3, force
          zz=zz-dp
        end do
* 桁
        do i=k+1,k+8
          write(1,*) i,cg
        end do
* 支承
        do i=ks,nodes
          write(1,*) i,cs
        end do

        zz=p
        do i=nodes+1,nodes+nn*4+1,4
          force=abs(pd2*(p-zz)/p)
          write(1,*) i    , force
          write(1,*) i+  1, force
          write(1,*) i+  2, force
          write(1,*) i+  3, force
          zz=zz-dp
        end do
        zz=p
        do i=nodes+(nn+1)*4+1,nodes+k-3,4
          force=abs(pd2*(p-zz)/p)
          write(1,*) i    , force
          write(1,*) i+  1, force
          write(1,*) i+  2, force
          write(1,*) i+  3, force
          zz=zz-dp
        end do
        do i=nodes+k+1,nodes+k+8
          write(1,*) i,cg
        end do
        do i=nodes+ks,2*nodes
          write(1,*) i,cs
        end do

      end do
      close(1)
      close(2)

      end