close all; clear; clc;

%% Constants
ge = 2.0023; gn = 5.585695; % Electron and proton g-values
h = 6.62607004e-34; % Planck's constant
muB = 9.27400968e-24; muN = 5.05078324e-27; % Joule per Tesla
a = 200e6; % Hertz, Hyperfine coupling constant

%% Solve for a range of fields
% Diagonal elements of the Hamiltonian
d1 = @(B) + ge*muB/2*B + gn*muN/2*B + a/4*h;
d2 = @(B) + ge*muB/2*B - gn*muN/2*B - a/4*h;
d3 = @(B) - ge*muB/2*B + gn*muN/2*B - a/4*h;
d4 = @(B) - ge*muB/2*B - gn*muN/2*B + a/4*h;
B0 = linspace(0,300,1e2); % B0 in Gauss
E1 = NaN(numel(B0),4); E2 = NaN(numel(B0),4);
for Bi = 1:numel(B0)
    B = B0(Bi)./1e4; % Field in Tesla
    H1 = [d1(B) 0     0     0     ;...
          0     d2(B) h*a/2 0     ;...
          0     h*a/2 d3(B) 0     ;...
          0     0     0     d4(B) ];
    H2 = [d1(B) 0     0     0     ;...
          0     d2(B) 0     0     ;...
          0     0     d3(B) 0     ;...
          0     0     0     d4(B) ];
    E1(Bi,:) = eig(H1)/h/1e6; % Energy in MHz
    E2(Bi,:) = eig(H2)/h/1e6; % Energy in MHz
end

% Field at which AlphaAlpha - BetaAlpha transition occurs
Bt = interp1(abs(E1(:,1)-E1(:,3)),B0,600);

%% Plot
figure(1); clf; hold on;
plot(B0,E2,'k--');
cm = plasma(5);
for j = 1:4
    plot(B0,E1(:,j),'color',cm(j,:));
end
xlabel('Magnetic field (B_0) / G');
ylabel('E/h / MHz');
grid on; box on; axis square;