c hari.f で単純支持された柱のオイラー座屈を解かせるための
c 入力データを作るプログラム。oiraa を実行させると oiraa.d という
c データファイルができるから、hari<oiraa.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  
        pi=2.d0*asin(1.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
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)
c 柱のオイラー座屈公式
         eul=(mode*pi/ell)**2*eiy
         write(*,*) ' オイラーの座屈解：',eul
c
         write(*,*)' 一発目の弧長を入れて下さい'
        read(*,*) r
        write(*,*)' 増分弧長を入れて下さい'
        read(*,*) rn
c  ### データの出力 ###
         open(7,file='oiraa.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) 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=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 局所系で与えたい境界条件はなし
         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
