(* Mathematica program for Problem 4.10

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

(* Ficheiro do Mathematica para a solucao do problema 6 *)

(* Input de variaveis *)
Clear[lbd,eta]

lbd=20;
eta=2;
c=-4.;
d=4.;

(* 
 Valor para comecar a encontrar a raiz com FindRoot. E preciso confirmar 
 com os  graficos para ver se faz sentido. O Mathematica encontra as 
 solucoes para y > yini
*)

yini=Pi/2 + 0.001   ;

(* Fim do input *)

(* Definicao das funcoes para a solucoes par, impar, eta >>1, e Delta y *)

f = Function[{y,lbd}, 2 y Sqrt[lbd-y^2]/(2 y^2 -lbd)];

feven = Function[{y, lbd, eta}, y Sqrt[lbd - y^2]*(1 + 
            Tanh[eta Sqrt[lbd - y^2]])/(y^2 - (lbd - y^2)*
              Tanh[eta*Sqrt[lbd - y^2]])];

fodd = Function[{y, lbd, eta}, y Sqrt[lbd - y^2]*(1 + 
            1/Tanh[eta Sqrt[lbd - y^2]])/(y^2 - (lbd - y^2)/
              Tanh[eta*Sqrt[lbd - y^2]])];

Deltay = Function[{y, lbd, eta}, 4 E^(-2*eta*Sqrt[lbd - y^2])*
      y*(lbd - y^2)/lbd/(2 + Sqrt[lbd - y^2])];

(* Calculo das solucoes: y0, yeven, yodd *)

r0 = FindRoot[f[y, lbd] == Tan[y], {y, yini}]; y0 = y /. r0 ;

reven = FindRoot[feven[y, lbd, eta] == Tan[y], {y, yini}]; yeven = y /. reven ;

rodd = FindRoot[fodd[y, lbd, eta] == Tan[y], {y, yini}]; yodd = y /. rodd ;

Deltayexact= yodd-yeven;

Print["yeven=",yeven,", yodd=",yodd,", Deltay=",Deltayexact] 

(* Grafico para observar as solucoes *)

ploteven=Plot[{Tan[y],feven[y,lbd,eta]},{y,0,Sqrt[lbd]},
PlotRange->{c,d}, AxesLabel->{"y","Caso Par"}]
plotodd=Plot[{Tan[y],fodd[y,lbd,eta]},{y,0,Sqrt[lbd]},
PlotRange->{c,d}, AxesLabel->{"y","Caso Impar"}]

(* Calculo das funcoes de onda par e impar *)

alfa=Function[{y,lbd},Sqrt[lbd-y^2]];

C4=1;
C1even=C4;
C1odd=-C4;
A2even=C4 E^(-(eta+1)*alfa[yeven,lbd])*(Cos[(eta+1)*yeven] +
		alfa[yeven,lbd]/yeven*Sin[(eta+1)*yeven]);
B2even=C4 E^(-(eta+1)*alfa[yeven,lbd])*(Sin[(eta+1)*yeven] -
		alfa[yeven,lbd]/yeven*Cos[(eta+1)*yeven]);
A1even=A2even;
B1even=-B2even;
C2even=(A2even*Cos[yeven*eta]+B2even*Sin[yeven*eta])/Cosh[eta*alfa[yeven,lbd]];
A2odd=C4 E^(-(eta+1)*alfa[yodd,lbd])*(Cos[(eta+1)*yodd] +
		alfa[yodd,lbd]/yodd*Sin[(eta+1)*yodd]);
B2odd=C4 E^(-(eta+1)*alfa[yodd,lbd])*(Sin[(eta+1)*yodd] -
		alfa[yodd,lbd]/yodd*Cos[(eta+1)*yodd]);
A1odd=-A2odd;
B1odd=B2odd;
C3odd=(A2odd*Cos[yodd*eta]+B2odd*Sin[yodd*eta])/Sinh[eta*alfa[yodd,lbd]];

(* Funcoes nao normalizadas *)

psievenaux=Function[{x,lbd,eta},
If[x<-(eta+1),C1even E^(alfa[yeven,lbd] x),
   If[x<-eta,A1even*Cos[yeven x]+B1even*Sin[yeven x],
      If[x<eta,C2even*Cosh[alfa[yeven,lbd] x],
	 If[x<eta+1,A2even*Cos[yeven x]+B2even*Sin[yeven x],
	    C4 E^(-alfa[yeven,lbd] x)]]]]];

psioddaux=Function[{x,lbd,eta},
If[x<-(eta+1),C1odd E^(alfa[yodd,lbd] x),
   If[x<-eta,A1odd*Cos[yodd x]+B1odd*Sin[yodd x],
      If[x<eta,C3odd*Sinh[alfa[yodd,lbd] x],
	 If[x<eta+1,A2odd*Cos[yodd x]+B2odd*Sin[yodd x],
        C4 E^(-alfa[yodd,lbd] x)]]]]];

(* Normalizacao *)

Neven=Integrate[psievenaux[x,lbd,eta]^2,{x,-10*eta,10*eta}];
Nodd=Integrate[psioddaux[x,lbd,eta]^2,{x,-10*eta,10*eta}];

(* Funcoes Normalizadas *)

psieven=Function[x,psievenaux[x,lbd,eta]/Sqrt[Neven]];
psiodd=Function[x,psioddaux[x,lbd,eta]/Sqrt[Nodd]];

Plot[psieven[x],{x,-3-eta,3+eta},AxesLabel->{"x","ψ"}, PlotRange->All]
Plot[psiodd[x],{x,-3-eta,3+eta},AxesLabel->{"x","ψ"}, PlotRange->All]