!CCXで計算した荷重−変位関係をgnuplot用にまとめる

! 実行例：
! $ ./pd<masatu.dat>pm.d
! gnuplot< pl 'pm.d' w lp

implicit none
!---入力---
integer :: nstep=82            !    :出力ステップ数(CGXで確認)
integer :: nx=30               !    :試験材の幅方向の分割数
integer :: ny=4                !    :試験材の高さ方向の分割数
integer :: nz=30               !    :試験材の引張方向の分割数
real    :: pz=2.d-6            ! MN :最大引張荷重
!---
real(8),allocatable :: w(:,:)  ! m  :引張方向(z)の節点変位
real(8),allocatable :: d(:)    ! m  :引張方向(z)の平均変位
real(8),allocatable :: p(:)    ! MN :荷重
integer    nitimen             !    :断面の節点数
integer    i,j,icut
real       rcut
character  ccut*20 
real(8)    step
nitimen=(nx+1)*(ny+1)
allocate(w(nitimen,nstep))
allocate(d(nstep),p(nstep))

!初期化
w=0.d0
d=0.d0
p=0.d0
do j=1,nstep
   !datファイルのデータの読み飛ばし
   read(*,'()')
   !datファイルの荷重の割合を読み込み
   read(*,*) ccut,ccut,ccut,ccut,ccut,ccut,ccut,ccut,ccut,step
   p(j)=(step-1.d0)*pz
   read(*,'()')
   !datファイルの節点変位の読み込み
   do i=1,nitimen
      read(*,*) icut,rcut,rcut,w(i,j)
   enddo
enddo
!横方向変位と荷重の計算
do j=1,nstep
   d(j)=sum(w(:,j))/size(w(:,j))
enddo
!gnuplot用のデータの出力
print*,'# delta (mm)    P (N)'
do j=1,nstep
   print'(1pe13.5,1pe13.5)',d(j)*1.d3,p(j)*1.d6
enddo

end
