# -*- coding: utf-8 -*-
"""
Created on Mon Feb 02 16:36:44 2015

Listing 13.6 : équation de transport par différences finies

@author: Ribot/Grivet
"""
from pylab import *

 # Vitesse
c = 1
  # Discrétisation en espace
xmin = 0.0; xmax = 4.0; npt = 81; nx = npt-2 ; h = (xmax-xmin)/(npt -1)
xx = linspace(xmin,xmax,npt); xx = xx.transpose()
xxint = xx[1:nx+1]
  # Matrices de discrétisation en espace 1D
K1Dpos = (eye(nx , nx)-diag(ones(nx-1),-1))/h;
K1Dneg = (diag(ones(nx-1),1) - eye(nx ,nx))/h;
K1Dcentre = (diag(ones(nx-1),1) - diag(ones(nx-1),-1))/h/2;
  # Discrétisation en temps
Tfin = 1; dt = 0.1*h/c ; ntps = int(Tfin/dt)
# ntps = 3, pour un ou ce
u = hstack((ones_like(xxint[xxint <=1]), zeros_like(xxint[xxint>1])))
meth = input('choix de la methode (chaine up,un ou ce): ')
for t in range(ntps):
    # Upwind- vitesse positive
    if meth == 'up':
        u -= c*dt*dot(K1Dpos,u)
    # Upwind- vitesse négative
    elif meth == 'un':
        u -= c*dt*dot(K1Dneg,u)
    #  Schéma centré
    elif meth == 'ce':
        u -= c*dt*dot(K1Dcentre,u)
uu = hstack([0,u,0]) 
plot(xx,uu);