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

# # import module

# In[1]:


import numpy as np
import holoviews as hv

hv.extension('bokeh', logo=False)


# # parameter

# In[2]:


rhosw = 1.65 # 水の水中比重
g = 9.8 # 重力加速度


# # 沈降速度式

# \begin{align}
# w_0 &= \sqrt{sgd}F \\
# \ F &=  
# \left\{ 
#    \begin{array}{ll}
#      \sqrt{ \dfrac{2}{3} + \dfrac{36 \nu^2}{sgd^3}}  -\sqrt{\dfrac{36 \nu^2}{sgd^3}} &(  d \leq 1 \textrm{mm}) \\
#      \sqrt{ \dfrac{2}{3} } &(  d > 1\textrm{mm}) 
# \end{array} \right. 
# \end{align}
# 
# \begin{align}
# \nu:水の動粘性係数
# \end{align}
# 

# In[3]:


def rubeyw0(dd, nu=10**(-6)):
    tmp = 36.0*nu**2/rhosw/g/dd**3
    F = np.sqrt(2.0/3.0 + tmp) - np.sqrt(tmp) if dd < 0.001 else np.sqrt(2.0/3.0)
    
    return np.sqrt(rhosw*g*dd)*F


# # 誤差関数

#  - 近似式を使う。
# https://en.wikipedia.org/wiki/Error_function#Polynomial
# 
#  - scipyにも関数があるのでそっちでも良い。

# \begin{align}
# \mathrm{erf(x)} &=  
# \left\{ 
#    \begin{array}{ll}
#     1-T &   &(  x \geq 0) \\
#     T-1 &  &(x < 0) 
# \end{array} \right. 
# \end{align}
# 
# \begin{align}
# T &= t \exp( -x^2    
# -1.26551223
# +1.00002368t
# +0.37409196t^2
# +0.09678418t^3
# -0.18628806t^4
# +0.27886807t^5 \\
# & \qquad \qquad -1.13520398t^6
# +1.48851587t^7
# -0.82215223t^8
# +0.17087277t^9
#   ) \\
#  t & = \dfrac{1}{1 + 0.5|x|}
# \end{align}

# In[4]:


def erfc(xi):
    
    def erf(xi):
        x = np.abs(xi)
        
        t = 1/(1+0.5*x)
        
        tt = t*np.exp( -x**2    
        -1.26551223*t**0
        +1.00002368*t**1
        +0.37409196*t**2
        +0.09678418*t**3
        -0.18628806*t**4
        +0.27886807*t**5
        -1.13520398*t**6
        +1.48851587*t**7
        -0.82215223*t**8
        +0.17087277*t**9
          )
        
        return  1.0 - tt if xi >= 0.0 else tt-1.0
    
    return 1.0 - erf(xi)


# # 基準面濃度式：板倉・岸の式
# Kは河床形態による係数. 0.0018を使用

# \begin{align}
#     C_a &= K \left( \dfrac{\alpha_*}{1+s} \dfrac{u_*}{w_0} \dfrac{\Omega}{\tau_*}-1 \right) \\
#     \Omega &= \dfrac{\tau_*}{B_*} \dfrac{\displaystyle \int^\infty_{a'} \xi \dfrac{1}{\sqrt{\pi}}  e^{-\xi^2} d\xi}
#     {\displaystyle \int^\infty_{a'}  \dfrac{1}{\sqrt{\pi}}  e^{-\xi^2} d\xi} 
#      + \dfrac{\tau_*}{B_* \eta_0} -1  \\
#      a' &= \dfrac{B_*}{\tau_*} - \dfrac{1}{\eta_0}
# \end{align}
# <!--     \frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x} = 0  -->

# \begin{align}
# \alpha_*=0.14,B_*=0.143,\eta_0=0.5,K=0.0018
# \end{align}

# 
# \begin{align}
#     \int^\infty_{a'} \xi \dfrac{1}{\sqrt{\pi}} e^{-\xi^2} d\xi &= \dfrac{1}{ 2\sqrt{\pi}} e^{-a'^2} \\
#     \int^\infty_{a'}  \dfrac{1}{\sqrt{\pi}} e^{-\xi^2} d\xi &= \dfrac{1}{2} (1-\mathrm{erf}(a'))
# \end{align}

# 
# \begin{align}
# \mathrm{erf}:誤差関数
# \end{align}
# 

# In[5]:


def ca(dd, taus, w0, K = 0.0018):
    alfa = 0.14
    Bs = 0.143
    eta = 0.5
    
    def omega(taus):
        aprime = Bs/taus - 1.0/eta
        
        c1 = 0.5 / np.sqrt(np.pi) * np.exp(-aprime**2) 
        c2 = 0.5 * erfc(aprime)
        
        if c2 <= 0.0 :
            rr = 0.0
        else:
            rr =taus/Bs*c1/c2+ taus/Bs/eta - 1.0
        
        return rr 
    
    us = np.sqrt(rhosw*g*dd*taus)
    r = K*( alfa/(rhosw+1)*us/w0*omega(taus)/taus - 1.0 )
    
    return 0.0 if r < 0.0 else r


# ## 濃度分布を確認

# In[6]:


r = np.arange(-1.3,2.1,0.1)
taus = 10**r

d = 0.173/1000
w0 = rubeyw0(d)
us = np.sqrt(rhosw*g*d*taus)
aaa = [ ca(d, tt, w0) for tt in taus]
hv.Curve((w0/us, aaa)).options(logx=True, logy=True)


# # Rouse分布

# \begin{align}
#     \dfrac{C}{C_a} &= \left( \dfrac{h-z}{z} \dfrac{a}{h-a} \right)^Z \\
#     Z &= \dfrac{w_0}{\kappa u_*}
# \end{align}
# 

# 
# \begin{align}
#     \kappa:カルマン定数, a:基準面の河床から高さ(=0.05h)
# \end{align}

# ## 図化

# In[7]:


h = 5.0
ib = 1.0/1000.0


# In[8]:


def cfig(h,ib,d):
    us = np.sqrt(g*h*ib)
    
    w0 = rubeyw0(d)
    taus = us**2/rhosw/g/d
    
    kappa=0.4
    zeta = w0/kappa/us
    
    y = np.arange(0.05,1.0,0.01) * h
    a = 0.05*h
    x = (y-a)/(h-a)
    
    CbyCa = (a*(h-y)/y/(h-a))**zeta
    Ca = ca(d, taus, w0)
    
    return CbyCa*Ca*10**6, x

g1 = hv.Curve(cfig(h,ib,0.1/1000),label='0.1mm',kdims=['C']).options(logx=True, logy=False)
g2 = hv.Curve(cfig(h,ib,0.2/1000),label='0.2mm',kdims=['C']).options(logx=True, logy=False)
g3 = hv.Curve(cfig(h,ib,0.5/1000),label='0.5mm',kdims=['C']).options(logx=True, logy=False)
g4 = hv.Curve(cfig(h,ib,1.0/1000),label='1.0mm',kdims=['C']).options(logx=True, logy=False)
g5 = hv.Curve(cfig(h,ib,2.0/1000),label='2.0mm',kdims=['C']).options(logx=True, logy=False)

f1 = (g1*g2*g3*g4*g5).options(xlabel='土砂濃度[ppm]', ylabel='相対水深',width=500, legend_position='bottom_left').redim.range(y=(0,1.0))


# In[9]:


f1


# In[10]:


hv.save(f1,'suspendProf.html')