# -*- coding: utf-8 -*-

import numpy as np
import math
import numpy.linalg as npl
import matplotlib.pyplot as plt
import scipy.stats
import scipy.fft as scifft

def u0(x):
	#return(np.sin(np.pi*x))
	return((x>0.25)*(x<0.75))

kappa = 1. # Coefficient de diffusion
J = 256 # Nombre de point à l'intérieur de ]0, 1[
dx = 1/(J+1)
X = np.linspace(dx,1-dx,J) # Vecteur des positions des degrés de liberté 
T = 0.01 # Temps final
U0 = u0(X)

U = np.copy(U0)
F = scifft.dst(U,type=1)

#plt.plot(X,U,'*-')
#plt.grid()
#plt.show()
#plt.plot(F,'*-')
#plt.grid()
#plt.show()

V = np.exp(-kappa*np.arange(1,J+1)*np.arange(1,J+1)*np.pi*np.pi*T)
F = V*F # Vecteur des amplitudes en fréquence au temps final 
U = scifft.idst(F,type=1) # Solution (approchée !) au temps final

X = np.hstack((0,X,1))
U0 = np.hstack((0,U0,0))
U = np.hstack((0,U,0))
plt.plot(X,U0,label="Donnée initiale")
plt.legend()
leg = "Solution approchée au temps T = " + str(T)
plt.plot(X,U,label=leg)
plt.legend()
plt.grid()
plt.show()

# En fait le phénomène de Gibbs n'est observable que pour des temps très très petits...