import numpy as np
from numpy import sin, cos, tan, arccos, arcsin, arctan, sqrt, exp, log, pi

# SI unit definition by CGPM on Nov 18-20, 2014
h           = 6.62607015e-34   # Js
h_bar       = 1.05457203e-34   # Js
hbar        = h_bar
e           = 1.602176634e-19  # C
kB          = 1.380649e-23     # JK-1
NA          = 6.02214076e23    # mol-1

c           = 2.99792458e8     # m/s

pi          = 3.14159265358979323846
pi2         = pi + pi
torad       = 0.01745329251944 # rad/deg
todeg       = 57.29577951472   # deg/rad
basee       = 2.71828183
me          = 9.1093897e-31    # kg
mp          = 1.6726231e-27    # kg
mn          = 1.67495e-27      # kg
u0          = 4.0 * 3.14*1e-7; # Ns^2C^-2
e0          = 8.854418782e-12; # C^2/N/m^2
e2_4pie0    = 2.30711e-28      # Nm<sup>2</sup>";
a0          = 5.29177e-11      # m";
R           = 8.31451          # J/K/mol";
F           = 96485.3          # C/mol";
g           = 9.81             # m/s2";

G                = 6.67259e-11        #Nm2/kg2 
DayToSecond      = 60 * 60 * 24       #s
SecondToDay      = 1.0 / DayToSecond
AstronomicalUnit = 1.49597870e11      #m
AU               = AstronomicalUnit

HartreeToeV = 27.2116          # eV";
# $me * $e*$e*$e / (4.0*$pi*$e0)^2 / $h_bar/$h_bar;
HTV         = HartreeToeV;
#HTV =~ s/\s+.$//;
RyToeV      = HartreeToeV / 2.0
KToeV       = kB / e
eVToK       = e / kB
JToeV       = 1.0 / e
eVToJ       = e
Debye       = 3.33564e-30       # Cm, for dipole (D)

T0_absolute = 273.15

log2pi_2 = 0.5 * log(2.0 * pi)

# Bernoulli number
Bn = []
Bn.append(1.0)              # B0
Bn.append(-1.0 / 2.0)       # B1
Bn.append( 1.0 / 6.0)       # B2
Bn.append( 0.0)
Bn.append(-1.0 / 30.0)      # B4
Bn.append( 0.0)
Bn.append( 1.0 / 42.0)      # B6
Bn.append(0.0)
Bn.append(-1.0 / 30.0)      # B8
Bn.append(0.0)
Bn.append(5.0 / 66.0)       # B10
Bn.append(0.0)
Bn.append(691.0 / 2730.0)   # B12
Bn.append(0.0)
Bn.append( 7.0 / 6.0)       # B14
Bn.append(0.0)
Bn.append(-3617.0 / 510.0)  # B16


def acos(rad):
    return arccos(rad)

def asin(rad):
    return arcsin(rad)

def atan(rad):
    return arctan(rad)

def degcos(deg):
    return cos(torad * deg)

def degsin(deg):
    return sin(torad * deg)

def degtan(deg):
    return tan(torad * deg)

def degacos(rad):
    return todeg * arccos(rad)

def degasin(rad):
    return todeg * arcsin(rad)

def degatan(rad):
    return todeg * arctan(rad)

def cosh(x):
    return (exp(x) + exp(-x))/2.0

def sinh(x):
    return (exp(x) - exp(-x))/2.0

def tanh(x):
    return sinh(x) / cosh(x)

def ipower(x, n):
    xn = 1
    while n > 0:
        if n % 2 == 0:
            x *= x
            n >>= 1
        else:
            xn *= x
            n -= 1
    return xn


def eVTonm(E):
    if(E == 0.0):
        return 0.0
    return (h / (E * e)) * c * 1.0e9

def nmToeV(wl):
    if(wl == 0.0):
        return 0.0
    return (h / (wl * e)) * c * 1.0e9

def kayserToHz(k): # k in cm-1
    return c * k / 1.0e-2

def mod(val, max):
    if val > max:
        val -= int(val / max) * max
    elif val < 0.0:
        val -= (-1 + int(val / max)) * max
    return val;

def log10(x):
    return log(x) * 0.4342944819

def factorial(n):
    if n <= 1:
        return 1
    y = 1
    for i in range(2, n+1):
        y *= i
    return y


def Reduce01(x):
    while x < 0.0:
        x += 1.0
    while x >= 1.0:
        x -= 1.0
    return x

def reduce01(x):
    if x < 0.0:
        x = x - int(x) + 1.0

# Next conditional if should not be 'elif'
# if x == -1, -2, -3 etc, x - int(x) + 1.0 is 1.0. Need to correct it to 0
    if x >= 1.0:
        x = x - int(x)

def round(val, digit = 4):
#    if type(val) is not float:
#        print("tksci.tksci.Round: val [{}] must be float".format(val))
#        return None

    k = pow(10.0, digit)
    return int(val * k + 0.1) / k

def Round(val, digit = 4):
    return round(val, digit)

def max_none(x):
    m = -1.0e100
    for v in x:
        if v is not None and m < v:
            m = v
    return m

def min_none(x):
    m = 1.0e100
    for v in x:
        if v is not None and m > v:
            m = v
    return m

