import sys
import numpy as np
from numpy import sin, cos, tan, sinh, cosh, tanh, exp, log, sqrt, pi
from pprint import pprint

#from tkutils import replaceif


xfrac = [
    [], # 0-point, not dfined
    [], # 1-point, not dfined
    [-1.0 / sqrt(3), + 1.0 / sqrt(3)],           # 2-point order
    [-sqrt(3.0/5.0), 0.0, +sqrt(3.0/5.0)],       # 3-point order
    [-0.861136311594052, -0.339981043584856,      
        +0.339981043584856, +0.861136311594052], # 4-point order
    [-0.906179845938663992797626878299, -0.538469310105683091036314420700, 0.0, 
        +0.538469310105683091036314420700, +0.906179845938663992797626878299] # 5-point order
    ]
weight = [
    [], # 0-point, not dfined
    [], # 1-point, not dfined
    [1.0, 1.0],                  # 2-point order
    [5.0/9.0, 8.0/9.0, 5.0/9.0], # 3-point order
    [0.347854845137453, 0.652145154862546,
        0.652145154862546, 0.347854845137453], # 4-point order
    [0.236926885056189087514264040719, 0.478628670499366468041291514835, 0.568888888888888888888888888888, 
        0.478628670499366468041291514835, 0.236926885056189087514264040719] # 5-point order
	]

def gauss_legendre(func, xmin, xmax, np):
    px = xfrac[np]
    pw = weight[np]
    S = 0.0
    n = len(px)
    xh = (xmin + xmax) / 2.0
    h = xmax - xmin
    for j in range(n):
        x = xh + px[j] * h / 2.0;
        S += pw[j] * func(x) * h;
    S *= 0.5
    return S, px, pw

Q = 0.00702985840661; # int(exp(-1/t - 1/(1-t) by Gauss-Legendre method
def IMTfunc(x):
    global Q
    return exp(-1.0 / x - 1.0 / (1.0-x)) / Q;

def imt(func, xmin, xmax, nblock, n_gl_order, ngl):
    global Q

    h = (1.0 - 0.0) / nblock;
    S = 0.0;
    for i in range(1, nblock):
# ui is in the converted range [0, 1]
        ui = i * h; 

# x' = phi(ui) is in [0,1]
        phiui = 0.0
        for j in range(ngl):
            hgl = (ui - 0.0) / ngl
            x0_gl = 0.0 + j * hgl
            x1_gl = x0_gl + hgl
            Sg, px, pw = gauss_legendre(IMTfunc, x0_gl, x1_gl, n_gl_order)
            phiui += Sg

# obtain the original x from x' = phi(ui)
        xi = xmin + (xmax - xmin) * phiui
        yi = func(xi);
        phi1 = exp(-1.0 / ui - 1.0 / (1.0 - ui)) / Q
        S += yi * phi1

    S *= (xmax - xmin) / nblock
    return S


def doubleexp(func, xmin, xmax, umin, umax, nblock):
    h = (umax - umin) / nblock;
    S = 0.0;
    for i in range(0, nblock+1):
# ui is in the converted range [-1, 1]
       ui = umin + i * h; 
# x' = phi(ui) is in [0,1]
       phiui = tanh(pi / 2.0 * sinh(ui))

# obtain the original x from x' = phi(ui)
       xi = xmin + (phiui + 1.0) / 2.0 * (xmax - xmin)
       yi = func(xi);
       coshsinh = cosh(pi / 2.0 * sinh(ui))
       coshnh = cosh(ui)
       phi1 = coshnh / (coshsinh*coshsinh)
       S += yi * phi1

    S *= pi / 2.0 * h * (xmax - xmin) / 2.0
    return S


def rhomberg(func, xmin, xmax, nhmax, eps, *, IsPrint = 0):
# create 2D array by specifing the size to be nhmax * nhmax
    S = np.empty((nhmax, nhmax), dtype=float)

    n = 1
    h0 = xmax - xmin
    y0 = func(xmin)
    y1 = func(xmax)
    S[0][0] = (y0 + y1) / 2.0 * h0;
    if IsPrint == 1:
        print("S00 = {}".format(S[0][0]))
    for k in range(1, nhmax+1):
        n *= 2
        h = h0 / n

# 台形則の計算
        S[k][0] = 0.0;
        for i in range(n):
            xi  = xmin + i * h;
            xi1 = xi + h;
            yi  = func(xi);
            yi1 = func(xi1);
            S[k][0] += (yi + yi1) / 2.0 * h;
            if IsPrint == 1:
                print("  S[{}][0] = {}".format(k, S[k][0]))

        for d in range(1, k+1):
            S[k][d] = (4.0**d * S[k][d-1] - S[k-1][d-1]) / (4.0**d - 1.0);
            if IsPrint == 1:
                 print("    S[{}][{}] = {}".format(k, d, S[k][d]))

        diff = S[k][k] - S[k-1][k-1]
        if IsPrint == 1:
            print("  diff = S[{}][{}] - S[{}][{}] = {} - {} = {}"
                    .format(k, k, k-1, k-1, S[k][k], S[k-1][k-1], diff))
        if abs(diff) < eps:
            return (S[k][k] + S[k-1][k-1]) / 2.0, diff, i

    return (S[k][k] + S[k-1][k-1]) / 2.0, diff, -1

