%----------------------------------
% This script plots Fig. 4b
%----------------------------------

% load data for \beta A/N - \beta U_0/N
% calculated from Einstein crystal simulations
d = load('ars.dat');
b = d(1:end,1);
ars = d(1:end,2);
p = 0.0001;

% functions to calculate U_0 and d U_0/ d\beta
sw = -50;
e = @(x) -1*( (59.9) + 2*sw * (2-1./x));
de = @(x) -2*sw./x.^2;

% load equation of state and <U_0>/N
% averages from NPT simulations of the crystal
dd = dlmread('DBI_isobar_P0.0001.dat','',1,0);

% beta
bb = dd(1:5,1);

% <U_0>/N
unp = dd(1:5,3);

% density
rho = dd(1:5,2);

% fit <U_0>/N for interpolation
[xData, yData] = prepareCurveData( bb, unp );
ft = 'pchipinterp';
[unpFit, gof] = fit( xData, yData, ft );

% fit rho for interpolation
[xData2, yData2] = prepareCurveData( bb, rho );
ft = fittype( 'smoothingspline' );
opts = fitoptions( 'Method', 'SmoothingSpline' );
opts.SmoothingParam = 0.999988546410448;
[rhoFit, gof] = fit( xData2, yData2, ft, opts );



%% \beta\mu assuming all bonds are active

mus = load('mus.dat');

% interpolate beta range
betas = (b(1):0.001:b(end));

% interpolate rho and U_0 at each beta value
rhoF = rhoFit(betas);
unpF = unpFit(betas);

% integrate starting from \mu at the lowest \beta
mu1=[];
start=1;
for i=1:numel(betas)
    if i ==1
        mu1 = [mu1;mus(1,2)];
    else
        % thermodynamic integration along isobar
        m = mus(1,2)+ trapz(betas(1:i),unpF(1:i))+trapz(betas(1:i),p./rhoF(1:i)) + trapz(betas(1:i),betas(1:i).*de(betas(1:i)));
        mu1=[mu1;m];
    end
end

% plot resulting \beta\mu_x
plot(betas,mu1,'lineWidth',2, 'color',[69/255, 171/255, 153/255])
hold on

%% take frac into account

betas = (0.3:0.0001:0.7);
% load average fraction of bonds active for Patch 4, as averaged from NPT
% simulations
fdata = load('c4frac.dat');
% betas
bf = fdata(:,1);
% fractions
frac = fdata(:,2);

% fit the fraction data to interpolate
[xData, yData] = prepareCurveData( bf, frac );
ft = 'pchipinterp';
[fitresult, gof] = fit( xData, yData, ft, 'Normalize', 'on' );
fracF = fitresult(betas);

rhoF = rhoFit(betas);
unpF = unpFit(betas);

% integrate starting from \mu at the lowest \beta
mu1=[];
start=1;
for i=1:numel(betas)
    if i ==1
        mu1 = [mu1;mus(1,2)];
    else
        % thermodynamic integration along isobar
        % with <d U_0/ d\beta>  = frac * dU_0/d\beta
        m = mus(1,2)+ trapz(betas(1:i),unpF(1:i))+trapz(betas(1:i),p./rhoF(1:i)) + trapz(betas(1:i),fracF(1:i)'.*betas(1:i).*de(betas(1:i)));
        mu1=[mu1;m];
    end
end

% plot resulting \beta\mu_x
plot(betas,mu1,'k','lineWidth',2);
hold on

% plot \beta\mu_x calculated from individual Einstein crystal simulations
scatter(mus(:,1),mus(:,2),'filled','MarkerEdgeColor',[0.79,0.4,0.47],'MarkerFaceColor',[0.79,0.4,0.47])


%% figure properties

xlabel('\beta')
ylabel('\beta\mu')
xline( 0.5181,'--','lineWidth',1.5)
set(gca,'FontSize',18)
annotation('textbox',[0.53,0.99,0,0], 'string', '\beta_{tp}')
saveas(gcf, 'fig4b','epsc');