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

import numpy as np
import numpy.linalg as npl
import matplotlib.pyplot as plt

def u0(x):
#    return np.exp(-50*(x-0.5)*(x-0.5))
#    return 0.*x
#    return (x>0.)
#    return ((x>1.)*(x<3)).astype(float)
#    return np.sin(2.*np.pi*x)
#	return(np.exp(1/(x-1)/(x-3))*(x>1).astype(float)*(x<3).astype(float)) # Pas super, avec beaucoup de mailles (4000 par exemple) ça plante : j'ai l'impression que là vectorize nous sauve... 
	if x > 1:
		if x < 3: 
			return(np.exp(1/(x-1)/(x-3)))
		else:
			return(0.)
	else : 
		return(0.)

L = 10.
J0 = 10
T = 3
C = 0.9 # Nombre de Courant
a = 1. 
K = 11

eup = np.zeros(K)
elf = np.zeros(K)
eru = np.zeros(K)
elw = np.zeros(K)
ebm = np.zeros(K)

for k in range(K):
	J = J0*2**k
	dx = L/(J+1)
	X = np.linspace(dx,L-dx,J) # Abscisses des points du maillage 
	Uup = np.vectorize(u0)(X)
	Ulf = np.vectorize(u0)(X)
	Uru = np.vectorize(u0)(X)
	Ulw = np.vectorize(u0)(X)
	Ubm = np.vectorize(u0)(X)
	n = 0
	t = 0.
	print("Nombre de points : ",J)
	plt.grid('on')
	while t<T:
		plt.clf()
		c = a # Coefficient le schéma de Rusanov : doit être supérieur à a
		dt = C*dx/c
		nu = dt/dx
		dt = min(dt,T-t)
		t = t + dt
		n = n + 1
		Uup = Uup - a*dt/dx*(Uup - np.concatenate(([0],Uup[:J-1]))) # Upwind
		Ulf = 0.5*(np.concatenate(([0],Ulf[:J-1])) + np.concatenate((Ulf[1:],[0]))) - 0.5*a*dt/dx*(np.concatenate((Ulf[1:],[0])) - np.concatenate(([0],Ulf[:J-1]))) # Lax-Friedrichs
		Uru = Uru - 0.5*a*dt/dx*(np.concatenate((Uru[1:],[0])) - np.concatenate(([0],Uup[0:J-1])))  + 0.5*c*dt/dx*(np.concatenate((Uru[1:],[0])) - 2.*Uru + np.concatenate(([0],Uup[0:J-1])))# Rusanov
		Ulw = Ulw - 0.5*a*dt/dx*(np.concatenate((Ulw[1:],[0])) - np.concatenate(([0],Ulw[0:J-1])))  + 0.5*a*a*dt*dt/dx/dx*(np.concatenate((Ulw[1:],[0])) - 2.*Ulw + np.concatenate(([0],Ulw[0:J-1])))# Lax-Wendroff
		Ubm = Ubm - a*dt/dx*(0.5*(3 - a*dt/dx)*Ubm - (2 - a*dt/dx)*np.concatenate(([0],Ubm[0:J-1])) + 0.5*(1 - a*dt/dx)*np.concatenate(([0],[0],Ubm[0:J-2]))) # Beam-Warming

	eup[k] = npl.norm(np.vectorize(u0)(X - a*T) - Uup,np.inf)
	elf[k] = npl.norm(np.vectorize(u0)(X - a*T) - Ulf,np.inf)
	eru[k] = npl.norm(np.vectorize(u0)(X - a*T) - Uru,np.inf)
	elw[k] = npl.norm(np.vectorize(u0)(X - a*T) - Ulw,np.inf)
	ebm[k] = npl.norm(np.vectorize(u0)(X - a*T) - Ubm,np.inf)
	
	plt.plot(X,Uup,label = "Upwind")
	plt.legend()
	plt.plot(X,Ulf,label="Lax-Friedrichs")
	plt.legend()
	plt.plot(X,Uru,label="Rusanov")
	plt.legend()
	plt.plot(X,Ulw,label="Lax-Wendroff")
	plt.legend()
	plt.plot(X,Ubm,label="Beam-Warming")
	plt.legend()
	plt.xlabel("x")
	plt.ylabel("U, solution approchée")
	plt.title("Solutions approchées au temps t = " + str(t))
	plt.grid()
	plt.show()
	
plt.plot(np.log(eup),marker = 'o')
plt.plot(np.log(elf),marker = 'x')
plt.plot(np.log(eru),marker = '*')
plt.plot(np.log(elw),marker = '+')
plt.plot(np.log(ebm),marker = 'o')
plt.grid()
plt.show()

ordreup = (np.log(eup[:-1]) - np.log(eup[1:]))/np.log(2.)
ordrelf = (np.log(elf[:-1]) - np.log(elf[1:]))/np.log(2.)
ordreru = (np.log(eru[:-1]) - np.log(eru[1:]))/np.log(2.)
ordrelw = (np.log(elw[:-1]) - np.log(elw[1:]))/np.log(2.)
ordrebm = (np.log(ebm[:-1]) - np.log(ebm[1:]))/np.log(2.)
print('Ordre du schéma upwind : ')
print(ordreup)
print('Ordre du schéma de Lax-Friedrichs : ')
print(ordrelf)
print('Ordre du schéma de Rusanov : ')
print(ordreru)
print('Ordre du schéma de Lax-Wendroff : ')
print(ordrelw)
print('Ordre du schéma de Beam-Warming : ')
print(ordrebm)

#############################