!CLTの弾塑性解析の計算結果(CCX)をgnuplot用にまとめる
! clt1so.f90 と clt8so.f90 に対応
! (荷重-たわみ(載荷点の) 関係)

! 実行例：
!　まずは、datファイルの載荷節点数とかを読み取ってから……、
! $ gf plc plcltso.f90
! $ ./plc<8so.dat>8so.d
! gnuplot< pl '8so.d' w lp 

implicit none
!#--入力---
integer :: nstep=50            !    :datファイルの出力ステップ数
integer :: nsaika=45           !    :datファイルに書かれた載荷節点数
real    :: pmax=1000.d-3       ! MN :最大荷重
!---
real(8),allocatable :: v(:,:)  ! m  :載荷点の載荷方向(y)の節点変位
real(8),allocatable :: d(:)    ! m  :載荷点の載荷方向(y)の平均変位
real(8),allocatable :: p(:)    ! MN :その時の荷重
integer    i,j,k,l,icut
real       rcut
character  ccut*20 
real(8)    step
allocate(v(nsaika,nstep))
allocate(d(nstep),p(nstep))

v=0.d0
d=0.d0
p=0.d0
do j=1,nstep
   do i=1,7!変位計1/2/3/7/4/5/6のデータ
      do k=1,4!4行分
         read(*,'()')!読み飛ばし
      enddo
   enddo
   read(*,'()')
   read(*,*) ccut,ccut,ccut,ccut,ccut,ccut,ccut,ccut,ccut,step!スッテプを読み取る
   p(j)=(pmax*step)*1.d3!(kN)
   read(*,'()')
   do i=1,nsaika!載荷節点数分
      read(*,*) icut,rcut,v(i,j),rcut!節点変位を読み取る
   enddo
enddo
!---断面の変位を平均化
do j=1,nstep
   d(j)=sum(v(:,j))/size(v(:,j))*1.d3
enddo
!---gnuplot用のデータの出力(荷重-たわみ)
print'(A)','#   d(mm)   P(kN)'
do j=1,nstep
      print'(1pe13.5,1pe13.5)',d(j),p(j)
enddo

end
