import math
import numpy as np
import matplotlib.pyplot as plt

# определим константы

R = 4
kk = 0.065
alpha = 0.001
T = [0, 15]
c = 1.84
e = 0.00001
n = 100
#gamma = alpha * R / kk


def in_gamma(const_c, const_kk, const_h_t):
    return const_kk * const_h_t / const_c


def in_gamma1(const_c, const_kk, const_h_t, const_h_r, const_r_i):
    return 1 + const_kk * const_h_t * (2 * const_r_i + const_h_r) / ((const_h_r ** 2) * const_r_i * const_c)


def in_gamma2(const_kk, const_h_t, const_c, const_h_r):
    return const_kk * const_h_t / (const_c * (const_h_r ** 2))


def in_gamma3(const_c, const_kk, const_h_t, const_h_r, const_r_i):
    return const_kk * const_h_t / (const_c * (const_h_r ** 2)) + in_gamma(const_c, const_kk, const_h_t) / (
            const_r_i * const_h_r)


def method(h_t, h_r, t):
    r = np.arange(0, R + h_r, h_r)
    values = np.exp(-(r / (0.1 * R)) ** 2)
    coef = [[0] * 2] * (len(r) - 1)
    for time in np.arange(1, len(np.arange(0, t, h_t)), 1):
        for i in np.arange(0, len(r), 1):
            if i == 0:
                coef[i] = [1, 0]
                values[i] = np.exp(-(r[i] / (0.1 * R)) ** 2)
            elif i == len(r) - 1:
                values[i] = coef[i - 1][1] / (1 - coef[i - 1][0] + alpha * h_r / kk)
            else:
                if i == 1:
                    A1 = (-1) * in_gamma3(c, kk, h_t, h_r, r[i])
                    B1 = in_gamma1(c, kk, h_t, h_r, r[i]) - in_gamma2(kk, h_t, c, h_r)
                    C1 = 0
                    F1 = values[i]
                    Alpha1 = - (A1 / (B1 + C1 * coef[i - 1][0]))
                    Betta1 = ((F1 - C1 * coef[i - 1][1]) / (B1 + C1 * coef[i - 1][0]))
                    coef[i] = [Alpha1, Betta1]
                elif i == len(r) - 2:
                    AI = in_gamma1(c, kk, h_t, h_r, r[i]) * (1 + h_r * alpha / kk) - in_gamma3(c, kk, h_t, h_r, r[i])
                    BI = 0
                    CI = (-1) * in_gamma2(kk, h_t, c, h_r)
                    FI = values[i]
                    AlphaI = - (AI / (BI + CI * coef[i - 1][0]))
                    BettaI = ((FI - CI * coef[i - 1][1]) / (BI + CI * coef[i - 1][0]))
                    coef[i] = [AlphaI, BettaI]
                else:
                    A = (-1) * in_gamma3(c, kk, h_t, h_r, r[i])
                    B = in_gamma1(c, kk, h_t, h_r, r[i])
                    C = (-1) * in_gamma2(kk, h_t, c, h_r)
                    F = values[i]
                    Alpha = - (A / (B + C * coef[i - 1][0]))
                    Betta = ((F - C * coef[i - 1][1]) / (B + C * coef[i - 1][0]))
                    coef[i] = [Alpha, Betta]

        for i in np.arange(len(r) - 2, -1, -1):
            values[i] = coef[i][0] * values[i + 1] + coef[i][1]
    return [r, values]


def f(x0, gamma):
    I = (math.tan(x0) - gamma / x0)
    return I


def f1(x0, gamma):
    I = 1 / (math.cos(x0) * math.cos(x0)) + gamma / (x0 * x0)
    return I


def newtons_method(x0, e, gamma):
    x0 = float(x0)
    while True:
        x1 = x0 - (f(x0, gamma) / f1(x0, gamma))
        if abs(x1 - x0) < e:
            return x1
        x0 = x1


""

# def sum():
#     Psi = []
#     resultsum = []
#     for i in range(1, 5 * n):
#     x.append4(newtons_method(i, e, gamma))
#     nu = sorted(list(set(x)))
#
#     rm = np.arange(0, R, R / n)
#
#     for i in range(1, n):
#         y = np.exp((-1) * (rm[i] / (0.1 * R)) ** 2) * math.cos(nu[i] * rm[i] / R)
#         coef = 4. * nu[i] / (2. * nu[i] * R + R * math.sin(2. * nu([i]))
#         Int = np.trapz(rm, y)
#
#         Psi.append(Int * coef)
#     sumprom = 0
#
#     for r in np.arange(0, R, R / n):
#         for i in range(1, n):
#             sumprom+=Psi(i) * np.exp((kk / (c) * (nu(i) ** 2) / R) * t * (-1)) * math.cos(nu[i] * r / R)
#         resultsum.append(sumprom)
#         sumprom = 0
#     return [rm, resultsum]


[x, y] = method(0.2, 0.2, 10)
# [a, b] = sum()
plt.plot(x, y, color='green')
plt.grid(True)
plt.show(True)