function allen8(sels)
% ALLEN8(sels)
%
% Magnitude prediction from scale 5 ONLY;
% as a function of the number of stations reporting
%
% INPUT:
%
% sels        0 No distance selection
%             1 Distance selection to within 100 km
%
% Last modified by fjsimons-at-alum.mit.edu, 18.07.2006

defval('sels',0)

stuff=load('/home/fjsimons/EALARMS/wtALL_diagnostics');

% Collect the data
mom=stuff(:,4);
dist=stuff(:,5);
evnu(:,1)=stuff(:,6);
skal(:,1)=stuff(:,7);
skal(:,2)=stuff(:,8);
skal(:,3)=stuff(:,9);
skal(:,4)=stuff(:,10);
skal(:,5)=stuff(:,11);

% Protect identical moments belonging to individual events
mom=mom+evnu/1e10;

% Thresholding division
thresh=stuff(:,3);
difer(diff(thresh))
disp(sprintf('Thresholding divisor %2.2f',thresh(1,1)))

if sels==1
  % Take out those outside of 100 km
  selkt=dist<100;
  mom=mom(selkt);
  skal=skal(selkt,:);
end

clf
ah=axes;
slev=20;
% Just calculate the magnitude errors at SCALE 5
index=5;

axes(ah)

% Get the statistics on the logarithmically scaled data
[s,g]=row2stats(mom,log10(skal(:,index)));
% Now the fit to the data for the low magnitudes
Gmeanl=polyfit(s(1:slev),g.mean(1:slev),1);
Gmeanl=[0.9638   -0.4546]
% Magnitude break point
disp(sprintf('Magnitude breakpoint %4.2f',s(slev)))
yfit=polyval(Gmeanl,s(1:slev));
% Figure out the error in the determination
% Actual seismic magnitude is in s(1:slev)
% Predicted moment is, from the averaged coefficients and the
% regression line; predict ALL with the LOW line
predmoml=1/Gmeanl(1)*g.mean-Gmeanl(2)/Gmeanl(1);
% Make sure this corresponds with what's plotted...
disp(sprintf('m_l = %4.2f logC_%i %+4.1f',...
			 1/Gmeanl(1),index,...
			 -Gmeanl(2)/Gmeanl(1)))
% Plot the error versus the TriNet magnitude
% Predicting the LOW magnitudes with the LOW regression
errll=s(1:slev)-predmoml(1:slev);
% Predicting the HIGH magnitudes with the LOW regression
errlh=s(slev+1:end)-predmoml(slev+1:end);
ell=plot(g.nonnans(1:slev),errll,'s');
hold on
elh=plot(g.nonnans(slev+1:end),errlh,'s');
% Now the fit to the data for the high magnitudes
Gmeanh=polyfit(s(slev+1:end),g.mean(slev+1:end),1);
Gmeanh=[0.6851    0.8416]

disp(sprintf('m_h= %4.2f logC_%i %+4.1f',...
			 1/Gmeanh(1),index,...
			 -Gmeanh(2)/Gmeanh(1)))
yfit=polyval(Gmeanh,s(slev+1:end));
% Figure out the error in the determination
% Predict ALL with the high line
predmomh=1/Gmeanh(1)*g.mean-Gmeanh(2)/Gmeanh(1);
% Plot the error; 
% Predicting the LOW magnitudes with the HIGH regression
errhl=s(1:slev)-predmomh(1:slev);
% Predicting the HIGH magnitudes with the HIGH regression
errhh=s(slev+1:end)-predmomh(slev+1:end);
ehl=plot(g.nonnans(1:slev),errhl,'^');
ehh=plot(g.nonnans(slev+1:end),errhh,'^');
% Labels, etc
xl=xlabel('number of reporting stations');
yl=ylabel('magnitude prediction error');
grid on

% Now the average of both regression lines
avhl=mean([predmoml(:) predmomh(:)],2)';
pav=plot(g.nonnans,s-avhl,'o','Markerf','k','MarkerS',3,'MarkerE','k');

% Cosmetics
hl=2;
set(ehl,'MarkerF',grey,'MarkerS',hl,'MarkerE',grey)
set(elh,'MarkerF',grey,'MarkerS',hl,'MarkerE',grey)
set(ehh,'MarkerF','k','MarkerS',hl,'MarkerE','k')
set(ell,'MarkerF','k','MarkerS',hl,'MarkerE','k')
%delete([ehl elh ehh ell])
%delete(pav)

% This error IS symmetric, it's about their average
[perb,k,l]=errorxy(g.nonnans,s-avhl,[],abs([errhl errhh]-(s-avhl)));
delete(perb)

shrink(gca,1,3)
set(ah(1),'ylim',[-1.3 1.3],'ytick',[-1:1],'xlim',[0 65])
longticks(ah,2)
hold off
% g.nonnans
figdisp
