% domain wall free energy of the ANNNI model

% constraint boundary condition

close all;

modeltype = '2DXAT';

widthpix = 400;
heightpix = 300;

leftmargin = 70;
bottommargin = 50;
topmargin = 10;
rightmargin = 5;

widthtotal = widthpix+leftmargin+rightmargin;
heighttotal = heightpix+topmargin+bottommargin;
left = leftmargin/widthtotal;
bottom = bottommargin/heighttotal;
width = widthpix/widthtotal;
height = heightpix/heighttotal;

psize = 18; % font size
msize = 8; % marker size

fig = figure('rend','painters','pos',[100 100 widthtotal heighttotal]);


sizes = [12, 16, 20, 24, 28];

kappa = 0.6;
hs = 0;

n=numel(sizes);
rpp=(1:n)';
corders = [(rpp-1)/(n-1), zeros(n,1), 1 - (rpp-1)/(n-1)];


criticalT = 2/log(1+sqrt(2));

lgds = {};
ngpoint = [];
allresults = {};

panelA = axes('Position', [left, bottom, width, height]);

for rp=1:length(sizes)
    asize = sizes(rp);
    fname = strcat('../data/',modeltype,'/bdf_', num2str(kappa), '_', num2str(hs) , '_', num2str(asize), '.dat');
    if ~exist(fname, 'file')
        continue
    else
        data = dlmread(fname);
    end
    fes = [1./data(:,2), ( data(:,5) - data(:,6) ) ./ data(:,2) * asize ];
    plot(fes(:,1), fes(:,2), 'color', corders(rp,:));
    allresults{end+1} = fes;
    lgds{end+1} = strcat('$L=', num2str(asize+2),'$');
    if 1==rp
        hold on
    end
end

% extrapolations
checkTs = linspace(0.6, 0.915).';
divideT = 0.89*0.95;

rst = ones(numel(sizes), 1);  % L^(-2)
rst2 = ones(numel(sizes), 1); % L^(-1)
extrapolatedFs = zeros(length(checkTs), length(sizes));

expns = 2*ones( numel(checkTs), 1);
expns( checkTs>divideT) = 2 - (checkTs(checkTs>divideT)-divideT)./(0.89-divideT); % fix exponent B

for rp=1:length(sizes)
    extrapolatedFs(:, rp) = interp1( allresults{rp}(:,1), allresults{rp}(:,2), checkTs);
end

for rp=1:length(checkTs)
    [b, bint] = regress(extrapolatedFs(rp,:).', [sizes.^(-2).', ones(numel(sizes), 1) ]);
    rst(rp) = b(2);
    [b, bint] = regress(extrapolatedFs(rp,:).', [sizes.^(-1).', ones(numel(sizes), 1) ]);
    rst2(rp) = b(2);
end
plot(checkTs, rst, '-k');
plot(checkTs(checkTs > 0.85 & checkTs <= 0.895), rst2(checkTs > 0.85 & checkTs <= 0.895), '--k');

% extrapolate peaks
peaks = [1.42, 0.1116; 1.2, 0.1038; 1.1, 0.08065; 1.04, 0.065; 0.89, 0]; % only four points does not fit. add zero point
valleys = [1.12, -0.1547; 1.035, -0.1003; 0.9925, -0.07096; 0.97, -0.05308; 0.955, -0.041];

ppoly = polyfit(peaks(:,1), peaks(:,2), 2);
vpply = polyfit(valleys(:,1), valleys(:,2), 2);
xs = linspace(0.89, 1.3);
pys = polyval(ppoly, xs);
vys = polyval(vpply, xs);
plot(xs(pys>=0), pys(pys>=0), ':k');
plot(xs(vys<=0), vys(vys<=0), ':k');
plot([0.6,1.3], [0,0], '-k');
xlabel('$T$', 'interpreter', 'latex');
ylabel('$\tau_L^{(+/-)}(T)$', 'interpreter', 'latex');
xlim([0.6,1.3]);
ylim([-0.2,0.4]);
xticks(0.6:0.2:1.3);
yticks(-0.2:0.2:0.4);

set(gca, 'fontname', 'times new roman', 'fontsize', psize);

%% add inset
insetA = axes('Position', [left+0.55*width, bottom+0.6*height, 0.45*width, 0.4*height]);
xspeak = 1./[18;22;24;30];
xsvalley = 1./[14;18;22;24;30];
firstzeros = [14, 0.9838; 18, 0.9542; 22, 0.9393; 26, 0.9303; 30, 0.9249];

plot( xspeak, peaks(1:end-1,1), 'xk', 'MarkerFaceColor', 'k', 'markersize', msize);
hold on;
plot( xsvalley, valleys(:,1), 'sk', 'MarkerFaceColor', 'k', 'markersize', msize);
plot( 1./firstzeros(:,1), firstzeros(:,2), 'ok', 'MarkerFaceColor', 'k', 'markersize', msize);


ppoly = polyfit([xspeak(2:end); 0.015; 0], [peaks(2:end-1,1); 0.93; 0.89], 2);
vpply = polyfit([xsvalley; 0], [valleys(:,1); 0.89], 2);
zpoly = polyfit(1./firstzeros(:,1), firstzeros(:,2), 2);

xs1 = linspace(0, 0.055);
xs2 = linspace(0, 0.071);
pys = polyval(ppoly, xs1);
vys = polyval(vpply, xs2);
zys = polyval(zpoly,xs2);

plot(xs1, pys, ':k');
plot(xs2, vys, ':k');
plot(xs2, zys, ':k');

xlim([0, 0.08]);
ylim([0.85, 1.5]);

xticks([0,0.03,0.06]);
yticks([0.9, 1.1, 1.3]);
ylabel('$T$', 'interpreter', 'latex');

xx = gca;
xx.XRuler.TickLength = [0.02, 0];
xx.YRuler.TickLength = [0.03, 0];

set(gca, 'fontname', 'times new roman', 'fontsize', psize-2);
annotation('textbox',[0.6, 0.47, 0.1, 0.1], 'String', '$1/L$', 'interpreter', 'latex', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);

print(strcat('../figures/artfig_ANNNIdomain', '.eps'), '-depsc');