(* Mathematica program to evaluate the transmission coefficient 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

 Are shown 4 different potentials and it is possible to add more. We use the 
 transfer matrix method to solve numerically non constant potentials.

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

(* Transfer Matrix. The notation is

 | C |    |MT(1,1)  MT(1,2)|  | A |
 |   | =  |                |  |   |
 | D |    |MT(2.1)  MT(2,2)|  | B |

where x0 is the coordinate of the separation and  K1, K2 are the wave 
numbers for x < x0 and x> x0, respectively, with the convention

          | Sqrt[2 m/\hbar^2 (V(x) - E)     if E < V(x0)
      K = |
          | I Sqrt[2 m/\hbar^2 (E - V(x))   if E > V(x0)

*)

MT=Function[{x0,K1,K2},m11=(K1+K2)/(2 K2) E^(x0 K1 -x0 K2);
	    m12=(-K1+K2)/(2 K2) E^(-x0 K1 -x0 K2);
	    m21=(-K1+K2)/(2 K2) E^(x0 K1 + x0 K2);
	    m22=(K1+K2)/(2 K2) E^(-x0 K1 + x0 K2);{{m11,m12},{m21,m22}}]

(* Transmission Coefficient *)

Tsq=Function[M,Abs[M[[1,1]]-M[[1,2]]*M[[2,1]]/M[[2,2]]]^2]

(* Atomic Units *)

uaE = 27.2 (* eV *)
uaL = 0.053 (* nm *)

(* Barrier V(x)  entre -a/2 < x < a/2 -> Width=a *)

V1=Function[{x,V0,a},If[x<-a/2,0,If[x<a/2,V0, 0]]] 

eps=0.001
V2=Function[{x,V0,a},If[x<-a/2,0,If[x<=-a/2 + a eps,V0/(a eps) (x+a/2),
	   If[x<a/2, V0/a/(1-eps) (a/2-x),0]]]]

V3=Function[{x,V0,a}, If[x<-a/2,0,If[x<0,V0 (2/a) (x+a/2),
            If[x<a/2,V0 (2/a) (a/2-x),If[x>=a/2,0]]]]]
V4=Function[{x,V0,a},If[x<-a/2,0,If[x<a/2,V0 (-4/a^2) (x-a/2)*(x+a/2),0]]]

Nmax = 5;
delta=2./Nmax;
xborders = Table[-1+delta*i, {i, -1, Nmax+1 }];
xmedaux=Table[(xborders[[i]]+delta/2),{i,1,Nmax+2}];
V1apr=Function[{V0,a},Table[V1[a/2*xmedaux[[i]],V0,a],{i,1,Nmax+2}]];
V2apr=Function[{V0,a},Table[V2[a/2*xmedaux[[i]],V0,a],{i,1,Nmax+2}]];
V3apr=Function[{V0,a},Table[V3[a/2*xmedaux[[i]],V0,a],{i,1,Nmax+2}]];
V4apr=Function[{V0,a},Table[V4[a/2*xmedaux[[i]],V0,a],{i,1,Nmax+2}]];
xaux=Table[xborders[[i]],{i,2,Nmax+2}];

T2V1=Function[{En,V0,a}, kappa =Table[If[V1apr[V0,a][[ik]]<En,
	I Sqrt[2 (En-V1apr[V0,a][[ik]])],Sqrt[2 (V1apr[V0,a][[ik]]-En)]],
	{ik,1,Nmax+2}];
    Res = IdentityMatrix[2];
    For[ir = 1, ir <=Nmax+1, 
	Res =  MT[a/2 xaux[[ir]], kappa[[ir]], kappa[[ir + 1]]] . Res; ir++];
    Tsq[Res]]

T2V2=Function[{En,V0,a}, kappa =Table[If[V2apr[V0,a][[ik]]<En,
	I Sqrt[2 (En-V2apr[V0,a][[ik]])],Sqrt[2 (V2apr[V0,a][[ik]]-En)]],
	{ik,1,Nmax+2}];
    Res = IdentityMatrix[2];
    For[ir = 1, ir <=Nmax+1, 
	Res =  MT[a/2 xaux[[ir]], kappa[[ir]], kappa[[ir + 1]]] . Res; ir++];
    Tsq[Res]]

T2V3=Function[{En,V0,a}, kappa =Table[If[V3apr[V0,a][[ik]]<En,
	I Sqrt[2 (En-V3apr[V0,a][[ik]])],Sqrt[2 (V3apr[V0,a][[ik]]-En)]],
	{ik,1,Nmax+2}];
    Res = IdentityMatrix[2];
    For[ir = 1, ir <=Nmax+1, 
	Res =  MT[a/2 xaux[[ir]], kappa[[ir]], kappa[[ir + 1]]] . Res; ir++];
    Tsq[Res]]

T2V4=Function[{En,V0,a}, kappa =Table[If[V4apr[V0,a][[ik]]<En,
	I Sqrt[2 (En-V4apr[V0,a][[ik]])],Sqrt[2 (V4apr[V0,a][[ik]]-En)]],
	{ik,1,Nmax+2}];
    Res = IdentityMatrix[2];
    For[ir = 1, ir <=Nmax+1, 
	Res =  MT[a/2 xaux[[ir]], kappa[[ir]], kappa[[ir + 1]]] . Res; ir++];
    Tsq[Res]]

(* Exact formula for the square barrier *)
T2Barrier=Function[{En,V0,a}, k=Sqrt[2 En]; q=Sqrt[2 (V0-En)];
		   (2*k*q)^2/((2*k*q)^2 + (k^2+q^2)^2*Sinh[ q a]^2)];
 

(* Approximate formula for the square barrier *)
T2BarrierApprox=Function[{En,V0,a}, k=Sqrt[2 En]; q=Sqrt[2 (V0-En)];
		   (4*k*q)^2/(k^2+q^2)^2*E^(-2 q a)];

(* Plot Potentials *)
(* Conventions:
Plot[{V1[x, V0, a], V2[x, V0, a], V3[x, V0, a], V4[x, V0, a]}, {x, xmin, xmax}]
 *)

Plot[{V1[x, 1, 1], V2[x, 1, 1], V3[x, 1, 1], V4[x, 1, 1]}, {x, -0.6,  0.6}]

(* Plot |T|^2 *)
(* Conventions:
Plot[{T2V1[En,V0,a],T2V2[En,V0,a],T2V3[En,V0,a],T2V4[En,V0,a],
	T2Barrier[En,V0,a]},{a,amin,amax}]
 *)

FunctionPlot=Plot[{T2V1[0.5,1,a], T2V2[0.5,1,a], T2V3[0.5,1,a], 
		   T2V4[0.5,1,a]}, {a,0.01,5}]

DataExact=Table[{a,T2Barrier[0.5,1,a]},{a,0.01,5,0.08}]
DataAprox=Table[{a,T2BarrierApprox[0.5,1,a]},{a,0.01,5,0.08}]
PlotDataExact=ListPlot[{DataExact,DataAprox}, PlotStyle->PointSize[0.01]]

Show[{FunctionPlot,PlotDataExact}]