!Southwell法で座屈荷重を計算
! 片持ち梁の横ねじれ座屈の非線形解析(CCX)シリーズに対応

! 実行例：
! $ ./sw<hyoko.dat>hzsw.d
! gnuplot< pl 'hzsw.d','sw.d' w l
!  hzsw.dに座屈荷重を出力

implicit none
!---入力---
integer :: nstep=96            !    :出力ステップ数(CGXで確認)
integer :: pl1=96              !    :線形回帰に使う終わりのプロット
integer :: pl2=20              !    :線形回帰に使わない終りのプロット
integer :: nx=8                !    :幅方向の分割数
integer :: ny=30               !    :高さ方向の分割数
real(8) :: zakutu=0.4682880d-1 ! MN :起こしたいモードの座屈荷重   
!---
real(8),allocatable :: u(:,:)  ! m  :横(x)方向の節点変位
real(8),allocatable :: d(:)    ! m  :横(x)方向の変位
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
character  ccut*20 
real(8)    step,pcr
nitimen=(nx+1)*(ny+1)
pcr=zakutu/0.8                 ! MN :ちょっと大きめに荷重をかけた
allocate(u(nitimen,nstep))
allocate(d(nstep),p(nstep))

!初期化
u=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*pcr
   read(*,'()')
   !datファイルの節点変位の読み込み
   do i=1,nitimen
      read(*,*) icut,u(i,j),rcut,rcut
   enddo
enddo
!横方向変位と荷重の計算
do j=1,nstep
   d(j)=abs(sum(u(:,j))/size(u(:,j)))
enddo
!gnuplot用のデータの出力
print*,'# delta (m)  delta/P (m/MN)'
do j=1,nstep
   print'(1pe13.5,1pe13.5)',d(j),d(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+d(j)
   y=y+d(j)/p(j)
   xy=xy+d(j)**2/p(j)
   xx=xx+d(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):',1.d0/a*1.d3
!gnuplot用のデータの出力
open(11,file='sw.d')
write(11,*) 0.d0,b
write(11,*) d(nstep-pl2),(a*d(nstep-pl2)+b)
close(11)

end
