c harisen.f で(単純支持された)柱のオイラー座屈を解かせるための
c 入力データを作るプログラム。oirasen を実行させると oirasen.d という
c データファイルができるから、harisen<oirasen.d と実行させればよい。
c
          program oirasen
        implicit real*8 (a-h,o-z)
        dimension alp(256),gmm(256),phi(256)
        write(*,*) '要素数を入れて下さい'
        read(*,*) nofe
        write(*,*)'座屈モードを入れて下さい'
        read(*,*) mode
        write(*,*)'外乱を入れる場合は、その倍率を入れて下さい'
        write(*,*)'(0 を入れた場合は座屈点を捜します)'
        read(*,*) gmode
        write(*,*)'モード計算のための出力が要りますか　要1／不要0 '
        read(*,*) modeo  
        pi=2.d0*asin(1.d0)
c  断面定数
c ヤング率：e
        e=200.d6
        e=1.d7 * 1.d3
c せん断弾性係数：g
        g=77.2d6
        g=7.d5 * 1.d3
c ポアソン比：poi
c 等方性材料モデルなら
        poi=e/2.d0/g-1.d0
c だけど木材みたいにE/G>10とかの材料でこれをやると
c ポアソン比が1を越えたりするので注意。
c つまりE/G>10とかの時点で等方性材料になり得ない。
        poi=0.d0
c せん断補正係数：pak
        pak=10.d0*(1.d0+poi)/(12.d0+11.d0*poi)
c せん断変形の影響を無視した答えが欲しいときは、
c  pak=pak*1.d10 みたいに1より極端に大きくする
c 軸長：ell
        ell=10.244d0
        ell=10.d0
c 断面の桁幅：b
        b=0.1d0
        b=5.d-2
c 断面の桁高：h
        h=0.5d0
        h=25.d-2
c 断面積：a
        a=b*h
c x軸（強軸）回りの断面二次モーメント：xi
        xi=b*h**3/12.d0
c Trahair の解と比較するために面内変位を無視するときは
c xi を仮想的にyiに対して極端に大きくする
c        xi=xi*1.d6
c そのようにxiが極端に大きい場合は面内にたわまないから
c その場合は曲げ面内(yz面）のせん断変形の影響も微少になるだろう
c
c y軸（弱軸）回りの断面二次モーメント：yi
        yi=h*b**3/12.d0
c 長方形断面のねじり定数:cj, そりねじり定数:wi
         call kukei(b,h,cj,wi)
         ea=e*a
         eix=e*xi
         eiy=e*yi
         gj=g*cj
         eiw=e*wi
         el=ell/float(nofe)
c
c 柱のオイラー座屈公式
         eul=(mode*pi/ell)**2*eiy
         write(*,*) ' オイラーの座屈解：',eul
c
         write(*,*)' 一発目の弧長を入れて下さい'
        read(*,*) r
        write(*,*)' 増分弧長を入れて下さい'
        read(*,*) rn
c  ### データの出力 ###
         open(7,file='oirasen.d')
         write(7,150) nofe
         write(7,150) mode
         write(7,300) gmode
         write(7,150) modeo
         do 20 i=1,nofe
          write(7,200) 0.d0, 0.d0, 0.d0
20       continue
         write(7,300) e
         write(7,300) g
         write(7,300) a
         write(7,300) xi
         write(7,300) yi
         write(7,300) cj
         write(7,300) wi
         write(7,300) el
         write(7,300) pak
         write(7,300) pak
c 比較対照解の種類：nhikak
c nhikak=1:圧縮単純梁のオイラー座屈の解
c nhikak=2：等曲げ円弧アーチの横ねじれ座屈のウラソフの解
c nhikak=3: 片持ち梁の横ねじれ座屈のトラヘアの解（面内変位無視）
            nhikak=1
               write(7,150) nhikak
c  載荷条件
c 載荷箇所は右端の1箇所
         write(7,150) 1
c 右端の節点 nofe+1 では
          write(7,150) nofe+1
c Fz に -1 を載荷
           write(7,400)  0.d0
           write(7,400)  0.d0
           write(7,400) -1.d0 
           write(7,400)  0.d0
           write(7,400)  0.d0
           write(7,400)  0.d0
           write(7,400)  0.d0
c 境界条件
c 両端の二箇所
         write(7,150) 2
c 左端では
          write(7,150) 1
c u,v,w,θzを拘束
           write(7,400)  0.d0
           write(7,400)  0.d0
           write(7,400)  0.d0
           write(7,400)  1.d0
           write(7,400)  1.d0
           write(7,400)  0.d0
           write(7,400)  1.d0
c 右端では
          write(7,150) nofe+1
c u,v,θzを拘束
           write(7,400)  0.d0
           write(7,400)  0.d0
           write(7,400)  1.d0
           write(7,400)  1.d0
           write(7,400)  1.d0
           write(7,400)  0.d0
           write(7,400)  1.d0
c θxもθyも拘束しなければ、普通は、弱軸回りで座屈すると思う。
c Ix=Iyみたいな断面の場合は、全節点で面外変位を拘束して
c (平面問題にして)座屈する方向を決めてやらないと解が出ないかも。
c そり(ねじれ率)を拘束すべきかどうかは、適宜 判断。
c
c
c 局所系で与えたい境界条件はなし
         write(7,150) 0
c
c 画面出力させたい要素
c 要素の左節点を出力させたいのが0箇所
         write(7,150) 0
c 要素の右節点を出力させたいのが1箇所
         write(7,150) 1
c 右端の要素
         write(7,150) nofe
c ファイル出力させたい要素（以下同様）
         write(7,150) 0
         write(7,150) 1
         write(7,150) nofe
c 弧長
         write(7,100) r
         write(7,100) rn
          close(7)
150      format(i3)
100      format(f15.5)
200      format(1p3d17.9)
300      format(1pd13.5)
400      format(1pd9.1)
500      format(4i2)
         stop
         end
c
c
c
c 長方形断面のねじり定数とそりねじり定数を求める
c 詳細は
c http://www.str.ce.akita-u.ac.jp/~gotou/programoj/#kukei
c 
          subroutine kukei(tan,tyou,cj,wi)
          implicit real*8 (a-h, o-z)
          dimension x(1000),y(1000)
c ねじりの応力関数、そり関数の級数の項数：kyuu
          kyuu=7
          a=tan/2.d0
          b=tyou/2.d0
c 短辺/2の分割数:na
c と長辺/2の分割数:nb（偶数）
          na=100
          nb=100
          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
       write(*,500) ' ねじり定数（Timoshenkoの長方形断面）：  ',cj
       write(*,500) ' ねじり定数（薄肉の近似：bt^3/3）：      ',bt33
       write(*,500) ' そりねじり定数（断面積分）：            ',wi
       write(*,500) ' そりねじり定数（薄肉近似：(bt)^3/144）：',
     &                (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
