function scatter_plot = scatter_plot(x, y, minax, maxax, head, head2)  
% matlab function to plot scatter plots, one plot per page
% and optionally print some basic statistics in the corner
% x = independent variable, y = dependent variable, min and max are for axes, 
% head = title string, head2 = second header
% written by Ian Renfrew, 1999
% function scatter_plot = scatter_plot(x, y, minax, maxax, head, head2)  

 plot(x, y, 'b+')	
hold on
 plot([minax maxax], [minax maxax], 'k-')
 axis([minax maxax minax maxax])
 set(gca,'DataAspectRatio', [1 1 1])
 set(gca, 'fontsize', 12)
% set(gca, 'xtick', [-14000:100:1400])
 set(gca, 'ytick', get(gca,'xtick'))          %equalises axes ticks
 xlabel('observed data')
 ylabel('model data')
 title([head head2])

 dhtext = (maxax-minax)/12;
% text(minax+dhtext, maxax-dhtext, '(a)','fontsize', 14); 

%correlation coefficients - results in a 2x2 matrix
ccxy = corrcoef([x y]);

%linear regression for (predicted, predictors), where predictors are indep.
%regressions of the form y = a + bx 
%reg(1) = intercept = a; reg(2) = slope = b
 reg_yx = regress_test(y,[ones(size(x)) x]);
 fit = reg_yx(1) + reg_yx(2)*x;
 plot(x, fit,'b- ', 'Linewidth', 1.5)

%peter's error statistics
 biaserr = mean(y-x); 
 slopeerr = std(fit-biaserr,1);
 ranerr = std(fit-(y-x),1);
 toterr = sqrt(biaserr.^2 + slopeerr.^2 + ranerr.^2);
 % root-mean-square error, added 16 May 2008
 rmserr = sqrt(mean((y-x).^2));                         

%add text table
 str2= {'obs.    model'};
 str3= [' mean  ' num2str(mean(x),'%8.3f') '  ' num2str(mean(y),'%8.3f')];
 str4= [' std    ' num2str(std(x),'%8.3f') '  ' num2str(std(y),'%8.3f')];
 str5= [' max   ' num2str(max(x),'%8.3f') '  ' num2str(max(y),'%8.3f')];
 str6= [' min    ' num2str(min(x),'%8.3f') '  ' num2str(min(y),'%8.3f')];
 str7= ['rcor       ' num2str(ccxy(1,2),'%8.3f')];
 str8= ['rsqrd      ' num2str(ccxy(1,2).^2,'%8.3f')];
 str9= ['intercept  ' num2str(reg_yx(1),'%8.3f')];
 str10= ['slope     ' num2str(reg_yx(2),'%8.3f')];
 str11= ['bias error ' num2str(biaserr,'%8.3f')];
 str12= ['slope error ' num2str(slopeerr,'%8.3f')];
 str13= ['random error ' num2str(ranerr,'%8.3f')];
 str14= ['total error  ' num2str(toterr,'%8.3f')];
 str15= ['r.m.s. error  ' num2str(rmserr,'%8.3f')];
 
print_text =1;
if (print_text ==1)
 dtext = (maxax-minax)/24;
 text(minax+3*dtext, maxax-1*dtext, str2)
 text(minax, maxax-2*dtext, str3)
 text(minax, maxax-3*dtext, str4)
 text(minax, maxax-4*dtext, str5)
 text(minax, maxax-5*dtext, str6)
 text(minax+16*dtext, minax+9*dtext, str7)
 text(minax+16*dtext, minax+8*dtext, str8)
 text(minax+16*dtext, minax+7*dtext, str9)
 text(minax+16*dtext, minax+6*dtext, str10)
 text(minax+16*dtext, minax+5*dtext, str11)
 text(minax+16*dtext, minax+4*dtext, str12)
 text(minax+16*dtext, minax+3*dtext, str13)
 text(minax+16*dtext, minax+2*dtext, str14)
 text(minax+16*dtext, minax+1*dtext, str15)
end;