def integ_rieman1(func, x, h):
    return func(x) * h

def rieman(func, x0, x1, ndiv):
    h = (x1 - x0) / ndiv
    S = sum([func(x0 + i * h) for i in range(ndiv)])
#    S = 0.0
#    for i in range(ndiv):
#        x = x0 + i * h
#        S += func(x)
    return S * h

def rieman_array(x0, dx, y, xmin, xmax):
    S = 0.0
    for i in range(len(y)):
        x = x0 + i * dx
        if x < xmin or xmax < x:
            continue
        S += y[i]
    return S * dx

def rieman_array_func(x, y):
    S = [0.0]
    for i in range(1, len(y)):
        S.append(S[i-1] + y[i] * (x[i] - x[i-1]))
    return S

def integ_trapezoid1(func, x, h):
    return (func(x) + func(x+h)) * 0.5 * h

def trapezoid(func, x0, x1, ndiv):
    h = (x1 - x0) / ndiv
    x = [x0 + i*h for i in range(ndiv+1)]
    S = sum([func(x[i]) + func(x[i+1]) for i in range(ndiv)])
    return S * 0.5 * h

def integ_simpson1(func, x, h):
    return (func(x) + 4.0 * func(x+0.5*h) + func(x+h)) / 6.0 * h

def integ_simpmson_by_list(y, h):
    n = len(y)
    ret = sum([4.0 * y[i] + 2.0 * y[i+1] for i in range(1, n-2, 2)])
    ret += y[0] + y[-1]
    ret *= dx / 3.0
    return ret

def simpson(func, x0, x1, ndiv, *, iPrint = 1):
    if ndiv % 2 != 0:
        ndiv0 = ndiv
        ndiv = int(ndiv / 2) * 2 + 2
        if iPrint == 1:
            print("")
            print("Warning in tkintegration.simpson:: " +
                  "ndiv must be an even: change from {} to {}".format(ndiv0, ndiv))
    h = (x1 - x0) / ndiv
    x = [x0 + i*h for i in range(ndiv+2)]
    S = sum([func(x[i]) + 4.0*func(x[i+1]) + func(x[i+2]) for i in range(0, ndiv, 2)])
#    h2 = h + h
#    S = 0.0
#    for i in range(0, ndiv, 2):
#        x = x0 + i * h
#        S += func(x) + 4.0 * func(x+h) + func(x+h2)
    return S * h / 3.0

def integ_simpson381(func, x, h):
    return (3.0*func(x) + 9.0*func(x+h/3.0) + 9.0*func(x+2.0*h/3.0) + 3.0*func(x+h)) / 24.0 * h

def simpson38(func, x0, x1, ndiv, *, iPrint = 1):
    if ndiv % 3 != 0:
        ndiv0 = ndiv
        ndiv = int(ndiv / 3) * 3 + 3
        if iPrint == 1:
            print("")
            print("Warning in tkintegration.simpson38:: " +
                  "ndiv must be a multiple of 3: change from {} to {}".format(ndiv0, ndiv))
    h = (x1 - x0) / ndiv
    x = [x0 + i*h for i in range(ndiv+3)]
    S = sum([3.0*func(x[i])   + 9.0*func(x[i+1]) 
           + 9.0*func(x[i+2]) + 3.0*func(x[i+3]) for i in range(0, ndiv, 3)])
#    h2 = h + h
#    h3 = h2 + h
#    S = 0.0
#    for i in range(0, ndiv, 3):
#        x = x0 + i * h
#        S += 3.0*func(x) + 9.0*func(x+h) + 9.0*func(x+h2) + 3.0*func(x+h3)
    return S * h / 8.0

def integ_bode1(func, x, h):
    return (14.0*func(x) + 64.0*func(x+h/4.0) + 24.0*func(x+h/2.0) 
          + 64.0*func(x+3.0*h/4.0) + 14.0*func(x+h)) / 180.0 * h

def integ_bode_by_list(y, h):
    n = len(y)
    ret = sum([64.0 * y[i] + 24.0 * y[i+1] + 64.0 * y[i+2] + 28.0 * y[i+3] for i in range(1, n-4, 4)])
    ret += 14.0 * (y[0] + y[-1])
    ret *= dx / 45.0
    return ret

def bode(func, x0, x1, ndiv, *, iPrint = 1):
    if ndiv % 4 != 0:
        ndiv0 = ndiv
        ndiv = int(ndiv / 4) * 4 + 4
        if iPrint == 1:
            print("")
            print("Warning in tkintegration.bode:: " +
                  "ndiv must be a multiple of 4: change from {} to {}".format(ndiv0, ndiv))
    h = (x1 - x0) / ndiv
    x = [x0 + i*h for i in range(ndiv+4)]
    S = sum([14.0*func(x[i])   + 64.0*func(x[i+1]) + 24.0*func(x[i+2])
           + 64.0*func(x[i+3]) + 14.0*func(x[i+4]) for i in range(0, ndiv, 4)])
#    h2 = h + h
#    h3 = h2 + h
#    h4 = h3 + h
#    S = 0.0
#    for i in range(0, ndiv, 4):
#        x = x0 + i * h
#        S += (14.0*func(x) + 64.0*func(x+h) + 24.0*func(x+h2)
#            + 64.0*func(x+h3) + 14.0*func(x+h4))
    return S * h / 45.0
