indx=(plotOutData(:,2)>=r0ini & plotOutData(:,2)<=r0end); 
plot2OutData=plotOutData(indx,:);

%figure(xaxis)
%hold off

Bms=log(plot2OutData(:,8)+plot2OutData(:,6)+plot2OutData(:,7))/log(10);

if xaxis==2
yBms=log(plot2OutData(:,xaxis))/log(10);
else
yBms=Bms;
end

Bms=bmscale*Bms;



yyBms=yBms(yBms>-Inf);
yyBms=sort(yyBms);
yBmsm=max([yyBms(1), min(yBms)]);
yBmsM=max(yBms);
yBinsize=Binsize;
yBinsize=(yBmsM-yBmsm)/Npts;


yBmsBins=yBmsm:yBinsize:yBmsM;

yzpms=[];
y1nmbr=[];
z1nmbr=[];
p1nmbr=[];
yzpstd=[];
y1std=[];
z1std=[];
p1std=[];
for i=yBmsBins
    ind1=(yBms>=i & yBms<i+yBinsize);
    yzpms=[yzpms, mean(Bms(ind1))];
    y1nmbr=[y1nmbr, mean(plot2OutData(ind1,9))];
    z1nmbr=[z1nmbr, mean(plot2OutData(ind1,10))];
    p1nmbr=[p1nmbr, mean(plot2OutData(ind1,11))];
    yzpstd=[yzpstd, std(Bms(ind1))];
    y1std=[y1std, std(plot2OutData(ind1,9))];
    z1std=[z1std, std(plot2OutData(ind1,10))];
    p1std=[p1std, std(plot2OutData(ind1,11))];
end

if xaxis==2
    plot(yBmsBins+yBinsize/2,yzpms,'ko',yBmsBins+yBinsize/2,yzpms,'y-','linewidth', 2.5)

    hold on
    for j=1:length(yBmsBins)
    plot([yBmsBins(j)+yBinsize/2 yBmsBins(j)+yBinsize/2],...
        [yzpms(j)+yzpstd(j) yzpms(j)-yzpstd(j)], 'ys')
    end

end

plot(yBmsBins+yBinsize/2,y1nmbr,'ko',yBmsBins+yBinsize/2,z1nmbr,'ko',...
    yBmsBins+yBinsize/2,p1nmbr,'ko',...
    yBmsBins+yBinsize/2,y1nmbr,'g-',yBmsBins+yBinsize/2,z1nmbr,'b-',...
    yBmsBins+yBinsize/2,p1nmbr,'r-', ...
    'linewidth', 2.5)

hold on
if stdon==1
hold on
for j=1:length(yBmsBins)
    plot([yBmsBins(j)+yBinsize/2 yBmsBins(j)+yBinsize/2],...
        [z1nmbr(j)+z1std(j) z1nmbr(j)-z1std(j)], 'b-','linewidth',1.5)
end

for j=1:length(yBmsBins)
    plot([yBmsBins(j)+yBinsize/2 yBmsBins(j)+yBinsize/2],...
        [y1nmbr(j)+y1std(j) y1nmbr(j)-y1std(j)], 'g.',...
        [yBmsBins(j)+yBinsize/2 yBmsBins(j)+yBinsize/2],...
        [p1nmbr(j)+p1std(j) p1nmbr(j)-p1std(j)], 'r-')
end

end

axis([yBmsm yBmsM 0 1.1*max([y1nmbr,z1nmbr,p1nmbr]+[y1std,z1std,p1std])])
