#!/usr/bin/env python
# coding: utf-8

# # 河川シンポの江頭先生の式の展開

# ベルヌーイの定理
# <br>  
# $$
#     \dfrac{v_1^2}{2g} + h_1 + z_{b1} = \dfrac{v_2^2}{2g} + h_2 + z_{b2} + i \Delta x
# $$  
# <br>  
# を変形すると
# <br>  
# $$
# \begin{align}
#     \dfrac{z_{b2}  - z_{b1}}{h_1} &= 1- \dfrac{h_2}{h_1} + \dfrac{1}{2}Fr_1^2\left\{1-\left(\dfrac{v_2}{v_1}\right)^2\right\} - \dfrac{i \Delta x}{h_1}\\
#  Fr_1 &= \dfrac{v_1}{\sqrt{gh_1}}
# \end{align}
# $$  
# 
# 

#  - 流れの連続式
# $$
#  v_1 h_1 B_1 =  v_2 h_2 B_2 
# $$
# 
#  - 流れの運動方程式：マニング式
# $$
# \begin{align}
#  v_1  &= \dfrac{1}{n_1}i_1^{1/2}h_1^{2/3} \\
#  v_2  &= \dfrac{1}{n_2}i_2^{1/2}h_2^{2/3}
# \end{align}
# $$
# 
#  - 流砂の連続式
#   
# $$
#  q_{b1}B_1 = q_{b2}B_2  
# $$
# 
#  - 流砂の運動方程式：流砂量式
# 
# $$
# \begin{align}
#  q_{b1}  & \sim \tau_1^m \sim (h_1i_1)^m \\
#  q_{b2}  & \sim \tau_2^m \sim (h_2i_2)^m
# \end{align}
# $$

# $
# X=\dfrac{h_2}{h_1}
# ,Y=\dfrac{v_2}{v_1}
# ,Z=\dfrac{i_2}{i_1}
# ,\xi=\dfrac{B_2}{B_1}
# ,\eta=\dfrac{n_2}{n_1}
# $
# 
# として，上式を変形すると
# 
# 流れの連続式より
# 
# $
# X=\dfrac{1}{\xi Y}
# $
# 
# 流れの運動方程式より
# 
# $
# Y=\dfrac{1}{\eta}Z^{1/2} X^{2/3}
# $
# 
# 流砂の連続式，運動方程式より
# 
# $
# Z =  \dfrac{1}{X \xi^{1/m}}
# $
# 
# これらを$X,Y,Z$について解くと
# $$
# \begin{align}
#     X &= \eta^{\frac{6}{7}} \xi^{-\frac{3(2 m - 1)}{7 m}} \\
#     Y &= \eta^{-\frac{6}{7}} \xi^{- \frac{m + 3}{7 m}} \\
#     Z &= \eta^{-\frac{6}{7}} \xi^{\frac{2(3 m - 5)}{7 m} }
# \end{align}
# $$

# $X,Y$を変形したベルヌーイの式に代入すると，
# 
# \begin{align}
#     \dfrac{z_{b2}  - z_{b1}}{h_1} &= 1 - \eta^{\frac{6}{7}} \xi^{-\frac{3(2 m - 1)}{7 m}} 
#     + \dfrac{1}{2}Fr_1^2\left\{1- \eta^{-\frac{12}{7}} \xi^{- \frac{2(m + 3)}{7 m}} \right\} - \dfrac{i \Delta x}{h_1}\\
#  Fr_1 &= \dfrac{v_1}{\sqrt{gh_1}}
# \end{align}

# ---
# 一箇所だけ乗数が違う．2(m+3)が2(2m+3)でしたが，多分先生の間違い.

# # 図化
# 
# $$
# \begin{align}
# %     dz_* &= \dfrac{z_{b2}  - z_{b1}}{h_1} \\
#     m &= 1.5 \\
#     \eta &= \dfrac{n_2}{n_1}=1 \\
#     \Delta x &= 1000 \\
#     h_1 &= 3.0 \\
#     i &= 0.001 \\
#     Fr_1 &= 0.6
# %     \dfrac{z_{b2}  - z_{b1}}{h_1} &= 1 - \eta^{\frac{6}{7}} \xi^{-\frac{3(2 m - 1)}{7 m}} 
# %     + \dfrac{1}{2}Fr_1^2\left\{1- \eta^{-\frac{12}{7}} \xi^{- \frac{2(m + 3)}{7 m}} \right\} - \dfrac{i \Delta x}{h_1}\\
# %  Fr_1 &= \dfrac{v_1}{\sqrt{gh_1}}
# \end{align}
# $$
# 
# 縦軸 $\dfrac{z_{b2}  - z_{b1}}{h_1}$，横軸$\dfrac{B_2}{B_1}$として図化

# In[1]:


import holoviews as hv
import numpy as np
hv.extension('bokeh')


# In[2]:


m = 1.5
dx = 1000.0
i = 0.001
h1 = 3.0
Fr = 0.6
y = lambda x : 1.0 - x**(-3*(2*m-1)/7/m) + 0.5*Fr**2*(1.0 - x**(-2*(m+3)/7/m)) - i * dx / h1
x = np.arange(1.0,2.01,0.1)
hv.Curve((x,y(x)), kdims='$B2/B1$', vdims='(zb2-zb1)/h1').options(width=400,height=400)


# この条件だと，$B_2/B_1$が1.7くらいでレベルになり，それより大きいと逆勾配になる感じです．

# # 参考:式展開-sympy

# In[3]:


from sympy import *
init_printing()

eta = Symbol('\eta',real=True, positive=True)
xi = Symbol('\\xi',real=True, positive=True)
X = Symbol('X',real=True, positive=True)
Y = Symbol('Y',real=True, positive=True)
Z = Symbol('Z',real=True, positive=True)
m = Symbol('m',real=True, positive=True)


# In[4]:


Ans = solve( [X-1/xi/Y, Y-Z**0.5*X**(2/3)/eta, Z - 1/X/xi**(1/m)], [X,Y,Z] ) 


# In[5]:


simplify(Ans[0])