implicit none

real :: csarea, compst, acompst, bcompst, ccompst, ps(20000,2), x1, x2, y1, y2, a, b, f1, f2
real :: c, d, v1, v2, w1, w2, xx, s1p3, smodulus, l, s1, s2, s50, x
real :: strain(100000),stress(100000)

integer :: i,ios,data
character(1000) :: filename
character(1000) :: filename2
character(1000) :: filecsv
character(1000) :: filecsv2
character(4)    :: csv
character(1000) :: linebuf
character(1000) :: linebuf2

 print *,'このプログラムはJIS A 1149に準処し、静弾性係数を求めます。'
 print *,"・開くファイルは一列目に荷重(MN),二から三列目にひずみ(με)となっていること。"
 print *,"・荷重、ひずみはは圧縮によるものが負の値であること。"

 print *,"開くファイル名を入力（拡張子不要）"
  read '(A)',filename

 print *,"供試体の断面積(mm^2)を入力"
  read *,csarea

 print *,"最大圧縮強度(N/mm^2)を入力"
 print *,"(データロガーとのずれが生じるため)"
  read *,compst


 csv='.csv'
 filecsv=trim(adjustl(filename))//csv
 
 acompst = (compst*10)/3
 bcompst = nint(acompst)
 ccompst = bcompst/10


  open(12,file = filecsv)
   do i=1,40000
    read(12,*,iostat=ios),l,s1,s2
    
     if(ios < 0)exit
      
     ps(i,1) = (-l*(10**6))/csarea
     ps(i,2) = -(s1+s2)/2
     data=i
   end do
  close(12)
  
  print *,"読み込んだデータ数：",i

 x1=0
 x2=0
 y1=0
 y2=0
 do i=1,20000

  if( (x1 <= ps(i,2)).and.(ps(i,2) <= 50 ))then
   x1=ps(i,1)
    x2=ps(i+1,1)
   y1=ps(i,2)
    y2=ps(i+1,2)

  else
   if(y1 > 50.0001) exit                                   !0.0001は安全マージン

  end if
 end do


    a=((x2-x1)/(y2-y1))
    b=(x1-a*y1)

   x=50
s50=f1(x,a,b)

!print *,s50

 v1=0
 do i=1,20000

  if( (v1 <= ps(i,1)).and.(ps(i,1) <= ccompst))then
   v1=ps(i,1)
    v2=ps(i+1,1)
   w1=ps(i,2)
    w2=ps(i+1,2)
  else
   if(v1 > ccompst+0.001)exit

  end if
 end do

     c=((v2-v1)/(w2-w1))
     d=(-(v1/c)+w1)

   xx=ccompst  
s1p3=f2(xx,c,d)


 smodulus=( (ccompst-s50)/(s1p3-50) )*1000


 print *,"静弾性係数(kN/mm^2)は",smodulus


end



 function f1(x,a,b)
  real :: a,b,x

   f1=a*x+b

  return

 end



 function f2(xx,c,d)
  real ::c,d,xx

   f2=(xx/c)+d

  return

 end


