---

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import j0, j1, jn_zeros
from scipy.optimize import curve_fit

# ==========================================
# 1. ПАРАМЕТРЫ ЗАДАЧИ (из методички)
# ==========================================
R = 0.1             # Радиус цилиндра, м
T0 = 100.0          # Начальная температура, °C
r0 = R / 2          # Точка замера, м
# Берем 1000 корней, чтобы линия начиналась ровно с 100, а не с 99.4
z_n = jn_zeros(0, 1000) 

# ==========================================
# 2. ДАННЫЕ ИЗ ЛОГОС (твой Excel)
# ==========================================
time_logos = np.array([0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000])
temp_logos = np.array([
    100.0, 
    99.999618305476005, 
    99.997252127994145, 
    99.989257146927443, 
    99.969682582904980, 
    99.930723936601183, 
    99.863634228334291, 
    99.759914934029354, 
    99.612150859932925, 
    99.414628053968144, 
    99.163586548299634
])

# ==========================================
# 3. МАТЕМАТИЧЕСКАЯ ФУНКЦИЯ (Формула 3)
# ==========================================
def get_temperature(t_val, a2_val):
    T = 0.0
    for z in z_n:
        term = (1.0 / (z * j1(z))) * j0(z * (r0 / R)) * np.exp(-((a2_val * z**2) / R**2) * t_val)
        T += term
    return 2 * T0 * T

# ==========================================
# 4. ПОДГОНКА ПАРАМЕТРА (чтобы графики совпали)
# ==========================================
# Мы говорим питону: "Найди такое значение a2, чтобы кривая прошла через точки ЛОГОСа"
# Начальное приближение ставим очень маленьким (1e-10), так как ЛОГОС остывает медленно
popt, pcov = curve_fit(
    lambda t, a2_fit: [get_temperature(ti, a2_fit) for ti in t], 
    time_logos, 
    temp_logos, 
    p0=[1e-9] 
)
a2_fitted = popt[0]

print(f"Табличная a2 (сталь) = 2.30e-05")
print(f"Эффективная a2 (которую использовал ЛОГОС) = {a2_fitted:.2e}")

# ==========================================
# 5. ПОСТРОЕНИЕ ГРАФИКОВ
# ==========================================
# Создаем плавную линию от 0 до 1000
time_python = np.linspace(0, 1000, 500)
# Считаем её с ПОДОГНАННЫМ параметром
temp_python_fitted = [get_temperature(t, a2_fitted) for t in time_python]

plt.figure(figsize=(10, 6))

# Идеальная аналитическая кривая (подогнанная под ЛОГОС)
plt.plot(time_python, temp_python_fitted, 'b-', linewidth=2, label='Аналитическое решение (Формула 3)')

# Твои точки из ЛОГОС
plt.plot(time_logos, temp_logos, 'ro', markersize=6, label='Данные ЛОГОС')

plt.title('Остывание цилиндра: Аналитика (Python) и ЛОГОС')
plt.xlabel('Время, t (шаг/секунды)')
plt.ylabel('Температура, T (°C)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.legend()

plt.tight_layout()
plt.show()

print(f"Температура при t=1000 по ЛОГОС: {temp_logos[-1]:.3f} °C")
print(f"Температура при t=1000 по формуле: {get_temperature(1000, a2_fitted):.3f} °C")