#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# Code fait par Sadek Al-Jibouri avec l'aide des widgets de Martin Verot.
# Le code est libre d'accès, à ne pas utiliser à des fins commerciaux
# Il permet de visualiser l'evolution de l'aimantation stable pour un ferro quelconque quand on s'approche de ou depasse sa temperature de Curie.

import matplotlib.pyplot as plt
import numpy as np
import widgets
import scipy.constants as constants
from matplotlib import rc
from numpy.random import random
from math import floor
from scipy.optimize import root_scalar
from scipy.special import cbrt

Tmax = 1800
N = 100
a0 = 2*1E-7
b0 = 1E-7
G0 = 1
Bmax = 1
T_c = 1500
V = 1e-2

M = np.linspace(-100,100)



parameters = {
    'T' : widgets.FloatSlider(value=T_c, description='Température', min=0, max=2*T_c),
    'B' : widgets.FloatSlider(value=0, description='Champ magnétique', min=-Bmax, max=Bmax)
    }

# def find_cubic_roots(a,b,c,d, bp = False):
#     # with a*x**3 + b*x**2 + c*x + d = 0
#     a,b,c,d = a+0j, b+0j, c+0j, d+0j
#     all_ = (a != np.pi)
# 
#     Q = (3*a*c - b**2)/ (9*a**2)
#     R = (9*a*b*c - 27*a**2*d - 2*b**3) / (54 * a**3)
#     D = Q**3 + R**2
#     S = 0 #NEW CALCULATION FOR S STARTS HERE
#     if np.isreal(R + np.sqrt(D)):
#         S = cbrt(np.real(R + np.sqrt(D)))
#     else:
#         S = (R + np.sqrt(D))**(1/3)
#     T = 0 #NEW CALCULATION FOR T STARTS HERE
#     if np.isreal(R - np.sqrt(D)):
#         T = cbrt(np.real(R - np.sqrt(D)))
#     else:
#         T = (R - np.sqrt(D))**(1/3)
# 
#     result = np.zeros(tuple([3])) + 0j
#     result[all_,0] = - b / (3*a) + (S+T)
#     result[all_,1] = - b / (3*a)  - (S+T) / 2 + 0.5j * np.sqrt(3) * (S - T)
#     result[all_,2] = - b / (3*a)  - (S+T) / 2 -  0.5j * np.sqrt(3) * (S - T)
#     #if bp:
#         #pdb.set_trace()
#     return result
    


def G(T,M,B):
    return(G0 + a0*(T-T_c)*M**2 + 0.5*b0*M**4 - V*B*M)

def Mmin(T,B):
    return(np.roots([2*b0, 0, 2*a0*(T-T_c), -V*B]))
    #return(find_cubic_roots(2*b0, 0, 2*a0*(T-T_c), -V*B))
    
def Dseconde(T,M,B):
    return(a0*(T-T_c)*2 + 3*2*b0*M**2)

Tlist = []
Tlist_instable = []
Mlist = []
Mlist_instable = []

def plot_data(T,B):
    root = []
    lines["G"].set_data(M,G(T,M,B))
    rootbrut = Mmin(T,B)
    for i in rootbrut :
        realroot = np.real(i)
        if np.abs(i - realroot) <=1e-7:
            if Dseconde(T,realroot,B) >= 0:
                Tlist.append(T)
                Mlist.append(realroot)
            else : 
                Tlist_instable.append(T)
                Mlist_instable.append(realroot)
            root.append(realroot)
    Lminima = [G(T,i,B) for i in root]
    lines["Bif"].set_data(Tlist,Mlist)
    lines["Bifinstable"].set_data(Tlist_instable,Mlist_instable)
    
    lines["Minima"].set_data(root,Lminima)
    lines["Bifinstant"].set_data([T]*len(root), root)
    fig.canvas.draw_idle()
    
fig = plt.figure(figsize=(12,6))
fig.suptitle("Illustration de la transition Ferro-Para".format(N))

ax = fig.add_axes([0.07, 0.3, 0.57, 0.6])
ax2 = fig.add_axes([0.72, 0.3, 0.25, 0.6])

lines = {}
lines['G'], = ax.plot([], [], color='red', label = "Potentiel enthalpie libre")
lines['Bif'], = ax2.plot([], [],"o", color='blue', label = "Equilibres stables")
lines['Bifinstable'], = ax2.plot([], [],"o", color='red', label = "Equilibres instables")
lines['Bifinstant'], = ax2.plot([], [],"o", color='red')
lines["Minima"], = ax.plot([], [], "o",color='red', label = "Minima")



xmin, xmax, ymin, ymax = -100,100,G0-1,G0+1



ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

ax2.set_xlim(0,2*T_c)
ax2.set_ylim(-100,100)

ax.set_xlabel('Aimantation M')
ax.set_ylabel('Enthalpie libre G')
def func(self):
    Tlist.clear()
    Mlist.clear()
    Tlist_instable.clear()
    Mlist_instable.clear()



ax2.set_xlabel("Température T (K)")
ax2.set_ylabel('Aimantation M')
ax2.legend(loc = "lower right")
ax.legend(loc = "lower right")
param_widgets = widgets.make_param_widgets(parameters, plot_data, slider_box=[0.17, 0.07, 0.42, 0.15])
#choose_widget = widgets.make_choose_plot(lines, box=[0.015, 0.25, 0.2, 0.15])
reset_button = widgets.make_reset_button(param_widgets)
reset_button.on_clicked(func)

if __name__=='__main__':
    plt.tight_layout()
    plt.show()