mhdDir='/gpfs/fs2/scratch/bliu17/wireTurbulence/V2plot/mhd/Hists2/';
hydroDir='/gpfs/fs2/scratch/bliu17/wireTurbulence/V2plot/hydro/Hists2/';

hydro_frameN=199;
mhd_frameN=194;

close all

fig=figure

Ind_plot=0;
for boxN=4:2:10
    
    %% hydro
    Ind_plot = Ind_plot + 1;
    
    fn =  strcat(hydroDir,'hist_logRho_frame',num2str(hydro_frameN),'_Box',num2str(boxN),'0000.curve');
    
    fid = fopen(fn,'r');

    formatSpec = '%f %f';

    A=textscan(fid,formatSpec,'CommentStyle','#');

    X=A{1};
    Y=A{2};

    B=[X Y];

     
    %% mhd
    fn2 =  strcat(mhdDir,'hist_logRho_frame',num2str(mhd_frameN),'_Box',num2str(boxN),'0000.curve');
    
    fid2 = fopen(fn2,'r');

    formatSpec = '%f %f';

    A2=textscan(fid2,formatSpec,'CommentStyle','#');

    X2=A2{1};
    Y2=A2{2};

    B2=[X2 Y2];
    
    C=B;
    C2=B2;

    %plot(C(:,1),C(:,2),'.');
    X=C(:,1)*log(10);
    Y=C(:,2);
    
    X2=C2(:,1)*log(10);
    Y2=C2(:,2);
    
    f = fit(X,Y,'gauss1');
    f2 = fit(X2,Y2,'gauss1');
    
    subplot(2,2,Ind_plot);
    plot(f,'r-',X,Y,'k+')
    
    hold on;
    
    plot(f2,'g--',X2,Y2,'bo')
    hold off;
    
    %plot(f,'r-',X,Y,'k.',f2,'g--',X2,Y2,'bo')
    
    %get coefficients
    Coeffs=coeffvalues(f);
    sigma=Coeffs(3)*0.707;

    Coeffs2=coeffvalues(f2);
    sigma2=Coeffs2(3)*0.707;
    
    ylim([0 0.18]);
    
    xlabel('$s=ln(\rho)$','Interpreter','Latex');
    ylabel('$p_{v}(s)$','Interpreter','Latex');
    
    h=legend('hydro',strcat('Gaussian 1 fit; $\sigma=$',num2str(sigma)),'mhd',strcat('Gaussian 1 fit; $\sigma=$',num2str(sigma2)));
    set(h,'Interpreter','latex');

    title(strcat('density PDF: box ',num2str(boxN)));
    

end

fclose('all');

