c （薄肉とは見なせない）長方形断面のねじり定数、そりねじり定数を
c 断面積分で求めるプログラム。
c ねじり応力関数は、
c ティモシェンコ・グーディア（金多 潔 監訳）『弾性論』（コロナ社）
c 4版p.325の応力関数を用いた。
c 正規化そり関数は、
c 下関正義・藤沼平一・茅野寛奈：反りを考慮した長方形断面
c アイソパラメトリック要素, ばね論文集, 第39号, 1994, p.33-39
c の p.34 の正規化反り関数を用いた。この正規そり関数の２乗を
c 断面積分して そりねじり定数を求めたけど、比較できる
c 正解を知らないので、当たってるかどうか分からない。
c とはいえ、長辺が短辺の10倍以上の断面だと薄肉近似（後述）
c の解に近づくから、まあ、そうおかしくもないのかも。
c （比較解を知ってる人がいたら教えて）
c ねじり定数もそりねじり定数も、分割数を十分に細かくすれば、
c 短辺と長辺を逆に与えても同じ値になる。
c 有効数字3桁以上の精度が欲しいときは、おおよそ
c 級数に関しては7項以上、
c 分割数に関しては 100*100 以上が目安でしょうか
c 
c
c 03/11/6版
c プログラムを作った人：後藤文彦
c 
c このプログラムは無保証です。
c このプログラムは改造／再配布して構いません。
c 但し、私的利用に留まらない場合は、もとのプログラムが
c http://www.str.ce.akita-u.ac.jp/~gotou/programoj/hari/kukei.f
c から取ってきたものと分かるように明示して下さい。
c 改造／再配布後も上記の条件を保持して下さい。
c バグの報告は、
c http://www.str.ce.akita-u.ac.jp/~gotou/jpg/meeru.jpg
c にお願いします。
c
           program kukei
          implicit real*8 (a-h, o-z)
          dimension x(5000),y(5000)
c ねじりの応力関数、そり関数の級数の項数：kyuu
c          kyuu=19
5         continue
        write(*,*) '-------------------------------------------------'
        write(*,*) 'ねじり応力関数の級数の項数（奇数）を入れて下さい'
        read(*,*) kyuu
        write(*,*) '長方形の短辺と長辺を入力して下さい'
        read(*,*) tanp,tyou
          a=tanp/2.d0
          b=tyou/2.d0
        write(*,*) '短辺の分割数と長辺の分割数を偶数で入力して下さい'
        read(*,*) na,nb
          na=na/2
          nb=nb/2
          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)
           yi=yi+xl**2*f
           xi=xi+ym**2*f
          call neziri(kyuu,a,b,xl,ym,pi,p)
           pj=pj+p*f
           call sori(kyuu,a,b,xl,ym,pi,w)
           ww=ww+w**2*f
c
c (0<x<a,0<y<b),(-a<x<0,0<y<b),(0<x<a,-b<y<0),(-a<x<0,-b<y<0)
c を別々に計算する場合（弾性なら符号が変わるだけだからw**2は
c どの象限でも同じになる
c          call sori(kyuu,a,b,xl,ym,pi,w1)
c          call sori(kyuu,a,b,-xl,ym,pi,w2)
c          call sori(kyuu,a,b,xl,-ym,pi,w3)
c          call sori(kyuu,a,b,-xl,-ym,pi,w4)
c          ww=ww+(w1**2+w2**2+w3**2+w4**2)*f
c
100      continue
c
c 各断面定数は、0<x<a,0<y<b の領域だけで積分したものを
c 4倍してるけど、ファイバーモデルなどで各断面分割要素の弾塑性を
c 考慮して弾性係数をかけながら積分するときは、ちゃんと
c -a<x<a, -b<y<b の領域で断面積分すること。
         xi=xi*4.d0
         yi=yi*4.d0
         pj=64.d0*a**2/pi**3*pj*4.d0
         ww=ww*4.d0
c
c
       write(*,500) ' 断面積：',tanp*tyou
       write(*,500) ' 強軸回りの断面二次モーメント（断面積分）：',xi
       write(*,500) ' 強軸回りの断面二次モーメント（ 正  解 ）：',
     &               tanp*tyou**3/12.d0
       write(*,500) ' 弱軸回りの断面二次モーメント（断面積分）：',yi
       write(*,500) ' 弱軸回りの断面二次モーメント（ 正  解 ）：',
     &               tyou*tanp**3/12.d0
       write(*,500) ' ねじり定数（断面積分）：                ',pj
          call timj(a,b,pi,tj)
               bt33=tyou*tanp**3/3.d0
       write(*,500) ' ねじり定数（ 正  解 ）：                ',tj*bt33
          call timj2(kyuu,a,b,pi,tj2)
       write(*,500) ' ねじり定数（積分は正解、級数は有限個）：',tj2
       write(*,500) ' ねじり定数（薄肉の近似：bt^3/3）：      ',bt33
       write(*,500) ' そりねじり定数（断面積分）：            ',ww
       write(*,500) ' そりねじり定数（薄肉近似：(bt)^3/144）：',
     &                (tyou*tanp)**3/144.d0
c (bt)^3/144 というのは _|_ 型断面の出っ張り | の厚さを
c 0 としたそりねじり定数
          goto 5
500      format(a,1pd15.7)
         end
c
c ******* ねじり定数 ********
c ねじり応力関数の級数について和を取る
c x,y について断面積分はメインでやる
         subroutine neziri(kyuu,a,b,x,y,pi,p)
          implicit real*8 (a-h, o-z)
           p=0.d0
          do 10 n=1,kyuu,2
           p=p+(-1.d0)**( (n-1)/2 )/n**3 *
     &     (        1.d0-
     &       cosh(n*pi*y/2.d0/a) / cosh(n*pi*b/2.d0/a)
     &     )        *cos(n*pi*x/2.d0/a) 
10        continue
          return
          end
c
c Timoshenko が計算した正解
c 手で積分して、級数も極限を取れる部分は取ってある
c 極限を取れない級数は、31項まで
         subroutine timj(a,b,pi,tj)
          implicit real*8 (a-h, o-z)
           tj=0.d0
           kyuu=31
          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 積分は（Timoshenkoが）手で計算した正解だけど、
c 級数は、指定した項数までしか取らない
         subroutine timj2(kyuu,a,b,pi,tj2)
          implicit real*8 (a-h, o-z)
           tn=0.d0
           tt=0.d0
           tj2=0.d0
          do 10 n=1,kyuu,2
           tn=tn+1.d0/n**4
           tt=tt+1.d0/n**5*tanh(n*pi*b/2.d0/a)
10        continue
           tj2=512.d0*a**3*b/pi**4*tn-1024.d0*a**4/pi**5*tt
          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