2020-10-23

Soil pressure coefficient considering earthquake by Mononobe-Okabe

Active soil pressure coefficient considering earthquake by Mononobe-Okabe can be obtained by following equation. (USACE EM 1110-2-2100 Appendix G)

$\begin{equation} P_a=K_a\cdot (1-k_v)\cdot \left\{\cfrac{1}{2}\cdot \gamma_s\cdot d^2+q\cdot d\right\} \end{equation}$

$\begin{equation} K_a=\cfrac{\cos^2 (\phi-\alpha-\theta)}{ \cos\theta\cdot \cos^2\alpha\cdot \cos(\alpha+\delta+\theta)\cdot \left\{ 1+\sqrt{\cfrac{\sin(\phi+\delta)\cdot \sin(\phi-\beta-\theta)}{\cos(\alpha+\beta+\theta)\cdot \cos(\alpha-\beta)}} \right\}^2 } \end{equation}$

$\begin{equation} \theta=\tan^{-1}\left(\cfrac{k_h}{1-k_v}\right) \end{equation}$

$\begin{align} & P_a & & \text{active soil pressure} & \quad & K_a & & \text{active poil pressure coefficient} \\ & \gamma_s & & \text{unit weight of soil} & \quad & d & & \text{soil depth} \\ & q & & \text{surcharge load} & \quad & \phi & & \text{internal friction angle of soil} \\ & \alpha & & \text{gradient of wall surface} & \quad & \beta & & \text{gradient of ground surface} \\ & \delta & & \text{friction angle of wall surface} & \quad & \theta & & \text{seismic inertia angle} \\ & k_h & & \text{horizontal seismic coefficient} & \quad & k_v & & \text{vertical seismic coefficient} \end{align}$

import numpy as np

kh=0.077 # horizontal seismic coefficient
#kh=0 # horizontal seismic coefficient
kv=0 # vertical seismic coefficient
phi=np.radians(30.0) # internal friction angle of soil
alpha=0 # gradient of wall surface
beta=0 # gradient of ground surface
delta=0 # friction angle of wall surface
theta=np.arctan(kh/(1-kv)) # seismic inertia angle

c1=np.cos(phi-alpha-theta)
c2=np.cos(theta) * (np.cos(alpha))**2 * np.cos(alpha+delta+theta)
c3=1 + np.sqrt(np.sin(phi+delta) * np.sin(phi-beta-theta) / np.cos(alpha+beta+theta) / np.cos(alpha-beta))
Ka=c1**2/c2/c3**2

print(Ka)

Calculation results are shown below.

Ka=0.381 # earthquake (kh=0.077, kv=0)
Ka=0.333 # normal (kh=0, kv=0)

That's all. Thak you.

2020-10-07

Python: calculation of required plate thickness of exposed type penstock against buckling pressure by Brent method

A critical buckling pressure of exposed type penstocks can be obtained by following equation.

