close all; clear; clc;

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

%% Solve for a range of fields
% Diagonal elements of the Hamiltonian
b1 = @(B) (g1-g2)*muB/2*B; b2 = @(B) (g1+g2)*muB/2*B;
A1 = a/4*h; A2 = sqrt(2)*A1;
B0 = linspace(0,300,1e2); % B0 in Gauss
E = NaN(numel(B0),8);
for Bi = 1:numel(B0)
    B = B0(Bi)./1e4; % Field in Tesla
    H = [0        0        b1(B)+A1 0        0       -A2       0        0       ;...
         0        0        0        b1(B)-A1 0        0        A2       0       ;...
         b1(B)+A1 0        0        0        0        A2       0        0       ;...
         0        b1(B)-A1 0        0        0        0        A2       0       ;...
         0        0        0        0        b2(B)+A1 0        0        0       ;...
         -A2      0        A2       0        0        b2(B)-A1 0        0       ;...
         0        A2       0        A2       0        0       -b2(B)-A1 0       ;...
         0        0        0        0        0        0        0       -b2(B)+A1 ];
    E(Bi,:) = eig(H)/h/1e6; % Energy in MHz
end

%% Plot
figure(1); clf; hold on;
cm = plasma(9);
cm = cm([1 5 2 6 3 7 4 8],:);
for j = 1:8
    plot(B0,E(:,j),'color',cm(j,:));
end
xlabel('Magnetic field (B_0) / G');
ylabel('E/h / MHz');

grid on; box on; axis square;