u0 = 9.274E-24; % Bohr magneton [J/Tesla]



% 87Rb numerical example

I = 3/2;

a = 1/2*6.835E9; % Hz

h = 6.626E-34; % J*s

gI = 0.0009954;

gJ = 2.002331;



N = 1000;

Bz = linspace(0,0.5,N); % Tesla

Fp = I+1/2;

x = (gI+gJ)*u0*Bz/(a*h*Fp);

Em = zeros(5,N); index = 1;



% m = -2 case requires some special treatment mathematically

m = -2;

Em(index,:) = -a*h/4 - m*gI*u0*Bz + a*h*Fp/2*(1-x);

% plot(Bz,Em(index,:)); hold on;

index = index + 1;

% Upper manifold

for m = -1:2

    Em(index,:) = -a*h/4 - m*gI*u0*Bz + a*h*Fp/2*sqrt(1+2*m*x/Fp+x.^2);

%     plot(Bz,Em(index,:)); hold on;

    index = index + 1;

end



for m = 1:4

    for n = (m+1):5

        plot(Bz,Em(m,:)-Em(n,:)); hold on;

    end

end