(* Mathematica program to evaluate the step potential in Quantum Mechanics.
Are used units where \hbar=1, m=1, so that the Schroedinger equation reads

      -1/2 (d^2 u)/(d x^2) + ( V(x) -E ) u(x) =0

  date : 6/10/2009
  author: Jorge C. Romão
  email : jorge.romao@ist.utl.pt
*)

(* Definition of basic variables *)
k = Function[En, Sqrt[2 En]]
q = Function[{En, V}, If[En>V,Sqrt[2 (En-V)],Sqrt[2 (V - En)]]]
alfa = Function[{En, V}, If[En<V,ArcTan[q[En, V]/k[En]],0]]
w = Function[En, En]
pot=Function[{x,V},If[x<0,0,V]]

(* Definition of R and T *)
R=Function[{En,V},If[En>V,(k[En]-q[En,V])/(k[En]+q[En,V]), 
E^(- 2 I alfa[En,V])]]
T=Function[{En,V},If[En>V,2 k[En]/(k[En]+q[En,V]),
	2 k[En]/Sqrt[k[En]^2+q[En,V]^2] E^(-I alfa[En,V])]]

(* Definition of the wave function *)
psi=Function[{x,t,En,V},If[x<0,(E^(I k[En] x) 
	+ R[En,V] E^(-I k[En] x)) E^(-I w[En] t),
	If[En<V,T[En,V] E^(-q[En, V] x) E^(-I w[En] t), 
	T[En,V] E^(I q[En, V] x) E^(-I w[En] t)]]]

(* Definition of the Real and Imaginary part of the wave function *)
Repsi=Function[{x,t,En,V},ComplexExpand[Re[psi[x,t,En,V]]]]
Impsi=Function[{x,t,En,V},ComplexExpand[Im[psi[x,t,En,V]]]]

(* Plots *)
(* Conventions  Repsi[x, t, En, V], Impsi[x, t, En, V], pot[x,V] *)
Animate[Plot[{Repsi[x, t, 1, 2], Impsi[x, t, 1, 2],pot[x,2]}, 
            {x, -10, 10},PlotRange->{-2,2}],{t,0,100,0.1}]