clear all: close all; numOfBoxes=5; N=16*(numOfBoxes^3) PdVector=[0.01,0.02,0.03]; PdLabel={'0.01','0.02','0.03'}; cM={'ko','ko','ko','ko','ko'}; cM1={'r-o','g-o','b-o','m-o','c-o'}; cM2={'r-s','g-s','b-s','m-s','c-s'}; cL={'k-','k-','k-','k-','k-'}; topDir = ['./gardnerHS1/']; Npds = numel(PdVector); data1 = []; data2 = []; data3 = []; idx = []; for i=1:Npds PdVector(i) datFile1 = [ topDir 'd_dab_betaP_phi_N' num2str(N) '_s' num2str(PdVector(i)) '.dat'] data1{i} = load(datFile1); datFile2 = [ topDir 'skewDab_phi_N' num2str(N) '_s' num2str(PdVector(i)) '.dat'] data2{i} = load(datFile2); data3{i} = [data1{i}(:,:), data2{i}(:,2)]; end data4=[]; for i=1:Npds idx{i}=(data3{i}(:,end)==max(data3{i}(:,5))); Pg = data3{i}(idx{i},2); idx_M = find(idx{i}==1); idx_L = find(idx{i}==1)-1; idx_U = find(idx{i}==1)+1; errorPg_L = ( data3{i}(idx_M,2) - data3{i}(idx_L,2) )/2; errorPg_U = ( data3{i}(idx_U,2) - data3{i}(idx_M,2) )/2; data4 = [ data4; Pg, errorPg_L, errorPg_U]; end PdCurve = min(PdVector):0.001:max(PdVector); PgData = 1./data4(:,1)'; PgFit = polyfit(PdVector,PgData,1); PgCurve = 1./polyval(PgFit,PdCurve); dim1=500; dim2=400; fig1=figure('position',[0,0,dim1,dim2]); errorbar(PdVector,data4(:,1),data4(:,2),data4(:,3),cM{i}); hold on; plot(PdCurve,PgCurve,cL{i}); title([ 'HS1 betaP_G vs s nParticles = ' num2str(nParticles)]); xlabel('s') ylabel('$${\beta P_{G} \bar{\sigma}^3}$$','Interpreter','Latex'); %axis equal xlim([0.0075,0.0325]) ylim([40,220])