c harisen.f で等曲げを受ける円弧アーチの横ねじれ座屈を解かせるための
c 入力データを作るプログラム。enkosen を実行させると enkosen.d という
c データファイルができるから、harisen<enkosen.d と実行させればよい。
c
          program enko
        implicit real*8 (a-h,o-z)
        dimension alp(256),gmm(256),phi(256)
        write(*,*) '要素数を入れて下さい'
        read(*,*) nofe
c        write(*,*)'座屈モードを入れて下さい'
c        read(*,*) mode
          mode=1
c        write(*,*)'外乱を入れる場合は、その倍率を入れて下さい'
c        write(*,*)'(0 を入れた場合は座屈点を捜します)'
c        read(*,*) gmode
          gmode=0.d0
c        write(*,*)'モード計算のための出力が要りますか　要1／不要0 '
c        read(*,*) modeo  
          modeo=0
c        write(*,*)'正の曲げなら+1、負の曲げなら-1を入力'
c        read(*,*) porn
          porn=1.d0
c        write(*,*)'中心角を「度」で入れて下さい'
c        read(*,*) th
         th=0.d0
        pi=2.d0*asin(1.d0)
        th=th*pi/180.d0
        alp(1)=th/2.d0
        alp(2)=alp(1)-3.d0*th/2.d0/float(nofe)
        do 5 i=3,nofe-1
        alp(i)=alp(i-1)-th/float(nofe)
5       continue
        alp(nofe)=-th/2.d0
c  断面定数
c ヤング率：e
        e=200.d6
        e=1.d7 * 1.d3
c せん断弾性係数：g
        g=77.2d6
        g=7.d5 * 1.d3
c ポアソン比：poi
        poi=e/2.d0/g-1.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
        print *, ' 梁の長さを入れて下さい'
        read *, ell
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 ウラソフの解析解
        call vlasov(mode,ell,eix,eiy,gj,eiw,th,pi,vln,vlp)
         write(*,*)' 修正Vlasov（負） =',vln
         print '(a,1pd11.3)','             ML/EIy=',vln*ell/eiy
         write(*,*)' 修正Vlasov（正） =',vlp
         print '(a,1pd11.3)','             ML/EIy=',vlp*ell/eiy
        eix2=eix*1.d6
        call vlasov(mode,ell,eix2,eiy,gj,eiw,th,pi,vln,vlp)
         write(*,*)' Vlasov（負）=',vln
         print '(a,1pd11.3)','        ML/EIy=',vln*ell/eiy
         write(*,*)' Vlasov（正）=',vlp
         print '(a,1pd11.3)','        ML/EIy=',vlp*ell/eiy
         write(*,*)' 一発目の弧長を入れて下さい'
        read(*,*) r
        write(*,*)' 増分弧長を入れて下さい'
        read(*,*) rn
c  ### データの出力 ###
         open(7,file='enkosen.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) alp(i), 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=2
               write(7,150) nhikak
c  載荷条件
c 載荷箇所は両端の二箇所
         write(7,150) 2
c 左端の節点1では
          write(7,150) 1
c Mxに -1*porn を載荷
           write(7,400)  0.d0
           write(7,400)  0.d0
           write(7,400)  0.d0
           write(7,400) -1.d0*porn
           write(7,400)  0.d0
           write(7,400)  0.d0
           write(7,400)  0.d0
c 右端の節点 nofe+1 では
          write(7,150) nofe+1
c Mx に +1*porn を載荷
           write(7,400)  0.d0
           write(7,400)  0.d0
           write(7,400)  0.d0
           write(7,400) +1.d0*porn
           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 ただ、これだと面外回転の回転軸が鉛直軸回りになって
c しまって、Vlasovとかの座屈モードと一致しない。で、
c 局所系で与えたい境界条件
c 両端の二箇所
         write(7,150) 2
c 左端の要素で
          write(7,150) 1
c 左側の回転変位を
           write(7,500) 0,1,0,0
c 右端の要素で
          write(7,150) nofe
c 右側の回転変位を
           write(7,500) 0,0,0,1
c 局所系で与える（アーチの半径軸回りに回転するようにする）
c
c 画面出力させたい要素
c 要素の左節点を出力させたいのが一箇所
         write(7,150) 1
c 真ん中（より一つ右の要素）
         write(7,150) int(float(nofe)/2.d0)+1
c 要素の右節点を出力させたいのはなし
         write(7,150) 0
c ファイル出力させたい要素（以下同様）
         write(7,150) 1
         write(7,150) int(float(nofe)/2.d0)+1
         write(7,150) 0
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 修正Vlasov の解析解を求めるサブルーチン
        subroutine vlasov(n,blen,eix,eiy,gj,eiw,th,pi,vln,vlp)
        implicit real*8(a-h,o-z)
        if(th.lt.1.d-6) th=1.d-6
        r=blen/th
        f=(n*pi/blen)**2*eiw+gj
        a=eiy/eix-1.d0+(1.d0/eix-eiy/eix/eix)*f
        b=(eiy+(1.d0-2.d0*eiy/eix)*f)/r
        c=((n*pi/blen)**2-1.d0/r/r)*eiy*f
        d=b**2-4.d0*a*c
        vln=(-b+sqrt(d))/2.d0/a
        vlp=(-b-sqrt(d))/2.d0/a
        return
        end
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
