        implicit real*8 (a-h,o-z)
        pi=2.d0*asin(1.d0)
        ex=210.d3!ヤング率
        ey=ex
        poi=0.3d0
        g=ex/2.d0/(1.d0+poi)!せん断弾性係数
!        g=80.8d3!せん断弾性係数
!        g=80.8d3!せん断弾性係数
!        ell=120. !梁の軸長
        ell=50.d0
        ell=40.d0
        do
        print*,'梁の長さを入れてください'
        read*,ell
        if(ell<0) exit
!        b=10. !フランジ幅
        b=300.d-2
        b=300.d-3
        !b=1.
!        tf=2. !フランジ厚
        tf=20.d-2
        tf=20.d-3
        !tf=0.05
!        h=12. !フランジ図心間距離
        h=568.d-2
        h=568.d-3
        !h=1.95
!        tw=2. !ウェブ厚
        tw=12.d-2
        tw=12.d-3
        !tw=1.
        a=2.d0*b*tf+(h-tf)*tw!断面積
        xi=(b*(h+tf)**3-(b-tw)*(h-tf)**3)/12.d0!x軸周りの断面2次モーメント
        yi=(2.d0*tf*b**3+(h-tf)*tw**3)/12.d0!y軸回りの断面2次モーメント
        ej=(2.d0*b*tf**3+(h-tf)*tw**3)/3.d0!ねじり定数
        ok=pi/ell*sqrt(ex*yi*h**2/4/g/ej)
!        oc1=2.462/sqrt(1+ok**2)+2.383*ok/sqrt(1+ok**2)
        oc1=2.437/sqrt(1+ok**2)+0.613*ok/sqrt(1+ok**2)
     &        -0.105*ok**2/sqrt(1+ok**2)
        oc2=1.
        oc3=1.
        kz=2.
        kw=1.
        zg=0.
        zj=0.
        wi=b**3*h**2*tf/24.d0
        a1=oc2*zg-oc3*zj
        a2=sqrt((kz/kw)**2*wi/yi+(kz*ell)**2*g*ej/pi**2/ex/yi+a1**2)
        pcr=oc1*pi**2*ex*yi/(kz*ell)**2*(a2-a1)
        p=pcr/ell 
        print*,'座屈荷重は',p,'MN'
!        print*,'座屈荷重は',p/1.d3,'KN'
!        print*,'座屈荷重は',p/1.d6,'MN'
!        
!
        print*,'A=',a
        print*,'Ix=',xi
        print*,'Iy=',yi
        print*,'J=',ej
        print*,'Iw=',wi
         eix=ex*xi
         eiy=ey*yi
         ea=a*sqrt(ex*ey)
         eiw=wi*sqrt(ex*ey)
         gj=g*ej
         el=ell/float(nofe)
c 片持ち梁の横ねじれ座屈のトラヘアの解（座屈前面内変位無視）
         a1=sqrt(eiy*gj)/ell/ell
         a2=sqrt(eiw/gj)*pi/ell
         tra=a1*(3.95d0+3.52d0*a2)
         write(*,*) 'Trahair=',tra,'MN'
         write(*,*) 'Trahair=',tra*1.d3,'kN'
!         write(*,*) 'Trahair=',tra/1.d6,'MN'
         print*,'L^3/EIy=',ell**3/eiy
!         
         end do
          end 
