! 引張試験のひずみを間引いて線形回帰＋ヤング率＋ポアソン比
! syokignu.fで初期不平衡値を直したデータを使用

implicit none
!---手入力---
integer :: ts=6,te=25           ! sec     : 載荷後ts~te秒のひずみデータを使用
real(8) :: bh=19.98*4.06/(10.**6.)   ! m^2     : 断面積
real(8) :: m=1.d0               ! kg      : おもり1個の重さ
!---
real(8) :: nama(1:540,4)         ! *10^6   : ひずみ(初期不平衡を直した)
real(8) :: g=9.8d0               ! m/sec^2 : 重力加速度
real(8) :: hei(0:17,4)           !         : ひずみの平均値
real(8) :: ou(0:17)              ! Pa      : 引張方向の応力
real(8) :: a(4)                  ! Pa      : ヤング率(y=ax+bのa)
real(8) :: b(4)                  ! Pa      : 回帰線の切片(y=ax+bのb)
real(8) :: x(4),y(4),xy(4),xx(4) !         : 最小二乗法に使用
integer i,j,k,n

!データの読み込み
do i=1,540!計測時間
read(*,*) n,(nama(i,j),j=1,4)
enddo
!ひずみと応力の算出
do i=0,17!重りの数
   do j=1,4
      hei(i,j)=0.d0
      do k=ts-1,te-1
         hei(i,j)=hei(i,j)+nama(30*i+k,j)
      enddo
      !ひずみの平均化と*10^-6
      hei(i,j)=hei(i,j)/((te-ts+1)*10.**6)
   enddo
   !応力の計算
   ou(i)=m*g*i/bh
enddo
!gnuplot用データ作成
open(11,file='t1')
open(12,file='t2')
open(13,file='y1')
open(14,file='y2')
do i=0,17
   do j=1,4
      write(10+j,*) hei(i,j),ou(i)
   enddo
enddo
close(11)
close(12)
close(13)
close(14)
!最小二乗法で線形回帰
print*, '       E縦1(GPa)                 E縦2(GPa)'
xx=0.d0
x=0.d0
do j=1,4
   do i=0,17
      x(j)=x(j)+hei(i,j)
      xx(j)=xx(j)+(hei(i,j))**2
      xy(j)=xy(j)+hei(i,j)*ou(i)
      y(j)=y(j)+ou(i)
   enddo
enddo
a(:)=((17.+1.)*xy(:)-x(:)*y(:))/((17.+1.)*xx(:)-(x(:))**2.)
b(:)=(xx(:)*y(:)-xy(:)*x(:))/((17.+1.)*xx(:)-(x(:))**2.)
print*, (a(i)/(10.**9.),i=1,2)
print*,'    Eave(GPa) :', (a(1)+a(2))/(2.*10.**9.)
print*,'  ポアソン比1 :',(-(a(1))/(a(3)))
print*,'  ポアソン比2 :',(-(a(2))/(a(4)))
print*,'ポアソン比ave :',(-(a(1)+a(2))/(a(3)+a(4)))
!せん断弾性係数は、等方性として算出
print*,'       G1(GPa):',a(1)/(1.+(-(a(1))/(a(3))))/(2.*10.**9.)
print*,'       G2(GPa):',a(2)/(1.+(-(a(2))/(a(4))))/(2.*10.**9.)
print*,'     Gave(GPa):',((a(1)+a(2))/2)/(1.+(-(a(1)+a(2))/(a(3)+a(4))))/(2.*10.**9.)
print*,'------------'
print*,'1 y=ax+b の a1,b1(GPa) :',a(1)/(10.**9.),b(1)/(10.**9.)
print*,'2 y=ax+b の a2,b2(GPa) :',a(2)/(10.**9.),b(2)/(10.**9.)
!gnuplot用データ作成
open(11,file='t1k')
open(12,file='t2k')
do j=1,2
   write(10+j,*) -1./a(j)*b(j),0.d0 
   write(10+j,*) maxval(hei(:,j)),(a(j)*maxval(hei(:,j))+b(j))
enddo
close(11)
close(12)
end
