#*************************************************************************
# programme de simulation de l'effet tunnel sur une barrière de potentiel
# rectangulaire
# Dominique Lefebvre  octobre 2012
# www.tangenteX.com
#*************************************************************************

# initialisation des paramètres de la simulation
me := 9.1e-31:	      # masse de l'électron (USI)
V0 := 10*1.6e-19:     # hauteur de la barrière de potentiel (10 eV en USI)
l := 3.6e-10:         # largeur de la barrière de potentiel (USI)
h := 1.055e-34:	      # hbarre


# définition de l'énergie de l'électron et de sa fréquence
En := 6*1.6e-19:      # énergie de l'électron incident (6 eV)
w := En/h:            # fréquence de l'onde associée à l'électron

# définition des coefficients des exponentielles
k1:=sqrt(2*me*En)/h: 
k2:=sqrt(2*me*(V0-En))/h:

# calcul des exponentielles
psi1 := x -> exp(I*k1*x) + B*exp(-I*k1*x):
psi2 := x -> F*exp(k2*x) + G*exp(-k2*x):
psi3 := x -> C*exp(I*k1*x):

# résolution de l'équation avec les conditions initiales définies
sol := solve({psi1(0) = psi2(0), psi2(l) = psi3(l), D(psi1)(0) = D(psi2)(0), D(psi2)(l)= D(psi3)(l)}, {B, C, F, G}):
assign(sol):

# calcul des fonctions
tpsi1(x,t) := psi1(x)*exp(-I*w*t):
tpsi2(x,t) := psi2(x)*exp(-I*w*t):
tpsi3(x,t) := psi3(x)*exp(-I*w*t):

#définition de la fonction "barrière de potentiel"
barriere := piecewise(0<x and x<l,V0/En, x<0,0, l<x,0):

# définition des courbes graphiques
with(plots):
c1 := animate(Re(tpsi1(x,t)), x=-4*l..0, t=0..2*Pi/w):
c2 := animate(Re(tpsi2(x,t)), x=0..l, t=0..2*Pi/w):
c3 := animate(Re(tpsi3(x,t)), x=l..l+4*l, t=0..2*Pi/w):

# définition de la barrière
cbarriere := plot(barriere, x=-4*l..4*l, color=blue, thickness=2):

# tracé des courbes et de la barrière
display([c1, c2, c3, cbarriere]);

# calcul du coefficient de transmission T
T := C*conjugate(C):
T := evalc(T):
T : simplify(T);