(* Mathematica program for Problem 12.5

  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

A=N[5.8*10^(-6),20];
a=N[1.16*10^(-4),20];
b=N[-1.75*10^(-7),20];

(* Exact Eigenvalues *)
E4=Function[B,A/4+(a/2+b/2)*B];
E2=Function[B,A/4-(a/2+b/2)*B];
E3=Function[B,-A/4 + Sqrt[A^2/4+(a/2-b/2)^2*B^2]];
E1=Function[B,-A/4 - Sqrt[A^2/4+(a/2-b/2)^2*B^2]];

(* Plots *)
Fig1=Plot[{E1[B] cst, E2[B] cst, E3[B] cst, E4[B] cst}, {B, 0, 2}, 
 PlotRange -> {-15, 15}, AspectRatio -> 0.8, 
 AxesLabel -> {"B (T)", "10⁵ E (eV)"}]

Fig2=Plot[{E1[B] cst, E2[B] cst, E3[B] cst, E4[B] cst}, {B, 0, 0.1}, 
 PlotRange -> {-1, 1}, AspectRatio -> 0.8, 
 AxesLabel -> {"B (T)", "10⁵ E (eV)"}]

(* Campo Fraco B < 0.01 T*)

E1Fraco=Function[B,-3/4 A ];
E2Fraco=Function[B,1/4 A -1/2 a B];
E3Fraco=Function[B,1/4 A  ];
E4Fraco=Function[B,1/4 A + 1/2 a B ];

BFraco1=Table[{B,E1Fraco[B] cst},{B,0,0.05,0.005}];
BFraco2=Table[{B,E2Fraco[B] cst},{B,0,0.05,0.005}];
BFraco3=Table[{B,E3Fraco[B] cst},{B,0,0.05,0.005}];
BFraco4=Table[{B,E4Fraco[B] cst},{B,0,0.05,0.005}];

FigBFraco=ListPlot[{BFraco1,BFraco2,BFraco3,BFraco4},
		   PlotStyle->PointSize->0.02];
Show[{Fig2,FigBFraco}]

(* Campo Forte  B > 0.1 T *)

E2Forte=Function[B,1/4 A - 1/2 a B -1/2 b B];
E1Forte=Function[B,-1/4 A -1/2 a B + 1/2 b B];
E3Forte=Function[B,-1/4 A +1/2 a B -1/2 b B];
E4Forte=Function[B,1/4 A + 1/2 a B + 1/2 b B];

BForte1=Table[{B,E1Forte[B] cst},{B,0,2,0.05}];
BForte2=Table[{B,E2Forte[B] cst},{B,0,2,0.05}];
BForte3=Table[{B,E3Forte[B] cst},{B,0,2,0.05}];
BForte4=Table[{B,E4Forte[B] cst},{B,0,2,0.05}];

Fig3=Plot[{E1[B] cst, E2[B] cst, E3[B] cst, E4[B] cst}, {B, 0.1, 2}, 
 PlotRange -> {{0,2},{-15,15}}, AspectRatio -> 0.8, 
	  AxesLabel -> {"B (T)", "10⁵ E (eV)"}];

Fig4=Plot[{E1[B] cst, E2[B] cst, E3[B] cst, E4[B] cst}, {B, 0., 0.3}, 
	  PlotRange -> {{0,0.3},{-2, 2}}, AspectRatio -> 0.8, 
 AxesLabel -> {"B (T)", "10⁵ E (eV)"}];

FigBForteA=ListPlot[{BForte1,BForte2,BForte3,BForte4},
		    PlotStyle->PointSize->0.01];
FigBForteB=ListPlot[{BForte1,BForte2,BForte3,BForte4},
		    PlotStyle->PointSize->0.02];
Show[{Fig3,FigBForteA}]
Show[{Fig4,FigBForteB}]