import numpy as np
import matplotlib.pyplot as plt

# ============================================================
# 1. PROPIEDADES DEL ACERO 1080
# ============================================================
k = 51              # W/(m·K)
rho = 7800          # kg/m³
cp = 470            # J/(kg·K)
alpha = k / (rho * cp)

# ============================================================
# 2. GEOMETRÍA Y MALLA DE 9 NODOS (1D a lo largo del eje Z)
# ============================================================
L = 0.1             # Longitud de la probeta (100 mm)
N = 9               # Número de nodos
dz = L / (N - 1)    # Distancia entre nodos (12.5 mm)

# Condiciones iniciales
T = np.full(N, 800.0) # Temperatura inicial 800 °C
T_agua = 25.0         # Temperatura del agua

# ============================================================
# 3. METODOLOGÍA INVERSA: EL VALOR DE H
# ============================================================
# El profesor dijo que el artículo sugiere un valor de ~7000 a 9000.
# ¡Aquí está la magia! Puedes cambiar este número para calibrar.
h_agua = 8000.0     # W/m²K (Valor del artículo)

# Estabilidad
dt = 0.5            # segundos (Paso de tiempo)
tiempo_total = 600  # 10 minutos

# ============================================================
# 4. SIMULACIÓN (Curva Teórica)
# ============================================================
tiempos = []
T_teorica = []      # Guardaremos la temperatura del nodo 1 (base)

for t in range(int(tiempo_total/dt)):
    T_new = np.copy(T)
    
    # Nodo 1 (Base - Agua): Tiene convección y conducción
    # Usamos el Área (A) y el Volumen (V=A*dz). Como A se cancela, queda:
    T_new[0] = T[0] + (dt / (rho * cp * dz)) * (h_agua * (T_agua - T[0]) + k * (T[1] - T[0]) / dz)
    
    # Nodos Interiores (2 al 8): Solo conducción
    for i in range(1, N-1):
        T_new[i] = T[i] + (alpha * dt / dz**2) * (T[i-1] - 2*T[i] + T[i+1])
        
    # Nodo 9 (Tapa - Aire): Asumimos que el aire casi no enfría (adiabático)
    T_new[N-1] = T_new[N-2] 
    
    T = np.copy(T_new)
    tiempos.append(t)
    T_teorica.append(T[0]) # Guardamos la temperatura de la base

# ============================================================
# 5. DATOS EXPERIMENTALES (Del artículo/tesis)
# ============================================================
# Estos datos simulan los que obtendrías del termopar del artículo.
# (Puedes ajustar estos números con los de la tesis exacta si los tienes).
tiempo_exp = [0, 5, 10, 20, 30, 50, 80, 120, 200, 300, 500, 600]
temp_exp = [800, 600, 450, 300, 220, 150, 100, 80, 60, 50, 40, 35]

# ============================================================
# 6. GRÁFICA COMPARATIVA (La que pide el profesor)
# ============================================================
plt.figure(figsize=(10, 6))

# Curva Teórica (Calculada con Python)
plt.plot(tiempos, T_teorica, color='blue', linewidth=2, label='Curva Teórica (Python)')

# Curva Experimental (Datos del Artículo)
plt.plot(tiempo_exp, temp_exp, 'o--', color='red', linewidth=2, label='Curva Experimental (Artículo)')

plt.title('Metodología Inversa: Comparación Térmica (Acero 1080)', fontsize=14)
plt.xlabel('Tiempo (s)', fontsize=12)
plt.ylabel('Temperatura (°C)', fontsize=12)
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend()
plt.xlim(0, 600)
plt.ylim(0, 850)

plt.tight_layout()
plt.show()