def normalize_none(l, Amin = 0.0, Amax = 1.0, vmin = None, vmax = None):
    if vmax is None:
        vmax = max_none(l)
    if vmin is None:
        vmin = min_none(l)
    if vmax - vmin == 0.0:
        vmax = vmin + 1.0

    for i in range(len(l)):
        if l[i] is None:
            continue

        l[i] = (l[i] - vmin) / (vmax - vmin) * (Amax - Amin) + Amin

    return l

def normalize(l, Amin = 0.0, Amax = 1.0, vmin = None, vmax = None):
    if vmax is None:
        vmax = max(l)
    if vmin is None:
        vmin = min(l)
    if vmax - vmin == 0.0:
        vmax = vmin + 1.0

    for i in range(len(l)):
        l[i] = (l[i] - vmin) / (vmax - vmin) * (Amax - Amin) + Amin

    return l

def IsFactor(I, idiv):
    if idiv <= 0:
        return 0
    res = I - int(I / idiv) * idiv
    if abs(res) < 1.0e-10:
        return 1
    return 0

def Factors(I):
    iMax = int(I / 2)
    e = [1]
    for ic in range(2, iMax+1):
        if IsFactor(I, ic):
            e.append(ic)
    e.append(I)
    return e

def Factorize(I):
    iMax = int(I / 2)
    check = np.empty(iMax+1)
    e = []
    ic = 2
    while 1:
        if I == 1:
            break
        if ic > iMax:
            break
        if check[ic]:
            ic += 1
            continue
        if IsFactor(I, ic):
            e.append(ic)
            I = I / ic
            continue
        j = 1
        while 1:
            if ic*j > I/2:
                break
            check[ic*j] = 1
            j += 1
            ic += 1

    if len(e) == 0:
        e.append(I)
    return e


def Gaussian(x, x0, whalf, A = None):
#A = 1/whalf * sqrt(ln2 / pi)
    if A is None:
        A = 0.469718639 / whalf
#a = whalf / sqrt(ln2)
    a = whalf / 0.832554611
    X = (x - x0) / a
    return A * exp(-X * X)

def Lorentzian(x, x0, whalf, A = None):
#A = 1/whalf/pi
    if A is None:
        A = 1.0 / whalf / pi
    X = (x - x0) / whalf
    return A / (1.0 + X * X)

def GaussLorentz(x, x0, whalf, C0 = 1.0, Gfraction = 0.5, Gwratio = 1.0, A = None):
    y = C0 * Gfraction * Gaussian(x, x0, whalf * Gwratio, A = A)
    y += C0 * (1.0 - Gfraction) * Lorentzian(x, x0, whalf, A = A)
    return y


def Ea_Arrhenius(xT, yP, eps = 1.0e-300, kTinv = 1000.0):
    nT = len(xT)
    xTinv = [kTinv / xT[iT] for iT in range(nT)]
    yEa   = []
    for iT in range(nT):
        if iT == 0:
            i0 = 0
            i1 = 1
        elif iT == nT - 1:
            i0 = nT - 2
            i1 = nT - 1
        else:
            i0 = iT - 1
            i1 = iT + 1

#        print("i=", iT, xTinv[i1], xTinv[i0], yP[i1], yP[i0])
        slope = (log(yP[i1]+eps) - log(yP[i0]+eps)) / (xTinv[i1] - xTinv[i0])
        Ea = -slope * kTinv * kB / e
#        print("i=", iT, xTinv[i1] - xTinv[i0], yP[i1], log(yP[i1]+eps), log(yP[i0]+eps), slope, Ea)

        yEa.append(Ea)

    return xTinv, yEa

def combination(N, k):
    return factorial(N) / factorial(k) / factorial(N - k)
    
def BernoulliNumber(n):
    if n == 0:
        return 1.0

    sum = 0.0
    for k in range(0, n):
        C = combination(n + 1, k)
        Bk = BernoulliNumber(k)
        sum += C * Bk
    return -1.0 / (n + 1.0) * sum

def loggamma(x):
    N = 8
    v = 1.0
    while x < N:
        v *= x
        x += 1.0

    w = 1.0 / x / x
    res = (((((((Bn[16] / 16.0 / 15.0) * w  
              + (Bn[14] / 14.0 / 13.0)) * w
              + (Bn[12] / 12.0 / 11.0)) * w 
              + (Bn[10] / 10.0 / 9.0)) * w
              + (Bn[8] / 8.0 / 7.0)) * w    
              + (Bn[6] / 6.0 / 5.0)) * w
              + (Bn[4] / 4.0 / 3.0)) * w    \
              + (Bn[2] / 2.0 / 1.0)
    
    return res / x + log2pi_2 - log(v) - x + (x - 0.5) * log(x)

# Gamma function: error << 1.0e-10
# 出典: C言語によるアルゴリズム辞典
def Gamma_C(x):
    if x < 0.0:
        return pi / (sin(pi * x) * exp(loggamma(1.0 - x)))

    return exp(loggamma(x))


# Gamma function: accuracy ~ 2e-5
def gamma_py(x):
    q = 1.0
    while x < 5.0:
        q *= x
        x += 1.0

    xs = x * x
    xx = x
    p = 0.0833333333333 / xx
    xx *= xs
    p -= 2.777777777778e-3 / xx
    xx *= xs
    p += 7.93651e-4 / xx
    y = (x - 0.5) * log(x) - x + p
    return 2.50663 * exp(y) / q

def gamma(x):
    return Gamma_C(x)

