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

! 実行例：
! $ ./meck<hz.dat>hzmeck.d
! gnuplot< pl 'hzmeck.d' using 1:2,'hzmeck.d' using 3:4,
!             'meck.d' using 1:2 w l,'meck.d' using 3:4 w l
!  hzmeck.dに座屈荷重を出力

implicit none
!---入力---
integer :: nstep=26            !     :出力ステップ数(CGXで確認)
integer :: pl1=26              !     :線形回帰に使う終わりのプロット
integer :: pl2=16              !     :線形回帰に使わない終りのプロット
integer :: nx=6                !     :幅方向の分割数
integer :: ny=30               !     :高さ方向の分割数
real(8) :: h=25d-2             ! m   :高さ
real(8) :: zakutu=0.1699666d-2 ! MN  :起こしたいモードの座屈荷重   
!---
real(8),allocatable :: u(:,:)  ! m   :横(x)方向の節点変位
real(8),allocatable :: v(:,:)  ! m   :荷重(y)方向の節点変位
real(8),allocatable :: d(:)    ! m   :横(x)方向の変位
real(8),allocatable :: phi(:)  ! rad :ねじれ角
real(8) :: p(100)              ! MN  :荷重
real(8)    xd,yd,xyd,xxd       !     :最小二乗法に使用
real(8)    xp,yp,xyp,xxp       !     :最小二乗法に使用
real(8)    ad,bd               !     :y=ax+bのa,b
real(8)    ap,bp               !     :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(v(nitimen,nstep))
allocate(d(nstep),phi(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),v(i,j),rcut
   enddo
enddo
!横方向変位の計算
do j=1,nstep
   d(j)=abs(sum(u(:,j))/size(u(:,j)))
enddo
!ねじれ角の計算
do j=1,nstep
   phi(j)=asin(abs(u(nx/2+1,j)-u(nitimen-nx/2,j))/h)
enddo
!gnuplot用のデータの出力
print*,'# delta (m)    phi/P (deg/MN)    phi (deg)    delta/P (m/MN)'
do j=1,nstep
   print'(1pe13.5,1pe13.5,1pe13.5,1pe13.5)',d(j),phi(j)/p(j),phi(j),d(j)/p(j)
enddo

!初期化
xd=0.d0
yd=0.d0
xyd=0.d0
xxd=0.d0
xp=0.d0
yp=0.d0
xyp=0.d0
xxp=0.d0
!最小二乗法で線形回帰
do j=nstep-pl1+1,nstep-pl2
   xd=xd+d(j)
   yd=yd+phi(j)/p(j)
   xyd=xyd+d(j)*phi(j)/p(j)
   xxd=xxd+d(j)**2.0
   xp=xp+phi(j)
   yp=yp+d(j)/p(j)
   xyp=xyp+phi(j)*d(j)/p(j)
   xxp=xxp+phi(j)**2.0
enddo
ad=((pl1-pl2)*xyd-xd*yd)/((pl1-pl2)*xxd-xd**2.)
bd=(xxd*yd-xyd*xd)/((pl1-pl2)*xxd-xd**2.)
ap=((pl1-pl2)*xyp-xp*yp)/((pl1-pl2)*xxp-xp**2.)
bp=(xxp*yp-xyp*xp)/((pl1-pl2)*xxp-xp**2.)
!座屈荷重の出力
print*,''
print*,'# 座屈荷重(kN):',1./sqrt(ad*ap)*10.**3.
!gnuplot用のデータの出力
open(11,file='meck.d')
write(11,*) 0.d0,bd,0.d0,bp
write(11,*) d(nstep-pl2),ad*d(nstep-pl2)+bd,phi(nstep-pl2),ap*phi(nstep-pl2)+bp
close(11)

end
