close all; clear; clc;

%% Constants
ge = 2.0023; gn = ge+0.3; %5.585695; % Electron and proton g-values
h = 6.62607004e-34; % Planck's constant
muB = 9.27400968e-24; muN = muB; % Joule per Tesla
j = -100e6; % Hertz, J-coupling

%% Solve for a range of fields
% Diagonal elements of the Hamiltonian
d1 = @(B) + ge*muB/2*B + gn*muN/2*B - j/4*h;
d2 = @(B) + ge*muB/2*B - gn*muN/2*B + j/4*h;
d3 = @(B) - ge*muB/2*B + gn*muN/2*B + j/4*h;
d4 = @(B) - ge*muB/2*B - gn*muN/2*B - j/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*j/2 0     ;...
          0     h*j/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

%% Plot
figure(1); clf; hold on;
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');
plot(B0,E2,'k--');

grid on; box on; axis square;