!P-δ^2法で座屈荷重を計算
!オンサイト木橋の局部座屈用に書き換えた

! 実行例：
! 一番下の最小二乗法のプロット用データファイルの名前を変えてから……
! $ ./pd2<kyo500.dat>kyo500.d
! gnuplot< pl 'kyo500.d','axby_kyo500.d' w l

implicit none
!---入力---
integer :: nstep=100           !    :出力ステップ数(CGXで確認)
integer :: pl1=25              !    :線形回帰に使う終わりのプロット
integer :: pl2=0               !    :線形回帰に使わない終りのプロット
integer :: nx=4                !    :鋼板の幅方向の分割数
integer :: ny1=8               !    :木材の高さ方向の分割数
integer :: ny2=16              !    :鋼板の上下木材間の高さ方向の分割数
real(8) :: pmax=500.d-3        ! MN :載荷した最大の荷重   
!---
real(8),allocatable :: u(:,:)  ! m  :横(x)方向の節点変位
real(8),allocatable :: u2(:,:) ! m  :横(x)方向の節点変位の2乗
real(8),allocatable :: d2(:)   ! m  :u2の平均値
real(8),allocatable :: p(:)    ! MN :荷重
real(8)    x,y,xy,xx           !    :最小二乗法に使用
real(8)    a,b                 !    :y=ax+bのa,b
integer    nitimen             !    :断面の節点数
integer    i,j,icut
real       rcut,ny
character  ccut*20 
real(8)    step,pcr
ny=ny1*2+ny2
nitimen=(nx+1)*(ny+1)
allocate(u(nitimen,nstep),u2(nitimen,nstep))
allocate(d2(nstep),p(nstep))

!---初期化
u=0.d0
u2=0.d0
d2=0.d0
p=0.d0
do j=1,nstep
   read(*,'()')!datファイルのデータの読み飛ばし
   read(*,*) ccut,ccut,ccut,ccut,ccut,ccut,ccut,ccut,ccut,step!datファイルの荷重の割合を読み込み
   p(j)=(step-1.d0)*pmax!プレストレスのステップtime=1.0分だけ引く
   read(*,'()')!datファイルの節点変位の読み込み
   do i=1,nitimen
      read(*,*) icut,u(i,j),rcut,rcut
      u2(i,j)=(u(i,j)-u(i,1))**2!プレストレス力による変位を基準に初期化し、節点ごとにそれを2乗
   enddo
enddo
do j=1,nstep
   d2(j)=sum(u2(:,j))/size(u2(:,j))!節点変位の2乗をスッテプごとに平均化
enddo
!---gnuplot用のデータの出力
print*,'# delta^2 (m)    P (MN)'
do j=1,nstep
   print'(1pe13.5,1pe13.5)',d2(j),p(j)
enddo
!---初期化
x=0.d0
y=0.d0
xy=0.d0
xx=0.d0
!---最小二乗法で線形回帰
do j=nstep-pl1+1,nstep-pl2
   x=x+d2(j)
   y=y+p(j)
   xy=xy+d2(j)*p(j)
   xx=xx+d2(j)**2
enddo
a=((pl1-pl2)*xy-x*y)/((pl1-pl2)*xx-x**2)
b=(xx*y-xy*x)/((pl1-pl2)*xx-x**2)
!--座屈荷重の出力
print*,''
print*,'# 座屈荷重(kN):',b*10.d0**3
!#--gnuplot用のデータの出力
open(11,file='axby_kyoa500.d')
write(11,*) 0.d0,b
write(11,*) d2(nstep-pl2),a*d2(nstep-pl2)+b
close(11)

end
