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


"""
  Numerial integration by double expotnential type conversion
"""


pi = 3.1415926535


# Integration parameters
range_type = "double inifite"
umin, umax = -2.0, 2.0
ndiv_max = 21


def func(x):
    if x < -1.0 or 1.0 < x: return 0.0
    return sqrt(1.0 - x*x)

def phi(u):
    return tanh(pi / 2.0 * sinh(u))

def dphi_du(u):
    v = cosh(pi / 2.0 * sinh(u))
    return pi / 2.0 * cosh(u) / v / v

#===================
# main routine
#===================
print("Numerical integration using by double exponential type conversion formula")
print(f"Integration range type: {range_type}")

print("")
print("Analytical values:")
Fx = pi / 2.0
print(f"  integ(f)={Fx}")

print("")
print(f"u range: {umin} - {umax}")
xmin = phi(umin)
xmax = phi(umax)
print(f"x range: {xmin} - {xmax}")

print("Integration by area-devided Gauss-Legendre integration")
print("ndiv\tS\terror")
for ndiv in range(3, ndiv_max + 1, 2):
    h = (umax - umin) / (ndiv - 1)
    u_list = np.arange(umin, umax + 1.0e-5, h)
    S = 0.0;
    for u in u_list:
        x = phi(u)
        y = func(x);
        w = dphi_du(u)
        S += y * w
    S *= h  
    print(f"{ndiv:5}\t{S:.8f}\t{S - Fx:12.6g}")


# 積分点と積分領域のグラフを描く
umin, umax = -2.0, 2.0
nrough = 41
h_rough = (umax - umin) / (nrough - 1)
u_list_rough = np.arange(umin, umax + 1.0e-5, h_rough)
x_list_rough = [phi(u) for u in u_list_rough]
y_list_rough = [func(x) for x in x_list_rough]

ndiv = 201
h = (umax - umin) / (ndiv - 1)
u_list   = np.arange(umin, umax + 1.0e-5, h)
x_list = [phi(u) for u in u_list]
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_between(x_list_rough, y_list_rough, 0.0, 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()

