! yokonezire.f90(yone.f90)を田村さんが座屈設計ガイドラインにのっている
! 各種の境界条件、荷重条件の式で解けるように改造したもの。下記は
! 両端固定、中央集中荷重の例。ただし、長辺/短辺>50ぐらいになると、
! sinh()がInfinityとなり計算できなくなるので、その場合は、
! J=ht**3/3, Iw=(ht)**3/144 の近似に置き換える。
        implicit real*8(a-h,o-z)
         pi=2.d0*asin(1.d0)
         b=19.5d0
         h=1000.d0
         ell=6000.d0
         e=206.d0*1.d3  !N/mm**2=MPa
         poi=0.3d0
         g=e/2.d0/(1.d0+poi)
		 cb1=0.938d0
		 cb2=0.715d0
		 cb3=4.8d0
		 hp=500.d0
          yi=h*b**3/12.d0
		  zi=b*h**3/12.d0
          eiy=e*yi
		  if(h/b>50) then
		   print*,"長辺/短辺>50となったので、JとIwは近似式を使います"
		   cj=h*b**3/3.d0
		   wi=(h*b)**3/144.d0
		  else
           call kukei(b,h,cj,wi)
		  end if
		  print*,"cj,wi=", cj, wi
		  print*,"h*b**3/3=",h*b**3/3.d0
		  print*,"(h*b)**2/144=",(h*b)**3/144.d0
          gj=g*cj
          eiw=e*wi
		  pk=0.5d0
		  !
		  cm1=(cb1*pi**2*e*zi)/(pk**2*ell**2)
		  cm2=cb2*hp
		  cm3=1.d0+((pk*ell)**2*gj)/(pi**2*eiw)
		  zakumo1=cm1*(cm2+sqrt(cm2**2+(wi/zi)*cm3))
		  zakumo2=cm1*(cm2-sqrt(cm2**2+(wi/zi)*cm3))
		  kaju1=(8*zakumo1)/ell
		  kaju2=(8*zakumo2)/ell
		  print*,'座屈モーメント+=',zakumo1
		  print*,'座屈モーメント-=',zakumo2
		  print*,'座屈荷重+=',kaju1
		  print*,'座屈荷重-=',kaju2
!
           a1=sqrt(eiy*gj)/ell/ell
           a2=sqrt(eiw/gj)*pi/ell
           tra=a1*(3.95d0+3.52d0*a2)
           print*,'tra=',tra,'N'
           print*,'tra=',tra/1.d3,'kN'
           print*,'tra=',tra/1.d6,'MN'
           print*,'tra=',tra/1.d9,'GN'
           ! tram=tra*ell**2/eiy
           !r=sqrt(yi/b/h)
           !rmd=ell/r
!          print '(a,1pd11.3)','wi=',wi
!          print '(f5.0, 3f8.2,1p4d11.3 )', ell,tra,pe,gosa,tram,pem,
!          print*, ell,tra,pe,gosa,tram,pem,gj,wi
!          write(*,100) ell,tra,pe,gosa,tram,pem,rmd
         end
! 
          subroutine kukei(tanp,tyou,cj,wi)
          implicit real*8 (a-h, o-z)
          dimension x(5000),y(5000)
! ねじりの応力関数、そり関数の級数の項数：kyuu
          kyuu=9
          a=tanp/2.d0
!		  print*,"a=",a
!		  read*, nanka
          b=tyou/2.d0
! 短辺/2の分割数:na
! と長辺/2の分割数:nb（偶数）
!          na=500
          na=4
          nb=500
          an=a/float(na)
          bn=b/float(nb)
          f=an*bn
          pi=2.d0*asin(1.d0)
        do i=1,na
         x(i)=an*i -an/2.d0
		end do
        do j=1,nb
         y(j)=bn*j -bn/2.d0
		end do
!
         xi=0.d0
         yi=0.d0
         pj=0.d0
         ww=0.d0
         do l=1,na
         do m=1,nb
           xl=x(l)
           ym=y(m)
!		   print*,"xl,ym",xl,ym
           call sori(kyuu,a,b,xl,ym,pi,w)
           ww=ww+w**2*f
		  end do
		  end do
!		  print*,"w,f=",w,f
           wi=ww*4.d0
!
          call timj(kyuu,a,b,pi,tj)
               bt33=tyou*tanp**3/3.d0
               cj=tj*bt33
!       write(*,500) ' ねじり定数（Timoshenkoの長方形断面）：  ',cj
!       write(*,500) ' ねじり定数（薄肉の近似：bt^3/3）：      ',bt33
!       write(*,500) ' そりねじり定数（断面積分）：            ',wi
!       write(*,500) ' そりねじり定数（薄肉近似：(bt)^3/144）：',
!     &                (tyou*tanp)**3/144.d0
! (bt)^3/144 というのは _|_ 型断面の出っ張り | の厚さを
! 0 としたそりねじり定数
! 500      format(a,1pd15.7)
         return
         end
!
! ******* ねじり定数 ********
! Timoshenko が計算した正解
! 手で積分して、級数も極限を取れる部分は取ってある
! 極限を取れない級数は、19項まで
         subroutine timj(kyuu,a,b,pi,tj)
          implicit real*8 (a-h, o-z)
           tj=0.d0
          do n=1,kyuu,2
           tj=tj+1.d0/n**5*tanh(n*pi*b/2.d0/a)
		  end do
           tj=1.d0-192.d0/pi**5*a/b*tj
          return
          end
!
!
! ******** そりねじり定数 **********
! 正規化そり関数の級数について和を取る
         subroutine sori(kyuu,a,b,x,y,pi,w)
          implicit real*8 (a-h, o-z)
          w=0.d0
          do n=1,kyuu,2
           w=w+(-1.d0)**( (n+1)/2 )/ n**3 *&
&       sinh(n*pi*y/2.d0/a) / cosh(n*pi*b/2.d0/a)&
&       *sin(n*pi*x/2.d0/a)
! sinhやcoshの()内が800程度を超えるとこれらの関数の値はInfinityとなる
! ので、目安として、()内が500以下程度に制限する必要がある。
!print*,"n*pi*y/2/a,n*pi*b/2/a=",n*pi*y/2.d0/a,n*pi*b/2.d0/a
!print*,"w,sinh,cosh,sin=",w,&
!&       sinh(n*pi*y/2.d0/a), cosh(n*pi*b/2.d0/a)&
!&       *sin(n*pi*x/2.d0/a)
!print*,"n,pi,y,a=",n,pi,y,a
           end do
           w=32.d0*a**2/pi**3* w + x*y
          return
          end
!

