"""
FTCS (Forward-Time Centered-Space) explicit finite difference solver
for the 1D heat equation with a movable point heat source:

    dT/dt = alpha * d²T/dx²  +  Q(t) * alpha / (k * dx) * delta(x - x_src)

Boundary conditions (both ends Robin / convective):
  - x=0  (Robin): convective exchange with ambient
  - x=L  (Robin): convective exchange with ambient

Returns T_matrix of shape (NX, NT).
"""

import sys
import os
import numpy as np

sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
import config


def solve():
    """Run the FTCS solver and return (T_matrix, x_vals, t_vals).

    Returns
    -------
    T_matrix : np.ndarray, shape (NX, NT)
    x_vals   : np.ndarray, shape (NX,)
    t_vals   : np.ndarray, shape (NT,)
    """
    alpha  = config.ALPHA
    k      = config.K
    L      = config.L
    T0     = config.T0
    Q_val  = config.Q_VAL
    t_step = config.T_STEP
    h      = config.H_CONV
    T_amb  = config.T_AMB
    T_end  = config.T_END
    NX     = config.NX
    NT     = config.NT
    x_src  = config.X_SRC

    x_vals = np.linspace(0.0, L, NX)
    t_vals = np.linspace(0.0, T_end, NT)
    dx = x_vals[1] - x_vals[0]
    dt = t_vals[1] - t_vals[0]

    r = alpha * dt / dx**2
    if r > 0.5:
        print(
            f"[FDM WARNING] Stability condition violated: "
            f"r={r:.4f} > 0.5  (dt={dt:.6g}, dx={dx:.6g}). "
            "Solution may diverge."
        )

    i_src = int(np.argmin(np.abs(x_vals - x_src)))

    robin_coeff = dx * h / k

    T_matrix = np.zeros((NX, NT), dtype=np.float64)
    T_matrix[:, 0] = T0

    T_cur = np.full(NX, T0, dtype=np.float64)

    for n in range(NT - 1):
        t_now = t_vals[n]

        T_new = T_cur.copy()

        # Interior FTCS (nodi 1 .. NX-2)
        T_new[1:-1] = T_cur[1:-1] + r * (T_cur[2:] - 2.0 * T_cur[1:-1] + T_cur[:-2])

        # Robin BC a x=0
        T_new[0] = (T_cur[1] + robin_coeff * T_amb) / (1.0 + robin_coeff)

        # Robin BC a x=L
        T_new[-1] = (T_cur[-2] + robin_coeff * T_amb) / (1.0 + robin_coeff)

        # Sorgente puntuale al nodo i_src — applicata dopo le BCs per non essere sovrascritta
        Q_now = Q_val if t_now >= t_step else 0.0
        T_new[i_src] += Q_now * alpha * dt / (k * dx)

        T_cur = T_new
        T_matrix[:, n + 1] = T_cur

    return T_matrix, x_vals, t_vals