---

import numpy as np
import matplotlib.pyplot as plt


# ==========================================================
# 1. ИСХОДНЫЕ ДАННЫЕ ИЗ ЗАДАНИЯ
# ==========================================================

# Таблица температур
T_table = np.array([2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000], dtype=float)

# Вариант 1
k_table_1 = np.array([
    8.200e-03,
    2.768e-02,
    6.560e-02,
    1.281e-01,
    2.214e-01,
    3.516e-01,
    5.248e-01,
    7.472e-01,
    1.025e+00
], dtype=float)

# Вариант 2
k_table_2 = np.array([
    1.600e+00,
    5.400e+00,
    1.280e+01,
    2.500e+01,
    4.320e+01,
    6.860e+01,
    1.024e+02,
    1.458e+02,
    2.000e+02
], dtype=float)

# Параметры для отладки из задания
R = 0.35       # см
T_w = 2000.0   # K
T_0 = 10000.0  # K
p = 4.0


# ==========================================================
# 2. ФУНКЦИИ T(r), u_p(r), k(T)
# ==========================================================

def T_of_r(r):
    z = r / R
    return (T_w - T_0) * z**p + T_0


def u_p_of_T(T):
    return 3.084e-4 / (np.exp(4.799e4 / T) - 1.0)


def u_p_of_r(r):
    return u_p_of_T(T_of_r(r))


def k_of_T_log_interp(T, T_nodes, k_nodes):
    """
    Интерполяция в логарифмах:
    xi = ln(T), eta = ln(k)
    """
    xi_nodes = np.log(T_nodes)
    eta_nodes = np.log(k_nodes)

    xi = np.log(T)
    eta = np.interp(xi, xi_nodes, eta_nodes)

    return np.exp(eta)


def make_k_function(k_table):
    def k_of_r(r):
        T = T_of_r(r)
        return k_of_T_log_interp(T, T_table, k_table)
    return k_of_r


# ==========================================================
# 3. СИСТЕМА УРАВНЕНИЙ ДЛЯ МЕТОДА СТРЕЛЬБЫ
# ==========================================================

def system(r, y, k_func):
    """
    y[0] = u
    y[1] = F

    Система:
        du/dr = -3 k(r) F
        dF/dr = k(r) (u_p(r) - u) - F/r
    """
    u = y[0]
    F = y[1]

    k = k_func(r)
    up = u_p_of_r(r)

    du = -3.0 * k * F
    dF = k * (up - u) - F / r

    return np.array([du, dF])


# ==========================================================
# 4. РУНГЕ-КУТТА 4 ПОРЯДКА
# ==========================================================

def rk4_step(fun, r, y, h, k_func):
    k1 = fun(r, y, k_func)
    k2 = fun(r + h / 2, y + h * k1 / 2, k_func)
    k3 = fun(r + h / 2, y + h * k2 / 2, k_func)
    k4 = fun(r + h, y + h * k3, k_func)

    return y + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0


# ==========================================================
# 5. РЕШЕНИЕ ЗАДАЧИ КОШИ ДЛЯ ЗАДАННОГО chi
# ==========================================================

def solve_for_chi(chi, k_table, n=4000):
    k_func = make_k_function(k_table)

    r = np.linspace(0.0, R, n + 1)
    h = r[1] - r[0]

    # Начальное значение
    up0 = u_p_of_r(0.0)
    u0 = chi * up0

    # Из условия симметрии:
    # F(0)=0
    # Для старта берём малое eps = h и используем асимптотику:
    # F(eps) ≈ 0.5 * k(0) * (u_p(0) - u(0)) * eps
    k0 = k_func(0.0)
    eps = h
    F_eps = 0.5 * k0 * (up0 - u0) * eps

    y = np.zeros((n + 1, 2))
    y[0] = [u0, 0.0]
    y[1] = [u0, F_eps]

    for i in range(1, n):
        y[i + 1] = rk4_step(system, r[i], y[i], h, k_func)

    u = y[:, 0]
    F = y[:, 1]
    up = u_p_of_r(r)

    # Невязка правого граничного условия:
    # F(R) = 0.39 * u(R)
    residual = F[-1] - 0.39 * u[-1]

    return r, u, F, up, residual


# ==========================================================
# 6. ПОДБОР chi МЕТОДОМ ДЕЛЕНИЯ ПОПОЛАМ
# ==========================================================

def find_chi(k_table, chi_left, chi_right, eps=1e-8):
    _, _, _, _, f_left = solve_for_chi(chi_left, k_table)
    _, _, _, _, f_right = solve_for_chi(chi_right, k_table)

    if f_left * f_right > 0:
        print("На концах интервала невязка одного знака.")
        print("Нужно выбрать другой интервал для chi.")
        return None

    while abs(chi_right - chi_left) > eps:
        chi_mid = (chi_left + chi_right) / 2
        _, _, _, _, f_mid = solve_for_chi(chi_mid, k_table)

        if f_left * f_mid <= 0:
            chi_right = chi_mid
            f_right = f_mid
        else:
            chi_left = chi_mid
            f_left = f_mid

    return (chi_left + chi_right) / 2


# ==========================================================
# 7. РАСЧЁТ И ГРАФИКИ
# ==========================================================

def solve_variant(k_table, variant_name, chi_left, chi_right):
    chi = find_chi(k_table, chi_left, chi_right)

    if chi is None:
        return

    r, u, F, up, residual = solve_for_chi(chi, k_table)

    print("==========================================")
    print(variant_name)
    print("chi =", chi)
    print("невязка =", residual)
    print("u(0) =", u[0])
    print("u(R) =", u[-1])
    print("F(R) =", F[-1])
    print("0.39*u(R) =", 0.39 * u[-1])
    print("==========================================")
    print()

    plt.figure(figsize=(9, 5))
    plt.plot(r, u, label="u(r)")
    plt.plot(r, up, label="u_p(r)")
    plt.plot(r, F, label="F(r)")
    plt.xlabel("r, см")
    plt.ylabel("значение")
    plt.title(f"Задача 2, {variant_name}")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()


# ==========================================================
# 8. ЗАПУСК
# ==========================================================

# Для варианта 1
solve_variant(k_table_1, "Вариант 1", 0.1, 0.3)

# Для варианта 2
solve_variant(k_table_2, "Вариант 2", 0.9, 1.0)