<br />
<b>Deprecated</b>:  The each() function is deprecated. This message will be suppressed on further calls in <b>/home/zhenxiangba/zhenxiangba.com/public_html/phproxy-improved-master/index.php</b> on line <b>456</b><br />
        implicit real*8(a-h,o-z)
         pi=2.d0*asin(1.d0)
         n=18
         do i=1,n
          read*, b,h,ell,e,g,pe
          yi=h*b**3/12.d0
          eiy=e*yi
          call kukei(b,h,cj,wi)
          gj=g*cj
          eiw=e*wi
           a1=sqrt(eiy*gj)/ell/ell
           a2=sqrt(eiw/gj)*pi/ell
           tra=a1*(3.95d0+3.52d0*a2)
           tram=tra*ell**2/eiy
           pem=pe*ell**2/eiy
           gosa=(pe-tra)/tra*100.d0
           r=sqrt(yi/b/h)
           rmd=ell/r
c          print '(a,1pd11.3)','wi=',wi
c          print '(f5.0, 3f8.2,1p4d11.3 )', ell,tra,pe,gosa,tram,pem,
c          print*, ell,tra,pe,gosa,tram,pem,gj,wi
          write(*,100) ell,tra,pe,gosa,tram,pem,rmd
100       format(f5.0, 3f8.2,1p2d11.3,f9.2)
         end do
         end
c 
          subroutine kukei(tan,tyou,cj,wi)
          implicit real*8 (a-h, o-z)
          dimension x(5000),y(5000)
c ねじりの応力関数、そり関数の級数の項数：kyuu
          kyuu=9
          a=tan/2.d0
          b=tyou/2.d0
c 短辺/2の分割数:na
c と長辺/2の分割数:nb（偶数）
          na=500
          nb=500
          an=a/float(na)
          bn=b/float(nb)
          f=an*bn
          pi=2.d0*asin(1.d0)
        do 10 i=1,na
         x(i)=an*i -an/2.d0
10      continue
        do 20 j=1,nb
         y(j)=bn*j -bn/2.d0
20      continue
c
         xi=0.d0
         yi=0.d0
         pj=0.d0
         ww=0.d0
         do 100 l=1,na
         do 100 m=1,nb
           xl=x(l)
           ym=y(m)
           call sori(kyuu,a,b,xl,ym,pi,w)
           ww=ww+w**2*f
100      continue
           wi=ww*4.d0
c
          call timj(kyuu,a,b,pi,tj)
               bt33=tyou*tan**3/3.d0
               cj=tj*bt33
c       write(*,500) ' ねじり定数（Timoshenkoの長方形断面）：  ',cj
c       write(*,500) ' ねじり定数（薄肉の近似：bt^3/3）：      ',bt33
c       write(*,500) ' そりねじり定数（断面積分）：            ',wi
c       write(*,500) ' そりねじり定数（薄肉近似：(bt)^3/144）：',
c     &                (tyou*tan)**3/144.d0
c (bt)^3/144 というのは _|_ 型断面の出っ張り | の厚さを
c 0 としたそりねじり定数
500      format(a,1pd15.7)
         return
         end
c
c ******* ねじり定数 ********
c Timoshenko が計算した正解
c 手で積分して、級数も極限を取れる部分は取ってある
c 極限を取れない級数は、19項まで
         subroutine timj(kyuu,a,b,pi,tj)
          implicit real*8 (a-h, o-z)
           tj=0.d0
          do 10 n=1,kyuu,2
           tj=tj+1.d0/n**5*tanh(n*pi*b/2.d0/a)
10        continue
           tj=1.d0-192.d0/pi**5*a/b*tj
          return
          end
c
c
c ******** そりねじり定数 **********
c 正規化そり関数の級数について和を取る
         subroutine sori(kyuu,a,b,x,y,pi,w)
          implicit real*8 (a-h, o-z)
          w=0.d0
          do 10 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)
10        continue
           w=32.d0*a**2/pi**3* w + x*y
          return
          end
c