%%This program is to compare DMSP boundary crossing data with Wiemer
%%boundary crossing data


for ihem=1:2

if ihem==1
%filename1 = '/home/emery/tieconvec/9807_9906Boundaries_sw.txt';
%filename1 = '/home/emery/tieconvec/1996Boundaries_sw.txt';
l=6;		% location of CRB mlat in older files (no maxVy mlat)
filename1 = '/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_Vsw_3yN.txt';
%filename2 = '/home/emery/tieconvec/1998Jul_1999Jun_Weimer05_CRBs.txt';
%filename2 = '/home/emery/tieconvec/1996Jan_Dec_Weimer01_CRBs.txt';
filename2 = '/home/emery/tieconvec/Boundaries_W05_3yN.txt';

else
filename1 = '/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_Vsw_3yS.txt';
filename2 = '/home/emery/tieconvec/Boundaries_W05_3yS.txt';
end

dataD = dlmread(filename1,'');
dataW = dlmread(filename2,'');
good = 0;
goodx = 0;
above = 0;
below = 0;

% mlat maxVy = 9, CRB=15; Vymax = 8
l = 9;
m = 8;

% i=1:2 is duskNH and 2:2 is dawnNH, 1:1 is all
for i = 1:1:length(dataD)
    if abs(dataW(i,l)) > 50
    if    (abs(dataD(i,l)) +2) >= abs(dataW(i,l)) && abs(dataW(i,l)) >= (abs(dataD(i,l)) -2)
            good = good +1;
	    VyWD = dataW(i,m)/dataD(i,m);
	  if (VyWD > 0.5 && VyWD < 2.0)
	    goodx = goodx + 1;
	  end
        elseif (abs(dataD(i,l)) + 2) < abs(dataW(i,l))
            above = above + 1;		% Weimer is polewards
    elseif abs(dataW(i,l)) < (abs(dataD(i,l)) -2)
            below = below +1;		% Weimer is equatorwards
    else
i
length(dataD)
good
goodx
above
below
        error('problem with algorithm')
            
    end
    end
end

ihem
length(dataD)
disp('good is:')
disp(good)
disp('good peak value (0.5<ratio<2) is:')
disp(goodx)
disp('above is:')
disp(above)
disp('below is:')
disp(below)
disp('total is:')
disp((i)/2)
disp('% good:')
%i+1/2 for dusk, i/2 for dawn, and i for all
disp((good/((i)/2)*100))
disp('% above:')
disp((above/((i)/2)*100))
disp('% below:')
disp((below/((i)/2)*100))

end	% ihem=1:2
