#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Jul 18 10:48:03 2022

@author: adomps
"""


##  VOUS NE DEVEZ MODIFIER QUE LES LIGNES PORTANT
## LE COMMENTAIRE << LIGNE_ELEVE >>.

##  Calculs en lien avec l'exercice << Explosion thermique >>
## tiré de Landau (méca flu).
# gamma est noté lambda par Landau.
# Selon la valeur de gamma, theta_1 existe ou non

import numpy as np
import matplotlib.pyplot as plt

## Valeurs numériques
## Pour déterminer l'existence d'un régime stationnaire,
## il suffit de se donner gamma.
## Les valeurs des paramètres sont données pour le
## cas où on résout le régime variable numériquement.


p0 = 240000  ## LIGNE ÉLÈVE ; valeur à modifier

alpha = 0.01
L = 1
ell = L / 2 #  ATTENTION : demi-longueur du réacteur
Lambda = 400 # j'ai pris les valeus du cuivre. Hum hum.
D = 1.2e-4
gamma = p0 * alpha * ell**2 / Lambda  # valeur critique gamma = 0.88
print('gamma', gamma)
T0 = 10



 

####################################################################
##      Régime stationnaire
####################################################################

# On trace une fonction f et on
# résout f = Cst. On n'a pas de solution si la Cste
# dépasse le max de f, égal à 0,6627, d'où
# le paramètre critique 2 * 0,6627**2 = 0,878

theta1 = np.linspace(0, 5, 100000)
f = np......#  LIGNE ÉLÈVE : calculer le tableau numpy des valeurs de la fonction f
cible = np.sqrt(gamma/2)


plt.figure(0)
plt.clf()
plt.title(" Recherche d'existence d'un régime stationnaire ")
plt.xlabel(r'$\theta_1$')
plt.plot(theta1, f, label = r'$f(\theta_1)$', color = 'black')
plt.hlines(cible, 0, np.max(theta1), linestyle = 'dashed')
plt.text(-0.5, cible, r'$\sqrt{\gamma/2}$')
plt.grid()
plt.legend(loc = 'lower right')


# Recherche de theta_1 (quand il exsite).
# De theta_1, on peut déduire la température
# T_1 au centre

def resolution_theta1(gamma) :
    theta1 = np.linspace(0, 2, 100000)
    cible = np.sqrt(gamma/2)
    f = np.exp(-theta1/2) * np.arccosh(np.exp(theta1/2))    
    indice_soluce = np.argmin(abs(f - cible))
    return(theta1[indice_soluce])    
    

   
if gamma > 0.88 : 
    print("Il n'existe pas de régime stationnaire.")
else :
    theta1 = resolution_theta1(gamma) 
    print('theta1 solution', resolution_theta1(gamma))    
    T_max_attendu = T0 + theta1 / alpha
    print('Régime staionnaire avec T_max attendu ', T_max_attendu)
 


######################################################################
#   Régime variable.
#######################################################################
#
#  Résolution de l'éq de la chaleur 1d par méthode explicite
#
#  Le terme source, les conditions de bords et tous les paramètres
#  du problème sont traités comme des variables globales.
#  Les conditions de bord sont exprimés dans une fonction dédiée.
#
#  Une fonction d'animation est importée depuis un fichier séparé.
# 
#
##############################################################################

# Paramètres du calcul numérique
N = 100
x = np.linspace(0, L, N)
delta_x = x[1] - x[0]
delta_t = 0.1
prefac = 0.4 # en dimension 1, la condition de stabilité est prefac < 1/2
delta_t = prefac * delta_x**2 / D
q = D * delta_t / Lambda


def pdt_vecto(T_anc) :
    T_nouv = np.zeros(N)
    T_nouv[1:-1] = T_anc[1:-1] + prefac * (T_anc[0:-2] + T_anc[2:] - 2 * T_anc[1:-1]) \
                    + q * p_source[1:-1]
    return T_nouv
                    

def bords_fixes(T_anc, T_nouv) :
    T_nouv[0] = T_anc[0]
    T_nouv[-1] = T_anc[-1]
    

def evolution(T_init, N_pas, fonction_pdt, fonction_bords) :
    champ_T = np.zeros((N_pas+1, N))
    champ_T[0] = T_init
    for i in range(N_pas) :
        champ_T[i+1] = fonction_pdt(champ_T[i])
        fonction_bords(champ_T[i], champ_T[i+1] )
    return champ_T

def representation(x, champ_T, modulo, num_fig) :
    ''' Représente le champ de température en fonction de x, 
    à diverses dates multiplies de modulo. '''
    plt.figure(num_fig)
    plt.xlabel('$x$ (m)')
    plt.ylabel('T (°C)')
    N_pas = len(champ_T)
    for j in range(0, N_pas, modulo) :
        plt.plot(x, champ_T[j])
    return


def T_mode(x, t, n, A ):
    ''' Solution en mode propre. Pour l'utiliser comme condition initiale,
    prendre t=0. On peut bien sûr faire la somme de plusieurs modes.'''
    return  A * np.sin(n * np.pi * x / L) * np.exp(- D * t * n** 2 * np.pi**2 / L**2)


def source_reaction_chi(T, p0, alpha, T0) :
    return  p0 * np.exp(alpha * (T - T0))


def evolution_reacteur(T_init) :
    global p_source
    champ_temp = np.zeros((N_pas+1, N))
    champ_temp[0] = T_init
    for i in range(N_pas) :
        p_source = source_reaction_chi(champ_temp[i], p0, alpha, T0)
        champ_temp[i+1] = pdt_vecto(champ_temp[i])
        bords_fixes(champ_temp[i], champ_temp[i+1])
    return champ_temp




T_init = T0 + T_mode(x, 0, 1, 0.6) + T_mode(x, 0, 3, -0.2)
N_pas = 30000

champ_T = evolution_reacteur(T_init)


if gamma < 0.88 :
    fin = N_pas   # on garde toute l'animation 
else : 
    fin = 3000 # on garde seulement le début, à la fin on a des NaN


representation(x, champ_T[0:fin], 100, 1)


from fonction_animation import studio_poinca
#  Décommenter et exécuter cette ligne séparément.    
# studio_poinca(x, champ_T[0:fin], delta_t)