import numpy as np
from math import exp, sqrt, pi
from scipy.integrate import quad
import matplotlib.pyplot as plt


"""
  Numerial integration by Imoto-Moriguchi-Takazawa (I.M.T.) formula
"""


# Integration parameters
xmin, xmax = -1.0, 1.0
nblock_max = 10

# define function to be integrated
def func(x):
    if x < -1.0 or 1.0 < x: return 0.0
    return sqrt(1.0 - x*x)

Q = 0.00702985840661; # int(exp(-1/t - 1/(1-t) by Gauss-Legendre method
#eps to aboid anomarly at x = 0 and 1.0. 
#  take 1.0/300 because the exponent of exp(x) is limited to approximately [-300, 300]
eps = 1.0 / 300.0  
def dphi_du(t):
    if t < eps or 1.0 - eps < t: return 0.0
    
    v = -1.0 / t - 1.0 / (1.0 - t)
    return exp(v) / Q;

def phi(u):
    ret = quad(dphi_du, 0.0, u)

    return ret[0]


#===================
# main routine
#===================
print("Numerical integration using by IMT formula")

fx = func(xmin)
#Fx = F(xmax) - F(xmin)
Fx = pi / 2.0
print("")
print("Analytical values:")
print(f"  f({xmin})={fx}")
print(f"  integ(f)[{xmin}, {xmax}]={Fx}")

print("")
print("Integration by area-devided Gauss-Legendre integration")
print("nBlock\tS\tError")
for nblock in range(2, nblock_max+1):
    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 = phi(ui)

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

    S *= (xmax - xmin) / nblock
    print(f"{nblock:5}\t{S:.8f}\t{S - Fx:12.6g}")


# 積分点と積分領域のグラフを描く
u_list_rough   = np.arange(0.0, 1.000001, 0.05)  # [0, 1]
x_list_rough = [phi(u) for u in u_list_rough]          # [0, 1]
x_list_rough = [-1.0 + 2.0 * x for x in x_list_rough]  # convert to [-1, 1]
y_list_rough = [func(x) for x in x_list_rough]

u_list   = np.arange(0.0, 1.000001, 0.02)  # [0, 1]
x_list = [phi(u) for u in u_list]          # [0, 1]
x_list = [-1.0 + 2.0 * x for x in x_list]  # convert to [-1, 1]
y_list = [func(x) for x in x_list]

# グラフ枠の作成
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, tight_layout = True, gridspec_kw={'hspace': 0})

# 上のグラフ
ax1.plot(x_list, u_list)
ax1.plot(x_list_rough, u_list_rough, linestyle = '', marker = 'o', markersize = 3.0)
# 下のグラフ
ax2.plot(x_list, y_list)
ax2.plot(x_list_rough, y_list_rough, marker = 'o', markersize = 3.0)
ax2.fill(x_list_rough, y_list_rough, color='skyblue', alpha=0.4)

ax1.set_ylabel('$u(x)$', fontsize=16)
ax1.tick_params(axis='both', which='major', labelsize=16)
# x目盛り文字列を削除
ax1.tick_params(axis='x', which='both', labelbottom=False) 
 
ax2.set_xlabel('$x$', fontsize=16)
ax2.set_ylabel('sqrt$(1-x^2)$', fontsize=16)
ax2.tick_params(axis='both', which='major', labelsize=16)

# 積分点について、下のグラフ枠から上のグラフ枠まで垂直線を引く
for x in x_list_rough:
    ax1.axvline(x, color='gray', linestyle='--', linewidth = 0.3)
    ax2.axvline(x, color='gray', linestyle='--', linewidth = 0.3)

plt.show()

