(* Mathematica program for Problem 12.7

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

(* Definitions of constants *)
Clear[alfa,gama];
alfa=1./137.036 ;
gama=(alfa/8)^2*13.6 ;
cst=10^5

(* Exact Eigenvalues *)
E1=Function[z,-gama +2 z gama];
E2=Function[z,-gama -2 z gama];
E3=Function[z,-3 gama + z gama/2 + Sqrt[4 gama^2 + 2/3*gama^2*z 
					+ gama^2*z^2/4]];
E4=Function[z,-3 gama + z gama/2 - Sqrt[4 gama^2 + 2/3*gama^2*z 
					+ gama^2*z^2/4]];
E5=Function[z,-3 gama - z gama/2 + Sqrt[4 gama^2 - 2/3*gama^2*z 
					+ gama^2*z^2/4]];
E6=Function[z,-3 gama - z gama/2 - Sqrt[4 gama^2 - 2/3*gama^2*z 
					+ gama^2*z^2/4]];

(* Plots *)
Fig1=Plot[{E1[z] cst, E2[z] cst, E3[z] cst, E4[z] cst, E5[z] cst, 
	   E6[z] cst},{z,0,10},PlotRange->{-30,30},AspectRatio->0.8,
           AxesLabel->{"β/γ","10⁵ E (eV)"}];
Show[Fig1]

(* Campo Fraco *)
E32BFraco1=Function[z,-gama + z*gama*4/3*3/2 ];
E32BFraco2=Function[z,-gama + z*gama*4/3*1/2 ];
E32BFraco3=Function[z,-gama - z*gama*4/3*1/2 ];
E32BFraco4=Function[z,-gama - z*gama*4/3*3/2 ];
E12BFraco1=Function[z,-5*gama + z*gama*2/3*1/2 ];
E12BFraco2=Function[z,-5*gama - z*gama*2/3*1/2 ];

TableE32BFraco1=Table[{z,E32BFraco1[z] cst},{z, 0, 1, 0.1}];
TableE32BFraco2=Table[{z,E32BFraco2[z] cst},{z, 0, 1, 0.1}];
TableE32BFraco3=Table[{z,E32BFraco3[z] cst},{z, 0, 1, 0.1}];
TableE32BFraco4=Table[{z,E32BFraco4[z] cst},{z, 0, 1, 0.1}];
TableE12BFraco1=Table[{z,E12BFraco1[z] cst},{z, 0, 1, 0.1}];
TableE12BFraco2=Table[{z,E12BFraco2[z] cst},{z, 0, 1, 0.1}];

PlotFraco=ListPlot[{TableE32BFraco1,TableE32BFraco2,TableE32BFraco3,
		    TableE32BFraco4,TableE12BFraco1,TableE12BFraco2},
                    PlotStyle->PointSize->0.011];

Fig2=Plot[{E1[z] cst, E2[z] cst, E3[z] cst, E4[z] cst, E5[z] cst, 
	   E6[z] cst},{z,0,1},PlotRange->{-8,2},AspectRatio->0.8,
           AxesLabel->{"β/γ","10⁵ E (eV)"}];

Show[{Fig1,PlotFraco}]
Show[{Fig2,PlotFraco}]

(* Campo Forte *)
EBForte=Function[{z,ml,ms},gama/3*(-7 + 8*ml*ms) + z*gama*(ml+2*ms)];

TableEBForte1=Table[{z, EBForte[z,1,1/2] cst   },{z, 1, 10, 0.1}];
TableEBForte2=Table[{z, EBForte[z,-1,-1/2] cst },{z, 1, 10, 0.1}];
TableEBForte3=Table[{z, EBForte[z,0,1/2]  cst  },{z, 1, 10, 0.1}];
TableEBForte4=Table[{z, EBForte[z,1,-1/2]  cst },{z, 1, 10, 0.1}];
TableEBForte5=Table[{z, EBForte[z,-1,1/2]  cst },{z, 1, 10, 0.1}];
TableEBForte6=Table[{z, EBForte[z,0,-1/2]  cst },{z, 1, 10, 0.1}];

PlotForte=ListPlot[{TableEBForte1,TableEBForte2,TableEBForte3,
		    TableEBForte4,TableEBForte5,TableEBForte6},
                    PlotStyle->PointSize->0.011];

Fig3=Plot[{E1[z] cst, E2[z] cst, E3[z] cst, E4[z] cst, E5[z] cst, 
	   E6[z] cst},{z,0,5},PlotRange->{-10,6},AspectRatio->0.8,
           AxesLabel->{"β/γ","10⁵ E (eV)"}];

Show[Fig1,PlotFraco,PlotForte]
Show[Fig3,PlotFraco,PlotForte]