%%This plotting code bins the Boundaries by B strength and then plots the
%%bins as well as saving them in separate files

%% Have 4484 NH Bndys in 3y and 3120 SH Bndys in 3y or 7604 total (or 3802 pairs)
fnameN = '/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_Vsw_3yN.txt';
fnameS = '/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_Vsw_3yS.txt';

colormap(jet)
dataAr = [];
Bvec = [];
mtot = 0
j=1;
for ihem = 1:2
  if (ihem==1)
    filename = fnameN;
  else
    filename = fnameS;
  end
    data = dlmread(filename,'');
    MLT = data(:,16);
    MLAT = data(:,15);
    B = sqrt((data(:,3).^2) + (data(:,4).^2) );
mtot=mtot+length(B)
    for k = 1:length(B)
%        if 0 <= B(k) <3	% This if did not work at all, but all were included
%  Bz 1nT bins: 6%>-1nT 17% 22% 20% 13% 8% 4% 3% 2% 1% 5%<-10nT
% Bz bins: 0to-2(22.4%),-2to-3(21.4%),-3to-4(20.3%),-4to-6(21.5%), <-6(14.4%)
% Mean+/-stdev: -1.30+/-0.46, -2.52+/-0.29, -3.45+/-0.28, -4.85+/-0.56, -9.53+/-3.38 (min -22.252nT)
% num if pairs (mtot=5650/2 = 2825 = 633 605 574 607 406)
% 50% Bz <-3.25
bsm=0.;
blg=-2.;
	if data(k,3)<=bsm && data(k,3)>blg	
	  for n=1:18
            dataAr(j,n) = data(k,n);
          end
%	  Bvec(j) = B(k);
	  Bvec(j) = data(k,3);
%yr325=data(k,1)
%dayn325=data(k,2)
            j=j+1;
        end
    end
    try
%    scatter(dataAr(:,15),dataAr(:,16),20,sqrt((dataAr(:,3).^2) + (dataAr(:,4).^2) ))
    scatter(dataAr(:,15),dataAr(:,16),20,abs(dataAr(:,3)))
%    dlmwrite('/hao/aim2/dmsp/dmsp_idm_4s/3yBoundariesBzSmall.txt', dataAr,'delimiter',' ','-append'); %write boundaries file
%    Bvec = [Bvec; (sqrt((dataAr(:,3).^2) + (dataAr(:,4).^2)))]; %create a vector with all B values
    hold on
    catch
    end
end
% Why is Mean(Bvec)>limits 0-3? - because if was not correct!
(j-1)/2
ratio = (j-1)/mtot
stdv = std(Bvec);
Mean = mean(Bvec);
maxBz = min(Bvec);
%caxis([bsm blg])		% This fails - need real numbers
caxis([0 2])
xlabel('MLAT')
ylabel('MLT')
%title(strcat('Boundaries 1996,9807-9906,2005 0 <= B(nT) < 3 ',' Standard Dev:', num2str(stdv),' Mean:', num2str(Mean)))
title(strcat('Boundaries 1996,9807-9906,2005 -2 < Bz(nT) < 0 ',' Standard Dev:', num2str(stdv),' Mean:', num2str(Mean)))
colorbar

%fplot= ['/hao/aim2/dmsp/dmsp_idm_4s/' '3yBoundariesBzSmall'];
%set(gcf,'PaperPositionMode','auto');
%print(gcf,'-dtiff',fplot);
