---

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

# ==========================================
# 1. ФИЗИЧЕСКИЕ ПАРАМЕТРЫ ЗАДАЧИ
# ==========================================
R = 0.1             # Радиус цилиндра, м (10 см)
a2 = 2.3e-5         # Температуропроводность стали, м^2/с
T0 = 100.0          # Начальная температура, °C
r0 = R / 2          # Точка замера (R/2 = 5 см от центра), м

# ==========================================
# 2. МАТЕМАТИЧЕСКОЕ ЯДРО (ФОРМУЛА 3)
# ==========================================
# Находим первые 1000 корней уравнения J0(z) = 0 для высокой точности
z_n = jn_zeros(0, 1000)

def get_temperature(t_val):
    """Аналитическое решение по формуле (3)"""
    T = 0.0
    for z in z_n:
        # Формула (3) из методички
        term = (1.0 / (z * j1(z))) * j0(z * (r0 / R)) * np.exp(-((a2 * z**2) / R**2) * t_val)
        T += term
    return 2 * T0 * T

# ==========================================
# 3. ДАННЫЕ ИЗ ЛОГОС
# ==========================================
# Берем данные из твоей таблицы 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
])

# ==========================================
# 4. РАСЧЕТ ТЕОРЕТИЧЕСКОЙ КРИВОЙ
# ==========================================
# Время от 0 до 1000 секунд (как в таблице ЛОГОСа)
time_python = np.linspace(0, 1000, 300)
# Считаем температуру в каждой точке
temp_python = [get_temperature(t) for t in time_python]

# ==========================================
# 5. ПОСТРОЕНИЕ ГРАФИКОВ
# ==========================================
plt.figure(figsize=(10, 6))

# Строим аналитическую кривую (синяя линия)
plt.plot(time_python, temp_python, 'b-', linewidth=2, label='Аналитическое решение (сталь, a²=2.3e-5)')

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

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

# Фиксируем оси, чтобы график выглядел красиво
plt.xlim(-50, 1050)
plt.ylim(99.0, 100.2)

plt.tight_layout()
plt.show()

# Дополнительный вывод в консоль
print(f"Характерное время tau = {R**2 / a2:.1f} сек")
print(f"Теоретическая температура при t=1000 сек: {get_temperature(1000):.3f} °C")
print(f"Температура из ЛОГОСа при t=1000 сек: {temp_logos[-1]:.3f} °C")