$\begin{equation} p_k=\cfrac{2\cdot E_s}{1-\nu_s{}^2}\cdot \left(\cfrac{t}{D_0'}\right)^3 \end{equation}$

$p_k$	critical buckling pressure of penstock
$t$	plate thickness excluding corrosion allowance
$D_0'$	design external diameter
$E_s$	elastic modulus of penstock (=206,000MPa)
$\nu_s$	Poisson's modulus of penstock (=0.3)

usually, a corrosion allowance is set to $\epsilon=1.5 mm$ , and the design plate thickness including corrosion allowance $t_0$ and the design external diameter $D_0'$ can be calculated as follows.

$\begin{equation} t_0=t+\epsilon \qquad D_0'=D_0+t_0 \end{equation}$

Where, $D_0$ means the design internal diameter.

In case of exposed type penstocks, it is required to withstand the buckling pressure of 0.02 MPa which includes a safety factor of 1.5, and the plate thickness with this capacity shall be defined as a lower limit of plate thickness.

The equation to be solved by Brent method is shown below. The equation has a shape of $f=0$ .

$\begin{equation} f=p_k-\cfrac{2\cdot E_s}{1-\nu_s{}^2}\cdot \left(\cfrac{t}{D_0+t+\epsilon}\right)^3=0 \end{equation}$

Required two initial values for Brent method are given as t1=1.0, t2=50.0 in the program.

Program code is shown below.

import numpy as np
from scipy import optimize


def func(t,d0,eps,pk):
    Es=206000 # elastic modulus of steel
    po=0.3    # Poisson's ratio of steel
    f=pk-2*Es/(1-po**2)*(t/(d0+t+eps))**3
    return f


def main():
    d0=4000.0 # internal diameter of penstock
    eps=1.5   # corrosion alloowance
    pk=0.02   # critical buckling pressure
    t1=1.0    # initial value for Brent method
    t2=50.0   # initial value for Brent method
    tt=optimize.brentq(func,t1,t2,args=(d0,eps,pk))
    t0=np.ceil(tt+eps) # required plate thickness
    print('t + eps=',tt+eps)
    print('t0=',t0)


#==============
# Execution
#==============
if __name__ == '__main__': main()

The calculation results are shown below.

t + eps= 15.695548237490522
t0= 16.0

That's all. Thank you.

2020-09-26

Non-uniform flow analysis by Python: calculation of water surface frofile of settling basin

A program to calculate water surface profile of settling basin which has an open channel with varied width and varied invert level. Since the condition is subcritical flow, the calculation will be carried out from downstream section to upstream section.

Basic equation of non-uniform flow analysis is shown below.

$\begin{equation*} \left(\cfrac{Q^2}{2 g A_2{}^2}+h_2+z_2\right) - \left(\cfrac{Q^2}{2 g A_1{}^2}+h_1+z_1\right) = \cfrac{1}{2}\left(\cfrac{n_2{}^2 Q^2}{R_1{}^{4/3} A_1{}^2} + \cfrac{n_1{}^2 Q^2}{R_2{}^{4/3} A_2{}^2}\right)\cdot(x_2-x_1) \end{equation*}$

$Q$	Discharge (constant value)
$x_2, z_2, A_2, R_2, n_2, h_2$	Properties of upstream section (distance, invert level, flow section area, hydraulic radius, Manning's roughness coefficient, water depth)
$x_1, z_1, A_1, R_1, n_1, h_1$	Properties of downstream section (distance, invert level, flow section area, hydraulic radius, Manning's roughness coefficient, water depth)

Where, subscripu '1' is for downstream section and subscript '2' is for upstream section.

The equations for calculations of section area $A$ and hydraulic radius $R$ for rectangular cross section are shown below.

$\begin{equation*} A=b\cdot h \qquad R=\cfrac{b\cdot h}{b+2\cdot h} \end{equation*}$

where, $b$ means channel width, $h$ means water depth.

Since this case is a condition under subcritical flow, the conditions of downstream with subscript of '1' are all known values, and calculation is carried out to obtain the water depth of $h_2$ using the know conditions.

The section properties are defined as a function of x in def sec(x,h). Since the calculation model has many varied section, the part of definition of section properties (def sec(x,h)) is long. However, the core part for bi-section method is simple.

Two required initial values for bi-section method are defined as ha=3.0 and hb=7.0 in the function of def bisection(h1,x1,x2,q). These initial values should be changed appropriately, depending on the problem to be solved.

A function def sec(x,h) has a fuction to output section properties (z, A, R, n) from input data of distance (x) and water depth (h). When this function def sec(x,h) is changed, this program can be applied to a non-uniform flow analysis of open channel with various cross section.

A program is shown below.

# Non-Uniform Flow Analysis (Subcritical flow)
import numpy as np


def sec(x, h):
    # definition of section property
    # x  : distance
    # h  : water depth
    # zz : invert level
    # aa : secion area
    # rr : hydraulic radius
    # nn : Manning's roughness coefficient
    n0=0.014 # roughness coefficient (normal value)
    zz,aa,rr,nn=0,0,0,0
    if 0.0 <= x < 11.0:
        ds=11
        nn=n0
        z1=562.2; b1=4.0
        z2=560.5; b2=26.0
        zz=z1-(z1-z2)/ds*x
        bb=b1+(b2-b1)/ds*x
        aa=bb*h
        rr=bb*h/(bb+2*h)
    if 11.0 <= x < 19.0:
        ds=8.0
        nn=n0
        z1=560.5; b1=12
        z2=560.5; b2=12
        zz=z1-(z1-z2)/ds*(x-11)
        bb=b1+(b2-b1)/ds*(x-11)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 19.0 <= x < 47.0:
        ds=29.0
        nn=n0
        z1=560.5     ; b1=12.0
        z2=z1+ds*0.02; b2=12.0
        zz=z1-(z1-z2)/ds*(x-19)
        bb=b1+(b2-b1)/ds*(x-19)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 47.0 <= x < 62.0:
        ds=15.0
        nn=n0
        z1=561.06    ; b1=12.0
        z2=z1+ds*0.02; b2=6.0
        zz=z1-(z1-z2)/ds*(x-47)
        bb=b1+(b2-b1)/ds*(x-47)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 62.0 <= x < 69.0:
        ds=7.0
        nn=n0
        z1=561.36    ; b1=6.0
        z2=z1+ds*0.02; b2=6.0
        zz=z1-(z1-z2)/ds*(x-62)
        bb=b1+(b2-b1)/ds*(x-62)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 69.0 <= x < 69.5:
        ds=0.5
        nn=n0
        z1=561.5; b1=6.0
        z2=562.0; b2=6.0
        zz=z1-(z1-z2)/ds*(x-69)
        bb=b1+(b2-b1)/ds*(x-69)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 69.5 <= x < 73.5:
        ds=4.0
        nn=n0
        z1=562.0; b1=6.0
        z2=562.0; b2=6.0
        zz=z1-(z1-z2)/ds*(x-69.5)
        bb=b1+(b2-b1)/ds*(x-69.5)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 73.5 <= x < 74.0:
        ds=0.5
        nn=n0
        z1=562.0; b1=6.0
        z2=561.5; b2=6.5
        zz=z1-(z1-z2)/ds*(x-73.5)
        bb=b1+(b2-b1)/ds*(x-73.5)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 74.0 <= x < 77.0:
        ds=3.0
        nn=n0
        z1=561.5; b1=6.5
        z2=561.5; b2=6.5
        zz=z1-(z1-z2)/ds*(x-74.0)
        bb=b1+(b2-b1)/ds*(x-74.0)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 77.0 <= x < 96.625:
        ds=19.625
        nn=n0
        z1=561.5; b1=6.5
        z2=563.0; b2=6.5
        zz=z1-(z1-z2)/ds*(x-77.0)
        bb=b1+(b2-b1)/ds*(x-77.0)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 96.625 <= x <= 102.125:
        ds=5.5
        nn=n0
        z1=563.0; b1=6.5
        z2=563.0; b2=6.5
        zz=z1-(z1-z2)/ds*(x-96.625)
        bb=b1+(b2-b1)/ds*(x-96.625)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    return zz,aa,rr,nn


def func(h2,h1,x1,x2,q):
    g=9.8
    z1,a1,r1,n1=sec(x1,h1)
    z2,a2,r2,n2=sec(x2,h2)
    e2=q**2/2/g/a2**2+h2+z2
    e1=q**2/2/g/a1**2+h1+z1
    e3=0.5*(n1**2*q**2/r1**(4/3)/a1**2 + n2**2*q**2/r2**(4/3)/a2**2)*(x2-x1)
    return e2-e1-e3


def bisection(h1,x1,x2,q):
    ha=3.0 # lower initial value for bisection method
    hb=7.0 # upper initial value for bisection method
    for k in range(100):
        hi=0.5*(ha+hb)
        fa=func(ha,h1,x1,x2,q)
        fb=func(hb,h1,x1,x2,q)
        fi=func(hi,h1,x1,x2,q)
        if fa*fi<0: hb=hi
        if fb*fi<0: ha=hi
        #print(fa,fi,fb)
        if np.abs(hb-ha)<1e-10: break
    return hi


def main():
    g=9.8    # gravity acceleration
    q=42.0   # discharge
    # starting point (sub-critical flow)
    h1=3.9 # water depth at starting point
    x1=0.0   # origin of distance coordinate
    z1,a1,r1,n1=sec(x1,h1) # section property
    v1=q/a1  # flow velocity
    print('{0:>8s}{1:>8s}{2:>8s}{3:>8s}{4:>10s}{5:>8s}'.format('x','z','h','z+h','Energy','v'))
    print('{0:8.3f}{1:8.3f}{2:8.3f}{3:8.3f}{4:10.5f}{5:8.3f}'.format(x1,z1,h1,z1+h1,z1+h1+v1**2/2/g,v1))

    # calculation point
    #xp=np.arange(1,103,1)
    xp=np.array([11,19,47,62,69,69.5,73.5,74,77,96.625,102.125])

    # water level calculation to upstream direction
    for x2 in xp:
        h2=bisection(h1,x1,x2,q)
        z2,a2,r2,n2=sec(x2,h2)
        v2=q/a2
        print('{0:8.3f}{1:8.3f}{2:8.3f}{3:8.3f}{4:10.5f}{5:8.3f}'.format(x2,z2,h2,z2+h2,z2+h2+v2**2/2/g,v2))
        x1=x2 # distance
        h1=h2 # water depth


#==============
# Execution
#==============
if __name__ == '__main__': main()

Calculated results are shown below. A water level drop is larger than expected.

       x       z       h     z+h    Energy       v
   0.000 562.200   3.900 566.100 566.46982   2.692
  11.000 560.500   5.971 566.471 566.47522   0.293
  19.000 560.500   5.971 566.471 566.47523   0.293
  47.000 561.060   5.410 566.470 566.47528   0.323
  62.000 561.360   5.091 566.451 566.47541   0.687
  69.000 561.500   4.950 566.450 566.47553   0.707
  69.500 562.000   4.444 566.444 566.47554   0.788
  73.500 562.000   4.444 566.444 566.47562   0.788
  74.000 561.500   4.954 566.454 566.47563   0.652
  77.000 561.500   4.954 566.454 566.47567   0.652
  96.625 563.000   3.431 566.431 566.47615   0.942
 102.125 563.000   3.431 566.431 566.47634   0.942

That's all. Thank you.

2020-09-26

Python : Calculation of normal depth of open channel with rectangular cross section by 'scipy.optimize.brentq'

Calculation program of normal depth of open channel with rectangular cross sectionby Brent's method is introduced in this post. Basic equations for normal depth calculation of open channel with rectangular cross section are shown below.

$\begin{gather} Q=A\cdot v \\ v=\cfrac{1}{n}\cdot R^{2/3}\cdot i^{1/2} \\ R=\cfrac{b\cdot h}{b+2\cdot h} \end{gather}$

$Q$	(known) discharge
$b$	(known) width of rectangular open channe
$n$	(known) Manning's roughness coefficient
$i$	(known) gradient of channel invert
$h$	(unknown) normal depth

Following equation is solved for $h$ .

$\begin{equation} f=Q-\cfrac{b\cdot h}{n}\cdot \left(\cfrac{b\cdot h}{b+2\cdot h}\right)^{2/3}\cdot i^{1/2}=0 \end{equation}$

The critical depth calculation of $h_c$ is extra.

# normal depth and critical depth of rectangular cross section
import numpy as np
from scipy import optimize


def cal_hc(q,b):
    # critical depth
    g=9.8
    hc=(q**2/g/b**2)**(1/3)
    return hc


def func(h,q,b,n,i):
    f=q-b*h/n*(b*h/(b+2*h))**(2/3)*i**(1/2)    
    return f    


def main():
    q=42.0  # discharge
    b=4.0   # channel width
    n=0.014 # Manning's roughness coefficient
    i=0.001 # invert gradient
    
    h1=0.0
    h2=10.0
    hh=optimize.brentq(func,h1,h2,args=(q,b,n,i))

    print('hn=',hh) # normal depth
    
    hc=cal_hc(q,b)
    print('hc=',hc) # critical depth
    
#==============
# Execution
#==============
if __name__ == '__main__': main()

Calculation results are shown below.

hn= 3.866645305835682
hc= 2.2407023732785825

That's all. Thank you.

2020-09-20

Jupyter の ipynb ファイルの中身を見たい

2020.09.20投稿

必要性

Jupyter notebook 上でコードを書いていて、「あのとき作った Jupyter のコードをコピーしたい！」ということがよくある。その ipynb ファイルが、たまたま今開いている Jupyter のディレクトリの下にあればいいのだが、そうでない場合は、一度 Jupyter notebook を終了・再起動し所要の ipynb ファイルを開いてエディタにコピー、もとの作業に戻るため Jupyter notebook を再起動するという手間をかけていた。なんとかならないかと思いながら調べていたら、効率よく他の ipynb ファイルの中身を確認する以下の３つの方法があることがわかったので、紹介しておく。もしここに記載の方法を使ったことにより不都合なことが起きても、筆者は知りません。自分の責任においてトライしてください。

VS code を使う

こちらのサイトに方法が記載されている。

最強のipynb-viewerはどれだ！jupyter-notebookを起動せずに簡単にipynbファイルを開く方法

私は VS Cod を使う方法を選択し使ってみた。

Mac の Get info の Open with: で、VS code をデフォルトにしておけば、ipynb ファイルをダブルクリックすることで開くことができる。

Jupyter lab を使う

Jupyter Lab なら複数のタブを開ける！そこでさっそくインストールして調子を見ているところです。私のマシンではなんとなく動きがもっさりしている感じ。

複数のブラウザで Jupyter を開く

Jupyter notebook でも　これ（http://localhost:8888/?token=xxxxxx......）を url に入れてやると違うブラウザで Jupyter notebook を開ける。例えば、私の場合、default は Safari であるが、terminal に出てくるこれ（http://localhost:8888/?token=xxxxxx......）を Firefox や Chrome の url にいれることにより、それぞれのブラウザで notebook を開ける。

以　上

2020-09-10

LG 27インチ４Kディスプレイ購入

９月７日注文、９月８日到着で、LGの27inch 4K display を購入しました。￥57,980也。ちょっと奮発しすぎたか。

４月より在宅勤務しています。会社に行くのは月に２〜３回ですかね。しかし７月よりお客さんから設計の仕事をいただき、今はかなり忙しい。まあ色々検討してレポートを書いているのですが、ディスプレイ上の作業領域が「もっと欲しい！」ということで、サブディスプレイを購入したのでした。LGにした理由は、USB-Cポートがついていたからというのが大きいです。このポートがあると、ディスプレイとiMacをUSB-Cケーブル１本でつなぐことができます。最初は１万円くらいのやつでいいかと思っていたのですが、色々Webを見ながら考えていたら、「目も悪くなってきたし、せっかくなら４Kにするか！」と思い立ったのでした。実際には27インチで4Kだと、100%のscaleでは文字が小さくて読めません。なので200%とかに拡大しています。「拡大したら4Kの意味ないじゃん」という向きもありますが、画面のドットそのものは精細なので、文字はとてもきれいです。買った当初は縦置きにしようと思い、試しに縦置きで使ってみたのですが、上の方がとても遠く感じます。そこで通常の横置きに。図面などを見るにも視野が広くなり便利です。このディスプレイは横から縦へ回転させることができるので、縦置きにするのは、南北に長い地形図やGoogle earthを見るときだけとし、普段は横置きとすることにしました。

はじめは縦置きにしてみましたが、あまり使いよくない。 f:id:damyarou:20200910154018j:plain

基本は横置きです。図面を見る時など視野が広くて便利。 f:id:damyarou:20200910154049j:plain

Google Earthなどで縦に画像を表示させたい場合は便利。 f:id:damyarou:20200910154122j:plain

とりあえず持っている機械類の集合写真。この写真はiPhoneで撮っているので、かわいそうに、iPhoneは写っていません。正面が新しいディスプレイ、左側のWin機は会社支給のPC. f:id:damyarou:20200910154154j:plain

つなぎ方のイメージ図。iMac以外はHDMIでつないでいます。 f:id:damyarou:20200910154221p:plain

以　上

2020-08-17

Python 逆調整池運用解析修正版（Flood Routine）

逆調整池運用解析プログラムの修正版

関数FLOODの中の収束計算で、前回貯水池水位のみを既知としていたが、本来貯水池水位と放流量はペアなので、前回放流量も既知とした。プログラム作成時になぜ収束計算の中で前回流量を既知として確定していなかったのかは思い出せない。
関連出力は、データフレームにして一気に書き出すことにした。

以下最新版プログラム（2020.08.08）

# ==============================================
# Flood Routine Analysis
# ==============================================
import sys
import numpy as np
import pandas as pd
from scipy import interpolate


def FLOOD(iOFC,fnameR1,fnameR2,fnameR3,ram,qrel):
    nr,rh,rv=INP_HV(fnameR1)                   # reservoir H-V
    nd,ti,q_in,ELini,Qoini=INP_INFLOW(fnameR2) # inflow time history
    no,elof,qof=INP_OUTFLOW(fnameR3)           # outflow capacity

    pEL=np.zeros(nd,dtype=np.float64) # storage of reervoir water level
    pQo=np.zeros(nd,dtype=np.float64) # storage of outflow
    pDE=np.zeros(nd,dtype=np.float64) # storage of error
    pVO=np.zeros(nd,dtype=np.float64) # storage of reservoir volume
    pRA=np.zeros(nd,dtype=np.float64) # storage of ramp up rate
    pUD=np.zeros(nd,dtype=int)        # storage of indicator od water level
    
    Vin=np.sum(q_in)*(ti[2]-ti[1])/2 # total volume of inflow
    Qimax=np.max(q_in)
    # Initial setting
    elv=ELini
#    elv=WLINI(q_in[0],elof,qof,no)     # Initial reservoir water level
    EL=elv
    VOL=RET_V(nr,rh,rv,elv)           # Initial reservoir volume
    q_out=Qoini
    qout0=Qoini
    print('q_in0={0:10.3f}'.format(q_in[0]))
    print('elv_ini={0:10.3f}'.format(EL))
    print('q_ini={0:10.3f}'.format(q_out))
    i=0
    iUD=0
    pUD[i]=iUD
    pEL[i]=EL
    pQo[i]=q_out
    pDE[i]=EL-elv
    pVO[i]=VOL
    pRA[i]=0
    # Iterative calculation
    iUD=1        # Assume reservoir water level increase
    dh=0.0005    # Water level increment for searching equilibrium point
    eps=0.001    # Convergence criteria
    itmax=int(1.0/dh)*100
    Vf=RET_V(nr,rh,rv,60.0)
    Vz=0
    Vin2=0
    Vot2=0
    q1=q_out
    for i in range(0,nd-1):
        dt=ti[i+1]-ti[i]             # Time interval
        Qin=0.5*(q_in[i]+q_in[i+1])  # Average inflow rate during dt
        hh=0.0                       # Zero set of water level increment
        k=0
        j=0
        while 1:
            k=k+1
            j=j+1;
            hh=hh+float(iUD)*dh      # Update waler level increment
            elv1=EL
            elv2=EL+hh
            if iOFC==8: q2=OFC8(q1,qout0,elv2,elof,qof,no,ti[i],ram,qrel)
            Qout=0.5*(q1+q2)         # Average outflow rate during dt
            dS=(Qin-Qout)*dt*3600.0  # Storage accumulated during dt
            R=VOL+dS                 # Total reservoir volume
            elv=RET_H(nr,rh,rv,R)    # Reservoir water level at reservoir volume R
            if iUD==1 and j==10:  # 10 times trial assuming reservoir water level increase
                if EL+hh>elv:     # if reservoir water level is less than assumed water level,
                    iUD=-1        # assume water level decrease
                    hh=0.0
                    j=0
            if iUD==-1 and j==10: # 10 times trial assuming reservoir water level decrease
                if EL+hh<elv:     # if reservoir water level is greater than assumed water level,
                    iUD=1         # assume water level decrease
                    hh=0.0
                    j=0
            if np.abs(EL+hh-elv)<eps: break  # Judgement of convergence
            if itmax<k: break
        VOL=R    # Cumulative volume
        EL=EL+hh # Elevation
        q_out=q2 # Outflow
        q1=q2
        
        pUD[i+1]=iUD
        pEL[i+1]=EL
        pQo[i+1]=q_out
        pDE[i+1]=EL-elv
        pVO[i+1]=VOL
        pRA[i+1]=(pQo[i+1]-pQo[i])/(ti[i+1]-ti[i])

    hmax=max(pEL)-ELini
    sys.stdout.write('Time: %10.3f\n'% np.max(ti))
    sys.stdout.write('h   : %10.3f\n'% hmax)
    sys.stdout.write('Qin : %10.3f %10.3f\n'% (np.min(q_in),np.max(q_in)))
    sys.stdout.write('Qout: %10.3f %10.3f\n'% (np.min(pQo),np.max(pQo)))
    sys.stdout.write('EL  : %10.3f %10.3f\n'% (np.min(pEL),np.max(pEL)))
    sys.stdout.write('\n')

    # TWC of 80m downstream of Re-Regulating Dam
    fnameR='inp_twc.txt'
    data=np.loadtxt(fnameR,skiprows=0,usecols=(0,1))
    qq=data[:,0]
    ee=data[:,1]
    ff=interpolate.interp1d(qq,ee)
    el_ds=ff(pQo)
    
    df = pd.DataFrame({ 'i'  : np.arange(0,nd), # numper
                        'iUD': pUD,   # indicator for water level
                        'ti' : ti,    # time
                        'EL' : pEL,   # water level of reservoir
                        'dEL': pDE,   # error
                        'VOL': pVO,   # reservoir volume
                        'Qi' : q_in,  # inflow to reservoir
                        'Qo' : pQo,   # outflow from reservoir
                        'RA' : pRA,   # ramp up rate
                        'ELD': el_ds  # downstream water level
                      })
    df = df.set_index('i')
    return df 


def WLINI(q,elof,qof,no):
    # To obtain initial water level at oitflow discharge q
    for i in range(0,no-1):
        if qof[i]<=q and q<=qof[i+1]: break
    if qof[no-1]<q: i=no-2
    x1=qof[i]
    y1=elof[i]
    x2=qof[i+1]
    y2=elof[i+1]
    a=(y2-y1)/(x2-x1)
    b=(x2*y1-x1*y2)/(x2-x1)
    elv=a*q+b
    if q<0.0: elv=elof[0]
    return elv


def OFC8(q1,qout0,elv2,elof,qof,no,tii,ram,qrel):
    Qlim=300.0
    if   5.0<tii and tii<= 29.0 :qout0=qrel
    if  29.0<tii and tii<= 53.0 :qout0=qrel
    if  53.0<tii and tii<= 77.0 :qout0=qrel
    if  77.0<tii and tii<=101.0 :qout0=qrel
    if 101.0<tii: qout0=39.0
    q2=OFC(elv2,elof,qof,no)
    rq2=q2
    if q2<qout0: qout0=q2
    q2=qout0
    st=0.01
    ts1=0.0
    ts2=24.0
    ts3=48.0
    ts4=72.0
    ts5=96.0
    if ts1+st<=tii and tii<=ts1+5.0+st:
        if elv2<68.0:
            q2=q1+ram*st
            if rq2<q2: q2=rq2
            if Qlim<=q2: q2=Qlim
        else:
            q2=rq2
    if ts2+st<=tii and tii<=ts2+5.0+st:
        if elv2<68.0:
            q2=q1+ram*st
            if rq2<q2: q2=rq2
            if Qlim<=q2: q2=Qlim
        else:
            q2=rq2
    if ts3+st<=tii and tii<=ts3+5.0+st:
        if elv2<68.0:
            q2=q1+ram*st
            if rq2<q2: q2=rq2
            if Qlim<=q2: q2=Qlim
        else:
            q2=rq2
    if ts4+st<=tii and tii<=ts4+5.0+st:
        if elv2<68.0:
            q2=q1+ram*st
            if rq2<q2: q2=rq2
            if Qlim<=q2: q2=Qlim
        else:
            q2=rq2
    if ts5+st<=tii and tii<=ts5+5.0+st:
        if elv2<68.0:
            q2=q1+ram*st
            if rq2<q2: q2=rq2
            if Qlim<=q2: q2=Qlim
        else:
            q2=rq2

    return q2


def OFC(elv,elof,qof,no):
    # Free overflow discharge
    for i in range(0,no-1):
        if elof[i]<=elv and elv<=elof[i+1]: break
    if elof[no-1]<elv: i=no-2
    x1=elof[i]
    y1=qof[i]
    x2=elof[i+1]
    y2=qof[i+1]
    a=(y2-y1)/(x2-x1)
    b=(x2*y1-x1*y2)/(x2-x1)
    q=a*elv+b
    if elv<=elof[0]: q=0.0
    return q

def RET_V(nr,rh,rv,elv):
    # To obtain reservoir cumulative volume from the water level
    for i in range(0,nr-1):
        if rh[i]<=elv and elv<=rh[i+1]: break
    if rh[nr-1]<elv: i=nr-2
    x1=rv[i]
    y1=rh[i]
    x2=rv[i+1]
    y2=rh[i+1]
    a=(y2-y1)/(x2-x1)
    b=(x2*y1-x1*y2)/(x2-x1)
    v=(elv-b)/a
    return v

def RET_H(nr,rh,rv,v):
    # To obtain reservoir water level from cumulative volume
    for i in range(0,nr-1):
        if rv[i]<=v and v<=rv[i+1]: break
    if rv[nr-1]<v: i=nr-2
    x1=rv[i]
    y1=rh[i]
    x2=rv[i+1]
    y2=rh[i+1]
    a=(y2-y1)/(x2-x1)
    b=(x2*y1-x1*y2)/(x2-x1)
    elv=a*v+b
    return elv

def INP_HV(fnameR1):
    # Input reservoir H-V data
    fin=open(fnameR1,'r')
    dat=fin.readlines()
    fin.close()
    n=len(dat)
    text=dat[0]
    text=text.strip()
    text=text.split()
    vcoef=float(text[0])
    rh=np.zeros(n-1,dtype=np.float64)
    rv=np.zeros(n-1,dtype=np.float64)
    for i in range(0,n-1):
        text=dat[i+1]
        text=text.strip()
        text=text.split()
        rh[i]=float(text[0])
        rv[i]=float(text[1])*vcoef
    nr=len(rh)
    return nr,rh,rv

def INP_INFLOW(fnameR2):
    # Input time sequence of inflow
    fin=open(fnameR2,'r')
    dat=fin.readlines()
    fin.close()
    n=len(dat)
    text=dat[0]
    text=text.strip()
    text=text.split()
    ELini=float(text[0])
    Qoini=float(text[1])
    ti  =np.zeros(n-1,dtype=np.float64)
    q_in=np.zeros(n-1,dtype=np.float64)
    for i in range(0,n-1):
        text=dat[i+1]
        text=text.strip()
        text=text.split()
        ti[i]  =float(text[0])
        q_in[i]=float(text[1])
    nd=len(ti)
    return nd,ti,q_in,ELini,Qoini

def INP_OUTFLOW(fnameR3):
    # Input outflow capacity (water level-discharge curve)
    fin=open(fnameR3,'r')
    dat=fin.readlines()
    fin.close()
    n=len(dat)
    elof=np.zeros(n,dtype=np.float64)
    qof =np.zeros(n,dtype=np.float64)
    for i in range(0,n):
        text=dat[i]
        text=text.strip()
        text=text.split()
        elof[i]=float(text[0])
        qof[i] =float(text[1])
    no=len(elof)
    return no,elof,qof


def main():
    lout=['out_150.csv','out_175.csv','out_200.csv','out_225.csv','out_250.csv']
    lram=[150,175,200,225,250]
    lqrel=[89.0,87.1,82.9,81.5,80.7]
    for out,ram,qrel in zip(lout,lram,lqrel):
        iOFC=8
        fnameR1='inp_RHV.txt' # H-V data of reservoir
        fnameR2='inp_c00.txt' # Inflow time history (Hydrograph)
        fnameR3='inp_OF8.txt' # Outflow capacity
        fnameW = out # Output file name
        df=FLOOD(iOFC,fnameR1,fnameR2,fnameR3,ram,qrel)
        df.to_csv(fnameW, sep=",")
    

#============================    
if __name__ == '__main__':
    main()

以　上