!オンサイト木橋のリラクセーションの計算結果(CCX)をgnuplot用にまとめる
! (幅方向-緊張力-たわみ 関係)

! 実行例：
! $ gf plr plonsirelax.f90
! $ ./plr<relax.dat
! gnuplot< set zrange[] reverse
! gnuplot< spl 'db.d' w pm3d 

implicit none
!#--入力---
integer :: nstep=18            !    :datファイルの出力ステップ数
integer :: nx1=8               !    :一本の木材の幅員方向(x)分割数
integer :: nx2=4               !    :鋼板の厚さ方向(x)分割数
integer :: ny1=8               !    :上下木材の高さ方向(y)分割数
integer :: ny2=16              !    :上下木材間隔の高さ方向(y)分割数
real    :: h1=120d-3           !    :一本の木材の高さ=幅(y or z)
real    :: b1=840d-3/2         ! m  :木材全体の幅(x)/2
real    :: b2=9.0d-3           ! m  :鋼板の厚さ(x)
real    :: py=70.d-3           ! MN :荷重
real    :: px0=70.d-3          ! MN :導入時のプレストレス力
real    :: px=40.d-3           ! MN :低下後の緊張力
!---
real(8),allocatable :: v0(:,:) ! m  :変位計のの載荷方向(y)の節点変位
real(8),allocatable :: v1(:,:) ! m  :載荷した上縁の載荷方向(y)の節点変位
real(8),allocatable :: d1(:)   ! m  :載荷した上縁の載荷方向(y)の平均変位
real(8),allocatable :: p(:)    ! MN :荷重
integer    nitimen             !    :断面の節点数
integer    njouen              !    :上縁の節点数
integer    nkaen               !    :下縁の節点数
integer    i,j,k,l,icut
real       rcut,dx1,dx2
character  ccut*20 
real(8)    step
nitimen=(nx2+1)*(ny1*2+ny2+1)+(nx1+1)*(ny1+1)*int(b1/h1)*2&
    &+(nx1/2+1)*(ny1+1)*2*int(mod(b1,h1)/h1+501d-3)
njouen=(nx2+1)+(nx1+1)*int(b1/h1)+(nx1/2+1)*int(mod(b1,h1)/h1+501d-3)
nkaen=(nx1+1)*int(b1/h1)+(nx1/2+1)*int(mod(b1,h1)/h1+501d-3)
dx1=h1/nx1
dx2=b2/nx2
allocate(v0(njouen,nstep))
allocate(v1(njouen,nstep))
allocate(d1(nstep),p(nstep))

!#--出力ファイル名
open(12,file='relax04.d',status='replace')

v0=0.d0
v1=0.d0
d1=0.d0
p=0.d0
do j=1,nstep
   !変位計(henni)のデータ
   read(*,'()')
   read(*,*) ccut,ccut,ccut,ccut,ccut,ccut,ccut,ccut,ccut,step
   p(j)=(step-1.d0)*py*1.d3
   read(*,'()')
   do i=1,njouen
      read(*,*) icut,rcut,v0(i,j),rcut
   enddo
   !上縁(saika)のデータ
   read(*,'()')
   read(*,*) ccut,ccut,ccut,ccut,ccut,ccut,ccut,ccut,ccut,step
   p(j)=(px0+(px-px0)*(step-2.d0))*1.d3
   read(*,'()')
   do i=1,njouen
      read(*,*) icut,rcut,v1(i,j),rcut
   enddo
enddo
!---断面の変位を平均化
do j=1,nstep
   d1(j)=sum(v1(:,j))/size(v1(:,j))*1.d3
enddo
!---gnuplot用のデータの出力(幅-緊張力-たわみ)
write(12,'(A)') '# haba (mm)   Px(kN)   tawami(mm)'
do j=1,nstep
   i=0
   do k=1,nx2+1!鋼板
      i=i+1
      write(12,'(1pe13.5,1pe13.5,1pe13.5)') (dx2*(k-1))*1.d3,p(j),v0(i,j)*1.d3
   enddo
   do l=1,int(b1/h1)!木材
      do k=1,nx1+1
         i=i+1
         write(12,'(1pe13.5,1pe13.5,1pe13.5)') (b2+h1*(l-1)+dx1*(k-1))*1.d3,p(j),v0(i,j)*1.d3
      enddo
   enddo
   if (mod(b1,h1)>h1/2-1.d-3.and.mod(b1,h1)<h1/2+1.d-3) then!木材の本数が奇数の時
      do k=1,nx1/2+1
         i=i+1
         write(12,'(1pe13.5,1pe13.5,1pe13.5)') (b2+h1*int(b1/h1)+dx1*(k-1))*1.d3,p(j),v0(i,j)*1.d3
      enddo
   endif
   if (mod(b1,h1)>h1/2-1.d-3.and.mod(b1,h1)<h1/2+1.d-3) then!折り返し(木材の本数が奇数の時)
      do k=nx1/2-1,0,-1
         i=i-1
         write(12,'(1pe13.5,1pe13.5,1pe13.5)') (b2+h1*int(b1/h1)+h1-dx1*k)*1.d3,p(j),v0(i,j)*1.d3
      enddo
   endif
   do l=int(b1/h1),1,-1!折り返し木材
      do k=nx1,0,-1
         i=i-1
         write(12,'(1pe13.5,1pe13.5,1pe13.5)') (b2+b1*2-h1*(l-1)-dx1*k)*1.d3,p(j),v0(i,j)*1.d3
      enddo
   enddo
   do k=nx2,0,-1!鋼板
      i=i-1
      write(12,'(1pe13.5,1pe13.5,1pe13.5)') (b2*2+b1*2-dx2*k)*1.d3,p(j),v0(i,j)*1.d3
   enddo
      write(12,*) ('')
enddo

end
