(* ::Package:: *)

(* *** 
Rechnernutzung in der Physik, WS23/24

Teilchen trifft auf Potentialbarrierie/Potentialsenke
*** *)

(* Definiere Potential: *)

V[x_,a_] := If[ Abs[x]<a, V0, 0];

(* Fragestellungen
1. Loesung fuer Wellenfunktion in den 3 Bereichen
2. Verwende Anschlussbedingungen, um die freien Konstanten 
   zu bestimmen.
3. Berechne Transmissions- und Reflexionskoeffizient
4. Betrachte die Energieabhaengigkeit des Transmissions- 
   und Reflexionskoeffizienten (sowohl fuer V0>0 und V0<0)
5. Ortsabhaengigkeit der Wellenfunktion und der 
   Wahrscheinlichkeitsdichte
6. Darstellung der Zeitabhaengigkeit des 
   Real- und Imaginaerteils der Wellenfunktion *)


(* zu 1. *)
(* ((hbar/i d/dx)^2/2/m + V) psi =  En psi
   Beachte: 
   Einlaufende Welle ist auf "1" normiert
   Setze hbar = 1 *)
psiL[x_] =     Exp[I*k*x]  + R * Exp[-I*k*x];
psiM[x_] = A * Exp[I*kM*x] + B * Exp[-I*kM*x];
psiR[x_] = T * Exp[I*k*x];

kRule = {kM -> Sqrt[2 * m * (En-V0)], k -> Sqrt[2 * En * m]};


(* zu 2. *)

eq = {psiL[-a]  == psiM[-a],
      psiL'[-a] == psiM'[-a],
      psiM[a]   == psiR[a],
      psiM'[a]  == psiR'[a]}

ABRTrule = Solve[ eq, {A,B,R,T} ][[1]]
ABRTrule = ABRTrule//FullSimplify


(* zu 3.: Transmissions- und Reflexionskoeffizient *)

(* 
  Transmissionskoeffizient: T = |j_trans/j_ein|
  Reflexionskoeffizient:    R = |j_refl/j_ein|
  Der Strom berechnet sich mittels
    j = 1 / (2*I*m) * ( psi^* d psi/dx - psi d psi^*/dx) 

  Definiere zunaechst die komplex konjugierte Wellenfunktion: *)

psiCrule = {Complex[a_,b_] :> Complex[a,-b],
            R -> RC,
            T -> TC,
            A -> AC,
            B -> BC}

ABRTruleC = {RC -> R,
             TC -> T,
             AC -> A,
             BC -> B} /. ABRTrule /. psiCrule


strom[psi_] := 1 / (2*I*m) * 
    ( (psi /. psiCrule) * D[psi,x] - D[(psi/.psiCrule),x] * psi )

Einlstrom = strom[ psiL[x] ] /. {R->0,RC->0}
(Rstrom   = Einlstrom - strom[ psiL[x] ] /. ABRTrule /. ABRTruleC 
            // ExpToTrig // Simplify)
(Tstrom   = strom[ psiR[x] ] /. ABRTrule /. ABRTruleC 
            // ExpToTrig // Simplify)

RR = Rstrom / Einlstrom
TT = Tstrom / Einlstrom

(* Check: RR+TT = 1: *)

Print["R + T = ", RR + TT // Simplify];


(* zu 4.: Energieabhaengigkeit des Transmissions- 
   und Reflexionskoeffizienten  *)

num = {m->1, a->1};

RRTT[En_, V0_] := Evaluate[ {RR, TT} /. kRule /. num ];

(* Test zur num. Auswertung *)
Print["{R, T}[En=8.5, V0=-10.2] = ",RRTT[8.5, -10.2]//N ];


(* Betrachte Potentialtopf der Tiefe V0 = -100 *)

pl1 = Plot[ Evaluate[ RRTT[en,-100]], {en,1,300}, 
	    PlotLabel -> "V0 = -100",
	    PlotStyle->{Red,Blue}];

(* Betrachte Potentialbarriere der Hoehe V0 = 2 *)
pl2 = Plot[ Evaluate[ RRTT[en,2]], {en,0,15},
	    PlotLabel -> "V0 = +2",
	    PlotStyle->{Red,Blue}];

GraphicsGrid[ { {pl1, pl2} } ]


(* zu 5.: Ortsabhaengigkeit der Wellenfunktion und der 
   Wahrscheinlichkeitsdichte *)

{psiL[x_,aa_], psiM[x_,aa_], psiR[x_,aa_]} = 
{psiL[x], psiM[x], psiR[x]} /. ABRTrule /. kRule /. {m->1,a->aa}

psi[x_,a_] := Which[      x <= -a, psiL[x,a], 
		     -a < x <= a,  psiM[x,a],
        	      a < x,       psiR[x,a]];

Print["psi[x=-1.5, a=1] = ", psi[-1.5,1]];
Print["psi[x=-1.3, a=2, En=1, V0=1.1] = ", Re[ psi[-1.3,2] /. {En->1,V0->1.1} ] ];


pl1 = Plot[ Re[ psi[x,2] /. {En->1,V0->-12} ], {x,-15,15},
	    Epilog -> {Thickness[0.03], Line[{{-2,0},{2,0}}] }];

pl2 = Plot[ psi[x,2] * Conjugate[psi[x,2]] /. {En->1,V0->-12}, {x,-15,15},
	    Epilog -> {Thickness[0.03],Line[{{-2,0},{2,0}}] }];

GraphicsGrid[ { {pl1, pl2} } ]


(* zu 6.: Zeitabhaengigkeit des 
   Real- und Imaginaerteils der Wellenfunktion *)

(* Um die Zeitentwicklung zu erhalten, muss die Wellenfunktion mit
   Exp[-I*En*t] multipliziert werden.
   Betrachte Real- und Imaginaerteil separat. *)

aa = 2.0;

psinum[t_,x_] := {Re[ Exp[-I*En*t]*psi[x,aa] ],
		  Im[ Exp[-I*En*t]*psi[x,aa] ]} /. {En ->5.0 (*2.0*), V0 -> 10.(*-12.0*)};

Table[ psinum[t,10], {t,0,10,1} ]

plotRe[t_] := Plot[ psinum[t,x][[1]], {x,-15,15},
		    PlotStyle -> {Blue},
		    PlotRange -> {-2,2},
		    Epilog -> {Thickness[0.03],Line[{{-aa,0},{aa,0}}]},
		    PlotLabel -> "t = " <> ToString[t]];

plotIm[t_] := Plot[ psinum[t,x][[2]], {x,-15,15},
		    PlotStyle -> {Red},
		    PlotRange -> {-2,2},
		    Epilog -> {Thickness[0.03],Line[{{-aa,0},{aa,0}}]}];

Animate[ {plotRe[t],plotIm[t]}, {t,0,10}, AnimationRate->0.5 ]