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

% load data for \beta A/N - \beta U_0/N
% calculated from Einstein crystal simulations
d = load('ars.dat');

% betas
b = d(:,1);
% \beta A/N - \beta U_0/N
ars = d(:,2);

% pressure
p = 0.0001;

% functions to calculate U_0 and d U_0/ d\beta
e = @(x) -1*(15.8/2*(1+tanh((1./x-1.9)/0.05)) + 59.9-15.8);
de = @(x) (1./x).*(15.8)/(2*0.05).*sech((1./x-1.9)/0.05).^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(3:15,1);

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

% density
rho = dd(3:15,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.999997275425068;
[rhoFit, gof] = fit( xData2, yData2, ft );


%% \beta\mu assuming all bonds are active

% load \beta\mu_x data
mus = load('mus.dat');

% interpolate beta range
betas = (mus(1,1):0.0001:mus(end,1));

% interpolate U_0 and density at each beta value
unpF = unpFit(betas);
rhoF = rhoFit(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),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

%% betamu using fraction of bonds active

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

% fit the fraction data to interpolate
[xData, yData] = prepareCurveData( bf, frac );
ft = 'pchipinterp';
[fracFit, gof] = fit( xData, yData, ft, 'Normalize', 'on' );
fracF = fracFit(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).*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.5291,'--','lineWidth',1.5)
set(gca,'FontSize',18)
annotation('textbox',[0.7,0.99,0,0], 'string', '\beta_{tp}')

saveas(gcf, 'fig4a','epsc');