(*  Gasiorowicz Exemplo 3.5

Este ficheiro destina-se a fazer o Exemplo 3.5 do Gasiorowicz utilizando 
o mathematica. Os comentarios devem ser auto-explicativos.

  date : 26/9/2009
  author: Jorge C. Romão
  email : jorge.romao@ist.utl.pt

*)

(* Definicoes das Funcoes da caixa *)

u=Function[{x,n,a},Sqrt[2/a] Sin[n Pi x/a] ]

(* Definicao da Funcao de Onda *)

psi=Function[{x,a},If[x<0||x>a,0,If[x<a/2,Sqrt[12/a] x/a,Sqrt[12/a] (1-x/a)]]]

(* Calculo dos Coeficientes Am*)

An = Function[n, Integrate[ u[x, n, a] psi[x,a], {x, 0, a}, 
			    Assumptions -> {n ∈ Integers,a>0}]]

Anexplicita=Function[n,Sqrt[96] Sin[n Pi/2]/n^2/Pi^2]

(* Funcao de Onda Aproximada. nmax e o maximo numero de termos *)

psiapprox=Function[{x,a,nmax}, Sum[Anexplicita[n] u[x,n,a],{n,1,nmax}]]

(* Testes *)
(* Expressao para os coeficientes *)

Teste1:= An[n] -Anexplicita[n]      

(* Normalizacao *)

Teste2:=Integrate[psi[x,a]^2,{x,0,a},Assumptions -> a>0]

(* Plots a=1 *)
(* So um termo *)

plot1=Plot[psiapprox[x,1,1],{x,0,1},PlotRange->{{0,1},{0,1.8}},
AxesLabel -> {"x/a","ψ"}];

(* 10 termos *)
     plot10=Plot[psiapprox[x,1,10],{x,0,1},PlotRange->{{0,1},{0,1.8}},
AxesLabel -> {"x/a","ψ"}] ;

(* 100 termos *)
     plot100=Plot[psiapprox[x,1,100],{x,0,1},PlotRange->{{0,1},{0,1.8}},
AxesLabel -> {"x/a","ψ"}];

Fig1=Show[{plot1,plot10,plot100}]

(* Valor Medio <x> *)

xmed=Integrate[psi[x,a] x psi[x,a],{x,0,a},Assumptions-> a>0]

(* Valor Medio <x^2> *)

x2med=Integrate[psi[x,a] x^2 psi[x,a],{x,0,a},Assumptions-> a>0]

(* Delta x *)

Deltax=Simplify[Sqrt[x2med -xmed^2],Assumptions-> a>0]

(* Funcao de Onda no espaco dos momentos *)

Cst=1/Sqrt[2 Pi hbar] 

phi=Function[{p,a},Cst Integrate[ psi[x,a] E^(-I/hbar p x),{x,0,a}, 
Assumptions-> a>0]]

(* Plot Distribuicao no espaco dos momentos *)

phipa = phi[p, a];

Fig2=Plot[Simplify[Abs[phipa /. {p -> hbar/a k}]^2  hbar/a,
  Assumptions -> {a > 0, hbar > 0}] , {k, -10, 10}, 
AxesLabel -> {"p a/ℏ", "|ϕ|^2"}]

(* Valor Medio <p> *)

pmed=Integrate[Abs[phipa]^2 p, {p, -Infinity, Infinity}, 
 Assumptions -> {a > 0, hbar > 0}]

(* Valor Medio <p^2> *)

p2med=Integrate[Abs[phipa]^2 p^2, {p, -Infinity, Infinity}, 
 Assumptions -> {a > 0, hbar > 0}]

(* Delta p *)

Deltap=Simplify[Sqrt[p2med-pmed^2],Assumptions->{a>0,hbar>0}]

(* Relacao de Heisenberg *)
(* Verifique que RIH > hbar/2 *)

RIH=Deltax Deltap