% This program illustrates the relation between parameters D, E, gperp, gpar and
% Delta/lambda and V for the Fe(IV) S=1 system. This is the loop script to make 
% figures 4a 4b and 4c. Bracci et al. (Journal Applied Magnetic Resonance 2020)
% This program calls gs_from_crystal_field.m. To work correctly you have to delete 
% (comment) V and Delta variables from gs_from_crystal_field program. In that 
% way V and Delta will be loaded directly from this program.


% Crystal field parameters that will be loaded into gs_from_crystal_field. Rember to delete (or comment) the lines V and Delta inside gs_from_crystal_field
V0=0.1; %Starting V
Delta0=20; %starting Delta
dDelta=10; % increment of Delta at each i loop

n=400; %number of points of x axis
Delta=Delta0;

%loop. At each i loop, Delta is changed. At each j loop, V is changed.
for i=1:n  
    Delta=Delta+dDelta;   
    Delta_array(i,1)=Delta;
    V=V0;
    gs_from_crystal_field
    dV=0.5*Delta/2; %increment of V
    for j=1:3
        V_array(i,j)=V;
        gs_from_crystal_field
        E_array(i,j)=E;
        D_array(i,j)=D;
        gx_array(i,j) = gxx;
        gy_array(i,j) = gyy;
        gz_array(i,j) = gzz;
        V=V+dV; 
    end
end

ratioDelta=Delta_array./lambda; % Delta/lambda. lambda as specified in gs_from_crystal_field

figure(401) %figure 4a
title('D and E vs \Delta/\lambda')
hold on
plot(ratioDelta,D_array(:,1),'color',[0, 0.4470, 0.7410],'LineWidth', 2);
plot(ratioDelta,D_array(:,2),'color',[0.4660, 0.6740, 0.1880],'LineWidth', 2);
plot(ratioDelta,D_array(:,3),'color',[0.8500, 0.3250, 0.0980],'LineWidth', 2);
plot(ratioDelta,E_array(:,1),'color',[0, 0.4470, 0.7410],'LineWidth', 2);
plot(ratioDelta,E_array(:,2),'color',[0.4660, 0.6740, 0.1880],'LineWidth', 2);
plot(ratioDelta,E_array(:,3),'color',[0.8500, 0.3250, 0.0980],'LineWidth', 2);
text(1,160, 'D','FontSize',15,'HandleVisibility','off')
text(1,20, 'E','FontSize',15,'HandleVisibility','off')
ylabel('[cm^-^1]')
xlabel('\Delta/\lambda')
axis tight
ylim ([-20 200]);
plots=get(gca, 'Children');
legend(plots([3,2,1]), {'V = 0','V = 0.25 \Delta','V = 0.5 \Delta'}); 
set(gca,'TickDir', 'out','box', 'off', 'Linewidth', 1,'FontSize',15,'xtick', [0.1 2 4 6 8 10]); 
legend('box', 'off', 'FontSize', 15);

figure(402) %figura 4b
title('E/D vs \Delta/\lambda')
hold on
plot(ratioDelta,abs(E_array(:,1)./D_array(:,1)),'color',[0, 0.4470, 0.7410],'LineWidth', 2);
plot(ratioDelta,abs(E_array(:,2)./D_array(:,2)),'color',[0.4660, 0.6740, 0.1880],'LineWidth', 2);
plot(ratioDelta,abs(E_array(:,3)./D_array(:,3)),'color',[0.8500, 0.3250, 0.0980],'LineWidth', 2);
ylabel('[cm^-^1]')
xlabel('\Delta/\lambda')
axis tight
ylim ([-0.01 0.23]);
plots=get(gca, 'Children');
legend(plots([3,2,1]), {'V = 0','V = 0.25 \Delta','V = 0.5 \Delta'}); 
set(gca,'TickDir', 'out','box', 'off', 'Linewidth', 1,'FontSize',15,'xtick', [0.1 2 4 6 8 10]); 
legend('box', 'off', 'FontSize', 15);


figure(403) %figura 4c
title('g_/_/ and g_\perp vs \Delta/\lambda')
plot(ratioDelta,gx_array(:,1),ratioDelta,gy_array(:,1),ratioDelta,gz_array(:,1),'color',[0, 0.4470, 0.7410],'LineWidth', 2)
hold on
plot(ratioDelta,gx_array(:,2),ratioDelta,gy_array(:,2),ratioDelta,gz_array(:,2),'color',[0.4660, 0.6740, 0.1880],'LineWidth', 2)
plot(ratioDelta,gx_array(:,3),ratioDelta,gy_array(:,3),ratioDelta,gz_array(:,3),'color',[0.8500, 0.3250, 0.0980],'LineWidth', 2)
text(4,1.70, 'g_/_/','FontSize',15,'HandleVisibility','off')
text(4,2.50, 'g_\perp','FontSize',15,'HandleVisibility','off')
ylabel('g-values')
xlabel('\Delta/\lambda')
xlim ([0.1 10])
ylim ([0.5 2.8]);
plots=get(gca, 'Children');
legend(plots([7,4,1]), {'V = 0','V = 0.25 \Delta','V = 0.5 \Delta'}); 
set(gca,'TickDir', 'out','box', 'off', 'Linewidth', 1,'FontSize',15,'xtick', [0.1 2 4 6 8 10]); 
legend('box', 'off', 'FontSize', 15);