%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Zeeman and Paschen-Back splitting of a fine structure level and shift of the hyperfine structure sublevels,
%  the interaction matrix is set up and diagonalized as a whole.
%
% Arbitrary spin structure and comparison with Breit-Rabi formula.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Breit
addpath C:\Users\Philippe\Physik\Appendix\Appendix;			Constant;
figure(1);							clf;
B=[0:.1:5]'*1e-4;

% I =1;		L =0;		J =1/2;		S =1/2;		Aj=152.137e6;			Bj=-.1;						%6Li 2S1/2
% %E1=BreitRabiFormel(B,I,J,L,S,Aj,Bj);
% E =PaschenBackHFSMatrix(B,I,J,L,S,Aj,Bj);
% 
% subplot(221),plot(1e4*B,1e-6*E','-k','LineWidth',1.3)%,1e4*B,1e-6*E1,'-r'
% set(gca,'xLim',[0 max(B*1e4)],'FontSize',12,'LineWidth',1.3);
% xlabel('\fontsize{12}{\itB}   (G)')
% ylabel('\fontsize{12}\Delta{\it\nu}   (MHz)')


I =3/2;		L =1;		J =3/2;		S =1/2;		Aj=-(I+.5)/2*1.17e6;		Bj=0;						%6Li 2P3/2
E =PaschenBackHFSMatrix(B,I,J,L,S,Aj,Bj);

subplot(222),plot(1e4*B,1e-6*E','-k','LineWidth',1.3)
set(gca,'xLim',[0 max(B*1e4)],'FontSize',16,'LineWidth',1.3);
xlabel('\fontsize{16}{\itB}   (G)')
ylabel('\fontsize{16}\Delta{\it\nu}   (MHz)')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Zeeman and Paschen-Back splitting of a fine structure level and shift of the hyperfine structure sublevels,
%  this routine uses the Breit-Rabi formula. The labeling is define by
%
% mImJ= 1  2   .. 2I+1 2I+2 .. 2I+1+2J+1
% mJ  =-J -J   .. -J   -J+1 .. +J
% mI  =-I -I+1 .. +I   -I	 .. +I

function vG = BreitRabiFormel(B,I,J,L,S,Aj,Bj)

Constant;
vHFS=2*Aj;
gJ=2;
F =1;
mF=[F:-1:-F];						[mF_,B_]=meshgrid(mF,B);
x=uB*gJ*B_/(2*pi*hP)/vHFS;			vG1=-vHFS/2/(2*I+1) - vHFS/2*sqrt(1+4*mF_/(2*I+1).*x+x.^2);
vG1=vG1;

F =2;
mF=[F:-1:-F];						[mF_,B_]=meshgrid(mF,B);
x=uB*gJ*B_/(2*pi*hP)/vHFS;			vG2=-vHFS/2/(2*I+1) + vHFS/2*sqrt(1+4*mF_/(2*I+1).*x+x.^2);
vG2=vG2;
vG=[vG1 vG2];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Zeeman and Paschen-Back splitting of a fine structure level and shift of the hyperfine structure sublevels,
%  this routine settles the matrix and gets its eigenvalues.

function E = PaschenBackHFSMatrix(B,I,J,L,S,Aj,Bj)

for k=1:length(B)
	for mImJ1=1:(2*I+1)*(2*J+1)
		for mImJ2=1:(2*I+1)*(2*J+1)
%			mJ1=floor((mImJ1-1)/(2*I+1))-J;			mI1=mImJ1-I-1-(J+mJ1)*(2*I+1);
%			mJ2=floor((mImJ2-1)/(2*I+1))-J;			mI2=mImJ2-I-1-(J+mJ2)*(2*I+1);
			mI1=floor((mImJ1-1)/(2*J+1))-I;			mJ1=mImJ1-J-1-(I+mI1)*(2*J+1);
			mI2=floor((mImJ2-1)/(2*J+1))-I;			mJ2=mImJ2-J-1-(I+mI2)*(2*J+1);
			Mhf(mImJ1,mImJ2)=PaschenBackHFSCoefficient(B(k),I,J,L,S,mI1,mJ1,mI2,mJ2,Aj,Bj);
		end
	end
	E(:,k)=sort(eig(Mhf));
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Zeeman and Paschen-Back splitting of a fine structure level and shift of the hyperfine structure sublevels,
%  this routine calculates one matrix element.

function Chf = PaschenBackHFSCoefficient(B,I,J,L,S,mI1,mJ1,mI2,mJ2,Aj,Bj)

Constant;
gJ=(3*J*(J+1)+S*(S+1)-L*(L+1))/2/J/(J+1);		gI=gJ*uK/uB;

denom=2*I*(2*I-1)*2*J*(2*J-1)+1e-9;
Im =sqrt(I*(I+1)- mI2   *(mI2-1));				Jm =sqrt(J*(J+1)- mJ2   *(mJ2-1));
Ip =sqrt(I*(I+1)- mI2   *(mI2+1));				Jp =sqrt(J*(J+1)- mJ2   *(mJ2+1));
Ipp=sqrt(I*(I+1)-(mI2+1)*(mI2+2));				Jpp=sqrt(J*(J+1)-(mJ2+1)*(mJ2+2));
Imm=sqrt(I*(I+1)-(mI2-1)*(mI2-2));				Jmm=sqrt(J*(J+1)-(mJ2-1)*(mJ2-2));
if mI1==mI2       & mJ1==mJ2
	Chf=(Aj+3*Bj/denom)*mI2*mJ2 - Bj*2*I*(I+1)*J*(J+1)/denom + (gJ*uB*mJ2-gJ*uK*mI1)*B/(2*pi*hP) + 6*Bj/denom*((mI2*mJ2)^2+(Im*Jp)^2/4+(Ip*Jm)^2/4);
elseif mI1==mI2+1 & mJ1==mJ2-1
	Chf=(Aj+3*Bj/denom)/2*Ip*Jm + 6*Bj/denom*Ip*Jm*(mI1*mJ1+mI2*mJ2)/2;
elseif mI1==mI2-1 & mJ1==mJ2+1
	Chf=(Aj+3*Bj/denom)/2*Im*Jp + 6*Bj/denom*Im*Jp*(mI1*mJ1+mI2*mJ2)/2;
elseif mI1==mI2+2 & mJ1==mJ2-2
	Chf=6*Bj/denom*Ip*Jm*Ipp*Jmm/4;
elseif mI1==mI2-2 & mJ1==mJ2+2
	Chf=6*Bj/denom*Im*Jp*Imm*Jpp/4;
else
	Chf=0;
end