close all; clear; clc;

%% Simulate the MARY curve for 2-proton system
t = linspace(0,1e-5,1e4+1); dt = t(2)-t(1); % Time vector
J = 0; % J-coupling
B0 = 10.^linspace(0,log10(300),50); % B0 in Gauss
S = NaN(numel(B0),numel(t));
for Bi = 1:numel(B0)
    fprintf('\nRunning point #%d (B0=%.1f)\n',Bi,B0(Bi));
    %S(Bi,:) = liouville_von_neumann_fun_1(t,B0(Bi));
    S(Bi,:) = liouville_von_neumann_fun_2(t,[3 27],J,B0(Bi));
end
fprintf('done\n')
E = exp(-t/2e-6); % Exponential envelope
S2 = bsxfun(@times,S,E);
% Normalised singlet probability as a function of magnetic field
mary = (sum(S2,2).*dt) ./ (sum(E).*dt);

%% Fit with single Lorentzian
lorentz = @(x,p) 1./(1+(x/p(2)*2).^2) .* p(1) + p(3);
w = zeros(size(mary)); w(1) = 1; w(B0>50)=1; % Fit the first point and high field
f = @(p) sqrt(mean(w.*(mary-lorentz(B0',p)).^2));
p0 = [-0.15 60 0.5];
pf = fminsearch(f,p0);

%% Plot
figure(1); clf; hold on;
plot(B0,mary,'o','markersize',2);
plot(B0,lorentz(B0,pf),'k-');
%plot(B0,yfit(B0)+mary(1),'k-');
grid on; box on;
ylabel('singlet probability'); xlabel('magnetic field (B_0) / G');
xlim([min(B0) max(B0)]);
axis square;
title('2-proton RP MARY');
drawnow;

%% Simulate the MARY curve for 1-proton system
% This one is known to behave abnormally
% See "Timmel, C. R. Magnetic Field Effects on Radical Pair Reactions
% (1998). Doctoral thesis. Print."

t = linspace(0,1e-5,1e4+1); dt = t(2)-t(1); % Time vector
%B0 = linspace(0,1000,100)*100; % B0 in Gauss
B0 = 10.^linspace(0,log10(1e5),50); % B0 in Gauss
S = NaN(numel(B0),numel(t));
for Bi = 1:numel(B0)
    fprintf('Running point #%d (B0=%.1f)\n',Bi,B0(Bi));
    S(Bi,:) = liouville_von_neumann_fun_1(t,B0(Bi));
    %S(Bi,:) = liouville_von_neumann_fun_4(t,[3 27],0,B0(Bi));
end
fprintf('done\n')
E = exp(-t/1e-6); % Exponential envelope
S2 = bsxfun(@times,S,E);
% Normalised singlet probability as a function of magnetic field
mary = (sum(S2,2).*dt) ./ (sum(E).*dt);

%% Plot
figure(2); clf; hold on;
plot(B0,mary,'-o','markersize',2);
%plot(B0,yfit(B0)+mary(1),'k-');
grid on; box on;
ylabel('singlet probability');
xlabel('magnetic field (B_0) / G');
axis square;
title('1-proton RP MARY');

%% Plot time evolution of singlet probability
figure(3); clf;
hold on;
plot(t*1e6,S2(1,:));
plot(t*1e6,S2(end,:));
plot(t*1e6,E,'k--')
grid on; box on;
xlabel('time / us'); ylabel('singlet probability');
a = gca; a.XScale = 'log';
title('Time evolution of singlet probability');