---

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

# ДАННЫЕ ЛОГОСА, ПЕРЕНЕСЁННЫЕ С ФОТО
time_logos = np.array([
    100, 100, 200, 300, 400, 500,
    600, 700, 800, 900, 1000
], dtype=float)

temperature_logos = np.array([
    100.000000000000000,
    99.998452691423466,
    99.989951118414226,
    99.964225981261563,
    99.907271292470710,
    99.803841029423637,
    99.639949022593072,
    99.405060423169289,
    99.092645256521450,
    98.700242827031630,
    98.228901765152543
])

# ПАРАМЕТРЫ ЗАДАЧИ
R = 0.1
T0 = 100.0
a2 = 2.3e-5
r0 = R / 2

print(f"Характерное время охлаждения: {R**2/a2:.2f} с")

# АНАЛИТИЧЕСКОЕ РЕШЕНИЕ
z = jn_zeros(0, 100)
coeff = 2 * T0 * j0(z * r0 / R) / (z * j1(z))

def analytical_temperature(t):
    t = np.atleast_1d(np.asarray(t, dtype=float))
    return np.exp(-a2 * np.outer(t, z**2) / R**2) @ coeff

temperature_python = analytical_temperature(time_logos)

# ТАБЛИЦА СРАВНЕНИЯ
results = pd.DataFrame({
    "Time_s": time_logos,
    "LOGOS_C": temperature_logos,
    "Python_analytical_C": temperature_python
})

results["Absolute_error_C"] = np.abs(
    results["LOGOS_C"] - results["Python_analytical_C"]
)

results["Relative_error_percent"] = (
    100 * results["Absolute_error_C"]
    / np.maximum(np.abs(results["Python_analytical_C"]), 1e-12)
)

print("\nСРАВНЕНИЕ РЕЗУЛЬТАТОВ")
print(results.round(5).to_string(index=False))

print("\nСредняя абсолютная ошибка:",
      round(results["Absolute_error_C"].mean(), 5), "°C")

print("Максимальная абсолютная ошибка:",
      round(results["Absolute_error_C"].max(), 5), "°C")

# ГРАФИК СРАВНЕНИЯ ТЕМПЕРАТУР
plt.figure(figsize=(10, 6))
plt.plot(time_logos, temperature_logos, "o-", label="ЛОГОС")
plt.plot(time_logos, temperature_python, "s--",
         label="Python: аналитическое решение")
plt.xlabel("Время, с (если TimeStep указан в секундах)")
plt.ylabel("Температура, °C")
plt.title("Сравнение температуры в точке r = 0.05 м")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig("temperature_comparison.png", dpi=150)
plt.show()

# ГРАФИК ОШИБКИ
plt.figure(figsize=(10, 5))
plt.plot(time_logos, results["Absolute_error_C"], "o-")
plt.xlabel("Время, с")
plt.ylabel("Абсолютная ошибка, °C")
plt.title("Расхождение ЛОГОСа и аналитического решения")
plt.grid(True)
plt.tight_layout()
plt.savefig("absolute_error.png", dpi=150)
plt.show()

# СОХРАНЕНИЕ ТАБЛИЦЫ
results.to_csv("comparison_LOGOS_analytical.csv",
               index=False, encoding="utf-8-sig")

print("\nГотово! Сохранены:")
print("comparison_LOGOS_analytical.csv")
print("temperature_comparison.png")
print("absolute_error.png")