c hari.f で等曲げを受ける円弧アーチの横ねじれ座屈を解かせるための
c 入力データを作るプログラム。enko を実行させると enko.d という
c データファイルができるから、hari<enko.d と実行させればよい。
c
          program enko
        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  
        write(*,*)'正の曲げなら+1、負の曲げなら-1を入力'
        read(*,*) porn
        write(*,*)'中心角を「度」で入れて下さい'
        read(*,*) th
        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 x軸回りの曲げ試験で測定されるヤング率：ex
        ex=200.d6
c y軸回りの曲げ試験で測定されるヤング率：ey
        ey=ex
c ねじり試験で測定されるせん断弾性係数：g
        g=77.2d6
c 軸長：ell
        ell=10.244
c 断面積：a
        a=9.288d-2
c x軸（強軸）回りの断面二次モーメント：xi
        xi=1.1363d-4
c 座屈前の面内変位の影響を無視したVlasovの解と比較するときは
c xi=xi*1.d6 などとして仮想的に面内剛性を極端に大きくする
c
c y軸（弱軸）回りの断面二次モーメント：yi
        yi=3.871d-5
c ねじり定数：ej
        ej=5.89d-7
c そりねじり定数：wi
        wi=5.5869d-7
c
         eix=ex*xi
         eiy=ey*yi
         ea=a*sqrt(ex*ey)
         eiw=wi*sqrt(ex*ey)
         gj=g*ej
         el=ell/float(nofe)
        call vlasov(mode,ell,eix,eiy,gj,eiw,th,pi,vln,vlp)
         write(*,*)' 修正Vlasov（負） =',vln
         write(*,*)' 修正Vlasov（正） =',vlp
        eix2=eix*1.d6
        call vlasov(mode,ell,eix2,eiy,gj,eiw,th,pi,vln,vlp)
         write(*,*)' Vlasov（負）=',vln
         write(*,*)' Vlasov（正）=',vlp
         write(*,*)' 一発目の弧長を入れて下さい'
        read(*,*) r
        write(*,*)' 増分弧長を入れて下さい'
        read(*,*) rn
c  ### データの出力 ###
         open(7,file='enko.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) ea
         write(7,300) eix
         write(7,300) eiy
         write(7,300) gj
         write(7,300) eiw
         write(7,300) el
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