mhdDir='/scratch/bliu17/wireTurbulence/Debug/run_dir/Hist_bar/';
hydroDir='/scratch/bliu17/wireTurbulence/Debug/hydro/Hist_bar/';

hydro_frameN=200;
mhd_frameN=200;

close all

fig=figure

Ind_plot=0;

fn =  strcat(hydroDir,'hist_logRho_frame',num2str(hydro_frameN),'_Box',num2str(9),'0000.curve');
fid = fopen(fn,'r');
formatSpec = '%f %f';
A=textscan(fid,formatSpec,'CommentStyle','#');

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

s=size(X);
msh=s(1);

Yh=zeros(10,msh);

fn =  strcat(mhdDir,'hist_logRho_frame',num2str(mhd_frameN),'_Box',num2str(9),'0000.curve');
fid = fopen(fn,'r');
formatSpec = '%f %f';
A=textscan(fid,formatSpec,'CommentStyle','#');

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

s=size(X);
msm=s(1);

Ym=zeros(5,msm);

for boxN=4:2:10
    Ind_plot = Ind_plot + 1;
    i=0;
    for hydro_frameN=190:200;
    
        i=i+1;
        
        %% hydro
        
        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','#');

        X1=A{1};
        Y=A{2};
        Yh(i,:)=Y;
        
        %B=[X Y];
    end
    
    hAv=sqrt(abs(mean(Yh(:,:))));
    hErr=std(sqrt(abs(Yh(:,:))));
    
    i=0;
    for mhd_frameN=190:200;
         i=i+1;
        %% 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};

        Ym(i,:)=Y2;
     end
     
     mAv=sqrt(abs(mean(Ym(:,:))));
     mErr=std(sqrt(abs(Ym(:,:))));
     
     f = fit(X1,hAv','gauss2');
     f2 = fit(X2,mAv','gauss2');
    
     X1=X1;%*log(10);
     X2=X2;%*log(10);
     
     subplot(2,2,Ind_plot);
     %plot(f,'r-',X1,hAv','k+')
     
     %plot(f2,'g--',X2,mAv','bo')
     
     errorbar(X1,hAv',hErr','--or')
     hold on;
     plot(f,'k-')
     hold on;
     errorbar(X2,mAv',mErr','--.b')     
     hold on; 
     plot(f2,'g-')
     hold off;
    
        %plot(f,'r-',X,Y,'k.',f2,'g--',X2,Y2,'bo')
    
     %get coefficients
     Coeffs=coeffvalues(f);
     sigma=Coeffs(3);%*0.707;
     sigma2=Coeffs(6);

     Coeffs2=coeffvalues(f2);
     sigma3=Coeffs2(3);%*0.707;
     sigma4=Coeffs(6);
    
     xlim([-3 4]);
     ylim([0 0.6]);
     
     xlabel('$s=ln(\rho)$','Interpreter','Latex');
     ylabel('$p_{v}(s)$','Interpreter','Latex');
    
     h=legend('hydro',strcat('Gaussian 2 fit; $\sigma_{1}=$',num2str(sigma),'$\sigma_{2}=$',num2str(sigma2)),'mhd',strcat('Gaussian 1 fit; $\sigma_{1}=$',num2str(sigma3),'$\sigma_{2}=$',num2str(sigma4)));
     set(h,'Interpreter','latex');

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


fclose('all');
