---

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

# -----------------------------
# Параметры задачи
# -----------------------------
CSV_PATH = "test.csv"

R = 0.1                 # м
a2 = 2.3e-5             # м^2/с, это a^2
T0 = 100.0              # °C
r0 = R / 2.0            # м, точка R/2

# ВАЖНО: в вашем CSV время, похоже, в миллисекундах.
# Если время в секундах, поставьте TIME_SCALE = 1.0
TIME_SCALE = 1e-3

# Число корней Бесселя
N = 300
z = jn_zeros(0, N)      # положительные корни J0(z)=0

# -----------------------------
# Аналитическая формула
# -----------------------------
def T_analytic(t, r=r0):
    """
    Аналитическое решение для бесконечного цилиндра:
    T(t,r)=2*T0*sum( J0(z_n*r/R)/(z_n*J1(z_n)) *
                     exp(-a2*z_n^2*t/R^2) )
    """
    if abs(t) < 1e-15:
        return T0
    terms = (
        j0(z * r / R)
        / (z * j1(z))
        * np.exp(-a2 * z**2 * t / R**2)
    )
    return 2.0 * T0 * np.sum(terms)

# -----------------------------
# Чтение CSV из ANSYS
# -----------------------------
def load_ansys_csv(path):
    rows = []
    with open(path, "r", encoding="utf-8-sig", newline="") as f:
        reader = csv.reader(f, delimiter=";")
        header = next(reader)

        for row in reader:
            if len(row) < 2:
                continue

            s_time = row[0].strip().replace(",", ".")
            s_temp = row[1].strip().replace(",", ".")

            if not s_time or not s_temp:
                continue

            rows.append([float(s_time), float(s_temp)])

    # В вашем файле первая строка, вероятно, соответствует t=0,
    # но в столбце времени стоит 100. Исправляем это.
    if len(rows) >= 2:
        if abs(rows[0][0] - rows[1][0]) < 1e-12 and abs(rows[0][1] - T0) < 1e-6:
            rows[0][0] = 0.0

    return rows

rows = load_ansys_csv(CSV_PATH)

# -----------------------------
# Сравнение
# -----------------------------
print("t_csv    t_phys,s    T_ANSYS     T_analytic   error       rel.error")
print("-" * 80)

results = []

for t_csv, T_ansys in rows:
    t = t_csv * TIME_SCALE
    T_an = T_analytic(t)

    err = T_ansys - T_an
    rel = err / T_ansys if abs(T_ansys) > 1e-12 else np.nan

    results.append([t_csv, t, T_ansys, T_an, err, rel])

    print(
        f"{t_csv:8.3f} {t:10.6f} "
        f"{T_ansys:11.6f} {T_an:12.6f} "
        f"{err:11.3e} {rel:11.3e}"
    )

# -----------------------------
# Сохранение результата сравнения
# -----------------------------
with open("comparison.csv", "w", encoding="utf-8-sig", newline="") as f:
    writer = csv.writer(f, delimiter=";")
    writer.writerow([
        "TimeStep_csv",
        "Time_phys_s",
        "T_ANSYS",
        "T_analytic",
        "error",
        "relative_error"
    ])
    writer.writerows(results)

# -----------------------------
# График
# -----------------------------
times = [row[1] for row in results]
T_ansys_list = [row[2] for row in results]
T_analytic_list = [row[3] for row in results]

plt.figure(figsize=(8, 5))
plt.plot(times, T_ansys_list, "o-", label="ANSYS")
plt.plot(times, T_analytic_list, "s--", label="Analytical")
plt.xlabel("t, s")
plt.ylabel("T, °C")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig("comparison_plot.png", dpi=200)
plt.show()