function [out] = camelplot(ages,errors)

% camelplot.m
%
% Syntax: [handles] = camelplot(ages,errors)
%
% Plots camel diagrams, i.e. synthetic probability ideograms, for 
% a set of data with errors, for example a collection of radiometric dates. 
% 
% for a discussion of the basic camel diagram concept see:
%
% Lowell, T. V., 1995, The application of radiocarbon age estimates to
% the dating of glacial sequences: An example from the Miami sublobe, Ohio,
% U.S.A.: Quaternary Science Reviews, v. 14, p. 85-99.
%
% Input arguments are vectors containing the values and errors, respectively, 
% of data values that are assumed to represent gaussian distributions around 
% the measured values. The error in each measurement is taken to be the S.D. 
% of its gaussian distribution. 
% 
% Returns a vector of handles to graphic objects: first the handles to the 
% gaussians corresponding to individual data points, then the handle to the 
% summary curve. 
%
% Greg Balco -- UW Cosmogenic Isotope Lab -- Feb. 2001
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License, version 2,
% as published by the Free Software Foundation (www.fsf.org).


numsamples = length(ages);

minAge = min(ages) - (4*errors(find(ages == min(ages))));
maxAge = max(ages) + (4*errors(find(ages == max(ages))));
ageRange = maxAge - minAge;


x = minAge:(ageRange/1000):maxAge;
totalPdf = zeros(size(x));

for a = 1:numsamples;
	
	mu = ages(a); sigma = errors(a);
	xn = (x - mu) ./ sigma;
	y = exp(-0.5 * xn .^2) ./ (sqrt(2*pi) .* sigma);
	eval(['pdf' num2str(a) ' = y;'])
	
	%eval(['pdf' num2str(a) ' = normpdf(x,ages(a),errors(a));']);
	eval(['totalPdf = totalPdf + pdf' num2str(a) ';']);
end;

hold on;
grid on;
b = [];

% plot individual curves, normalised to 1/n
for a = 1:numsamples;
	eval(['h' int2str(a) ' = plot(x,pdf' num2str(a) '/numsamples,''r'');']);
	eval(['b = [b h' int2str(a) '];']);
end;

% total curve, normalised to 1
hx = plot(x,totalPdf/numsamples,'k');

out = [b hx];