#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 29 10:43:43 2025

@author: louisjeanmart
"""


#NYST chemical database

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from scipy.optimize import newton
from CoolProp.CoolProp import PropsSI
from scipy.signal import savgol_filter
import pandas as pd




np.set_printoptions(threshold=np.inf) 
#------------------Database-----------

# Charger le fichier Excel
df = pd.read_excel('/Users/louisjeanmart/Desktop/Memoire/resultats_oxydation.xlsx')  # Remplace par le bon chemin

# Transformer chaque colonne en liste
List_tga = []
for col in df.columns:
    
    List_tga.append(df[col].tolist())

Fe_540_Xtga = List_tga[0]
Fe_540_Ytga = List_tga[4]
Fe_540_Ytga2 = List_tga[5]
        
        
 


def enthalpy_Fe(T):
    
    
    t = T/1000
    
    # A = 23.97449
    # B = 8.367750
    # C = 0.000277
    # D = -0.000086
    # E = -0.000005
    # F = 0.268027
    # G = 62.06336
    # H = 0 
    # return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H

    if T >= 298 and T < 700:
        
        A = 18.42868
        B = 24.64301
        C = -8.913720
        D = 9.664706
        E = -0.012643
        F = -6.573022
        G =  42.51488
        H = 0.000000
        
        
        return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H
    

    
    if T >= 700 and T <= 1042:
        
        A = -57767.65
        B = 137919.7
        C = -122773.2
        D = 38682.42
        E = 3993.080
        F = 24078.67
        G =  -87364.01
        H = 0
        
        
        return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H
 
    
    if T > 1042 and T < 1100:
        
        A = -325.8859
        B = 28.92876
        C = 0
        D = 0
        E = 411.9629
        F = 745.8231
        G =  241.8766
        H = 0
        
        
        return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H
    
    
    if T >= 1100 and T < 1809:
        
        A = -776.7387
        B = 919.4005
        C = -383.7184
        D = 57.08148
        E = 242.1369
        F = 697.6234
        G = -558.3674
        H = 0
        
        
        return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H
    
    if T >= 1809 and T < 3133.345:
        
        A = 46.02400
        B = -1.884667*1e-8
        C = 6.094750*1e-9
        D = -6.640301*1e-10
        E = -8.246121*1e-9
        F = -10.80543
        G = 72.54094
        H = 0
        
        
        return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H
    

def enthalpy_Fe0(T):
    
    t = T/1000
    
    if T >= 298 and T < 1650:
        
        A = 45.75120
        B = 18.78553
        C = -5.952201
        D = 0.852779
        E = -0.081265
        F = -286.7429
        G =  110.3120
        H = -272.0441
        
        
        return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H
    
    
    if T >= 1650 and T < 5000:
        
        A = 68.19920
        B = -4.501232*1e-10
        C = 1.195227*1e-10
        D = -1.064302*1e-11
        E = -3.092680*1e-10
        F = -281.4326
        G =  137.8377
        H = -249.5321
        
        
        return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H


def enthalpy_Fe3O4(T):
    
    t = T/1000
    
    if T >= 298 and T < 900:
        
        A = 104.2096
        B = 178.5108
        C = 10.61510
        D = 1.132534
        E = -0.994202
        F = -1163.336
        G = 212.0585
        H = -1120.894
    
        return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H
    
    if T >= 900 and T < 3000:
        
        A = 200.8320
        B = 1.586435*10-7
        C = -6.661682*10-8
        D = 9.452452*10-9
        E = 3.186020*10-8
        F = -1174.135
        G = 388.0790
        H = -1120.894
    
        return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H
    
def enthalpy_O2(T):
    
    t = T/1000
    
    if T >= 100 and T < 700:
        
        A = 31.32234
        B = -20.23531
        C = 57.86644
        D = -36.50624
        E = -0.007374
        F = -8.903471
        G = 246.7945
        H = 0
   
    if T >= 700 and T < 2000:
        
        A = 30.03235
        B = 8.772972
        C = -3.988133
        D = 0.788313
        E = -0.741599
        F = -11.32468
        G = 236.1663
        H = 0
        
    if T >= 2000 and T < 6000:
        
        A = 20.91111
        B = 10.72071
        C = -2.020498
        D = 0.146449
        E = 9.245722
        F = 5.337651
        G = 237.6185
        H = 0
    
    return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H


def enthalpy_SiO2(T):
    
    t = T/1000
    
    if T >= 298 and T < 1996:
        
        A = 72.77482
        B = 1.293543
        C = -0.004360
        D = 0.000798
        E = -4.140645
        F = -941.1364
        G = 114.4370
        H = -905.4887
    
        return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t + F #- H
    
    


Matrice_H = []
T_liste = np.arange(300,1995,1) #rafiner ca 


for i in range (len(T_liste)):
    h_liste = [T_liste[i],enthalpy_Fe(T_liste[i]),enthalpy_Fe0(T_liste[i]),enthalpy_Fe3O4(T_liste[i]),enthalpy_O2(T_liste[i]),enthalpy_SiO2(T_liste[i])]
    #modifier la division par les masse mmolaire
    Matrice_H.append(h_liste)



#-------------------------------------
pi = np.pi
R = 8.314

####param initialisation

Tg = 1219#1222



#Source papier
#k0_FeO = 0.03862188626860003 #grow rate parabolic FeO
k0_FeO = 2.3245865931379544*1e21



k0_Fe3O4 = 2.616009223322519*1e20 #grow rate parabolic Fe3O4
#Ta_FeO = 44912.60066810567 #jsp FeO (je pense energie d'activation)
Ta_FeO = 96589.02989590589


Ta_Fe3O4 = 95076.94922500897  


rho_FeO = 5745 #density Fe0 [kg/m3]
rho_Fe3O4 = 5170 #density  Fe3O4 [kg/m3]
rho_Fe = 7874
rho_SiO2 = 2196 # ou 2648 si quartz

rho_fe7Si = 0.07*2330 + 0.93*7874

#valeur intial
rp_0 = 10*1e-6
#rp_0 = 31.5*1e-6 #pour poudre Mathieu

X_0 = 10*1e-9 #A modifier  





#rapport massique stochio
M_Si = 28.0855/1000
M_FeO = (55.845 + 15.999)/1000
M_O2 = (15.999*2)/1000
M_Fe = (55.845)/1000
M_FeO3 = 3*M_Fe + (4*15.999)/1000
M_SiO2 = M_Si + M_O2

nu_Fe_FeO = M_Fe/M_FeO
nu_Fe_Fe3O4 = 3*M_Fe/M_FeO3
nu_02_SiO2 = M_O2/M_Si

m0 = 4/3*pi*rp_0**3 * rho_fe7Si

m_Fe = m0 - 0.07*m0
m_tot = m0 + nu_02_SiO2*0.07*m0
m_SiO2 = m_tot - m_Fe



r_Fe = (m_Fe/(4/3*pi*rho_Fe)) **(1/3)
X_SiO2 = m_SiO2/(4*pi*r_Fe**2*rho_SiO2)


Ap = 4 * pi * (r_Fe + X_SiO2)**2 #Aire jusqu'a couche FeO




#### Second EDO bilan d'énergie

nu_O2_FeO = M_O2/M_FeO #rapport massique O2/FeO stoich
nu_O2_Fe3O4 = 2*M_O2/M_FeO3 #rapport massique O2/Fe2O3 stoich

#radiative heat losses
epsilon = 0.88 #emissivity
sigma = 5.670374419 * 1e-8 #boltzman cst

#Covective heat losses
lambda_air = 6.763*1e-2  #considéré cst
h_p = 2*lambda_air/(20*1e-6) #Re considéré = 0
air_elem = ['N2','O2']
air_mol =  [0.79,0.21]


#####Energie release
q_Fe3O4 = 4.841*1e6 #source papers
q_FeO = 3.787*1e6


h_O2 = enthalpy_O2(Tg)*1000


W_O2 = M_O2 #molar mass = molar weight pas en unité mais en valeur 

Matrice_h = []


debit = 46*1e-3/60
speed = debit/(0.13**2*pi)
#print(speed)



lambda_air = 0
Pr_air = 0
Re_air = 0
rho_air = 0
cp_air = 0
mu_air = 0

for i in range(len(air_elem)):
    
    lambda_air += (PropsSI('L', 'T', 298, 'P', 101325, air_elem[i])*air_mol[i]) 
    lambda_air += 1/(air_mol[i]/PropsSI('L', 'T', 298, 'P', 101325, air_elem[i]))
    lambda_air *=0.5
    Pr_air += (PropsSI('PRANDTL', 'T', 298, 'P', 101325, air_elem[i])*air_mol[i])
    rho_air += (PropsSI('D', 'T', 298, 'P', 101325, air_elem[i])*air_mol[i])
    cp_air += (PropsSI('C', 'T', 298, 'P', 101325, air_elem[i])*air_mol[i])
    mu_air += (PropsSI('V', 'T', 298, 'P', 101325, air_elem[i])*air_mol[i])

   # print(rho_air)
Re_air = speed*rho_air*2*rp_0/mu_air

Nu = 2 + 0.552*Re_air**0.5*Pr_air**(1/3)
    
h_test = Nu*lambda_air/(20*1e-6)

def dHp_dt(dm_FeO_dt, dm_O2_dt, Tp, Tg, Ap, q_FeO, h_O2, W_O2, h_p, sigma, epsilon):
   
    
    term_chimique = q_FeO * dm_FeO_dt  - (h_O2 / W_O2) * dm_O2_dt
    term_convection = h_p * (Tp - Tg)
    term_radiation = sigma * epsilon * (Tp**4 - Tg**4)
    
    test_perte = h_test*(Tp - 298)
    #term_radiation = 0 #negligeable
    
    #print(term_convection)
    
    return term_chimique - Ap * (term_convection + term_radiation ) #- Ap*test_perte




Tp = 300


H0 = (m_Fe/M_Fe*Matrice_H[Tp-300][1] + m_SiO2/M_SiO2*Matrice_H[Tp-300][5])*1000#+mFeO_0*Matrice_H[Tp-300][2]+m_Fe3O4_0*Matrice_H[Tp-300][3]





test = []
X=np.zeros(2)
X1 = []
X2 = []
Rayon= []


def reaction_odes(t, m, params,isoT,Tg):
    # Initialisation des valeurs initiales
    (nu_Fe_FeO, rho_FeO, k0_FeO, Ta_FeO) = params
    
    
   
    
    if isoT == True:
        T_guess = 1173
        
    else:
        
        if t==0:
            T_guess = 300
        else:
        
            T_guess = test[-1]
        
   
     
    Tp = T_guess

   
    r_Fe = (m[1]/(rho_Fe*4/3*pi))** (1/3) #me parait bon mais a voir
    Rayon.append(r_Fe)
    X[0] = m[0]/(4*pi*((r_Fe+X_SiO2)**2)*rho_FeO)
    A_FeO = 4 * np.pi * (r_Fe + X_SiO2 + X[0])**2     
 
    
    
    X1.append(X[0])
    #X2.append(X[1])
    
    
    dm_FeO_dt = (rho_FeO * A_FeO * k0_FeO / X[0]) * np.exp(-Ta_FeO / Tp)    
    dm_Fe_dt = -nu_Fe_FeO * dm_FeO_dt 
    dm_O2_dt = -nu_O2_FeO * dm_FeO_dt
    
    
    
    Tsf = (2*Tp + Tg)/3
    lambda_air = 0
    Pr_air = 0
    Re_air = 0
    rho_air = 0
    cp_air = 0
    mu_air = 0
    for i in range(len(air_elem)):
        
        lambda_air += (PropsSI('L', 'T', Tsf, 'P', 101325, air_elem[i])*air_mol[i]) 
        lambda_air += 1/(air_mol[i]/PropsSI('L', 'T', Tsf, 'P', 101325, air_elem[i]))
        lambda_air *=0.5
        Pr_air += (PropsSI('PRANDTL', 'T', Tsf, 'P', 101325, air_elem[i])*air_mol[i])
        rho_air += (PropsSI('D', 'T', Tsf, 'P', 101325, air_elem[i])*air_mol[i])
        cp_air += (PropsSI('C', 'T', Tsf, 'P', 101325, air_elem[i])*air_mol[i])
        mu_air += (PropsSI('V', 'T', Tsf, 'P', 101325, air_elem[i])*air_mol[i])
        

    Re_air = speed*rho_air*0.13/mu_air
    
    h_p = 2*lambda_air/(20*1e-6) #Re considéré = 0
    

    
    Hp = dHp_dt(dm_FeO_dt, dm_O2_dt, Tp, Tg, A_FeO, q_FeO, h_O2, W_O2, h_p, sigma, epsilon)
    
    
  
    l = []
    z = 0
    mem = 0
    for i in range(len(Matrice_H)):
     
    
         Hp_new = ((m[1])/M_Fe*Matrice_H[i][1] + (m[0]/M_FeO)*Matrice_H[i][2] + m_SiO2/M_SiO2*Matrice_H[i][5] )*1000
         
        
         
         l.append(Hp_new)
         
         
    z = 10000 # Initialiser avec une valeur infinie
    T_test = Tp
    for i in range(len(l)):
        
        diff = abs(l[i] - m[3])  # Calculer la différence entre H_new et H0

        if diff < z:
           

            z = diff
        
            mem = i
            T_test = Matrice_H[i][0]
            T_guess = Tp
            

    Tp = T_test
 
    Matrice_h.append(l[-1])
    
    test.append(Tp)
    
    

    return [dm_FeO_dt,dm_Fe_dt,dm_O2_dt,Hp]



def euler_explicite(f, y0, t0, tf, h, params,Tg):
    n = int((tf - t0) / h) + 1
    t = np.linspace(t0, tf, n)
    y = np.zeros((n, len(y0)))
    y[0] = y0

    for i in range(n - 1):
        dydt = f(t[i], y[i], params,False,Tg)
        y[i + 1] = y[i] + h * np.array(dydt)
        

    return t, y


# Paramètres d'exemple

params = (nu_Fe_FeO, rho_FeO, k0_FeO, Ta_FeO)
m_0 = [1e-14,m_Fe,0,H0]

t_span = (0, 0.008)
h = 0.0001
t_list = []
y_list =[]

temp = []
    
t, y = euler_explicite(reaction_odes, m_0, t_span[0], t_span[1], h, params,Tg)





def moving_average(data, window_size):
    smoothed = np.convolve(data, np.ones(window_size)/window_size, mode='valid')
    pad_len = (window_size - 1) // 2
    smoothed_padded = np.pad(smoothed, pad_width=pad_len, mode='edge')
    return smoothed_padded

y_smooth = moving_average(test, window_size=15)



# for i in range(len(test)):
#     if test[i] > 1173:
        
#         print(test[i])
        
plt.plot(t[:-1], test, label='T')
plt.show()
#print(test)




def moving_average(data, window_size):
    smoothed = np.convolve(data, np.ones(window_size)/window_size, mode='valid')
    pad_len = (window_size - 1) // 2
    smoothed_padded = np.pad(smoothed, pad_width=pad_len, mode='edge')
    return smoothed_padded



# T_liste = [1210,1215,1219,1220]
# for i in range(len(T_liste)):
#     test = []
#     t, z = euler_explicite(reaction_odes, m_0, t_span[0], t_span[1], h, params,T_liste[i])
    
#     t_list.append(t)
#     y_list.append(z)
#     temp.append(test)
    
    

# plt.plot(t_list[3][:-1], moving_average(temp[3], window_size=11), label='r')
# #plt.plot(t_list[3][:-1], moving_average(temp[0], window_size=11), label='r')
# #plt.plot(t_list[3][:-1], moving_average(temp[1], window_size=11), label='r')
# plt.plot(t_list[3][:-1], moving_average(temp[2], window_size=11), label='r')



# #plt.plot([t_list[3][548], t_list[3][548]], [1100, temp[3][549]], linestyle='--', color='gray')  # ligne verticale
# #plt.plot([t_list[3][0], t_list[3][548]], [temp[3][549], temp[3][549]], linestyle='--', color='gray')  # ligne horizontale

# plt.xlim(0.002,0.008)
# plt.ylim(1175,1300)
# plt.xticks([])
# plt.yticks([])
# ax = plt.gca()
# ax.spines['top'].set_visible(False)
# ax.spines['right'].set_visible(False)
# #plt.xlabel('Time (s)')
# plt.savefig("Termal_runawaySi.pdf")
# #plt.ylabel('Temperature ')
# plt.show()