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

# Schéma upwind pour les ondes avec Dirichlet homogène. 

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

def u0(x): # donnée initiale
#    return np.exp(-50*(x-0.5)*(x-0.5))
#    return 2.*abs(x - 0.5) - 1
#    return ((x>0.25)*(x<0.75)).astype(float)
    return np.sin(2.*np.pi*x)

def u1(x): # dérivée en temps initiale
	return 0.*x
#	return (x>(0.5 - 0.5*dx))*(x<=(0.5 + 0.5*dx)).astype(float)/dt

L = 1.
J = 100
T = 3
c = 1 # vitesse des ondes
C = 0.5 # Nombre de Courant

dx = L/(J+1)
X = np.linspace(dx,L-dx,J) # Abscisses des points
dt = C*dx/c

U = u0(X)
V = u0(X) + dt*u1(X)

#V = (u1(X) - u0(X))/dt
#W = (U - np.concatenate(([0],U[:J-1])))/dx
#A = V + c*W
#B = V - c*W

n = 0
t = 0.

plt.grid('on')

while t<T:
	plt.clf()
	dt = min(dt,T-t)
	t = t + dt
	n = n + 1

	W = np.copy(V)	
	V = 2*V - U + c*c*dt*dt/dx/dx*(np.concatenate((V[1:],[0])) - 2*V + np.concatenate(([0],V[:J-1])))
	U = W

	plt.plot(X,U)
	plt.xlabel("x")
	plt.ylabel("U, solution approchée")
	plt.title("Solution approchée au temps t = " + str(t))
	plt.ylim(-1,1)
	plt.grid()
	plt.draw()
	plt.pause(dt)

	print ("itération  ", n)
	print("t = ", t)

plt.plot(X,U)
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()

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