function ageprob = probdens(samplename, ae, startage, endage, interval, histbin, minerr, fontsize)

% function computes the probability distribution function for
% a set of geochronological ages and their associated uncertainty
% Probability function from Sircombe (2000)
%
% SYNTAX: ageprob = probdens(samplename, ae, startage, endage, interval,
% histbin, minerr, fontsize)
%
% samplename: Sample name as a string ex-'AM04-01'
% ae: 2 x j matrix of  ages and their uncertainties in Ma (j- # of analyses)
% startage: minimum age in Ma to compute summed probability
% endage: maximum age in Ma to compute summed probability
% interval: step from min to max age using this interval in Ma
% histbin: binsize for histogram in Ma
% minerr: minimum error estimate to be associated with a single kernel set (35 Ma)
% fontsize: fontsize for figures in points (optional)
%
% 1/07- function by R. Mapes, russ_mapes(a)hotmail.com


% make sure correct number of arguments have been input
if nargin < 7
    error('Not enough input aguments')
end

if nargin == 7
    fontsize = 10; % default size
else
    % use input fontsize
end

% figure out how many steps to use between the start age and the end
steps = 1+((endage-startage)/interval);

% get rid of NaN values in ae matrix
ae = removeNaN(ae)

% split age-error matrix into two columns
a = ae(:,1);
e = ae(:,2);

% create matrix for individual probabilities
ageprobmat=zeros(length(a),steps);

% fill the probability matrix for each age/error at each time
for i = 1:length(a) % # of analyses
    %start the time clock
    time=startage;
    
    if e(i) < minerr % this reassigns the uncertainty estimate if it is below a user defined level
        e(i) = minerr;
    end
    
    for j=1:steps % # of intervals
	% Function for the probability of an analysis with normally distributed error
	% at the given time step from Sircombe (2000).
        ageprobmat(i,j)=(1/(e(i)*(2*pi)))*exp((-(time-a(i))^2)/((2*e(i))^2));

	% Step forward to the next time interval
        time = time+interval; % step time forward one step
    end
end

% Create a matrix of the ages tested to tack onto the probability matrix
time = startage;
for i = 1:steps
    times(i)=time;
    time=time+interval;
end

% create a 2 x j matrix of ages and summed probability 
ageprob = [times' sum(ageprobmat,1)'];

%sumageprob = sum(ageprob(:,2)); %figure out what the total area under the prob curve is
ageprob = [ageprob(:,1) (ageprob(:,2)/length(a))]; % divide by the total so total probability sums to 1

% plot raw ages as histogram
set(gcf,'DefaultAxesColorOrder',[0 0 1; 1 0 0])
set(gca,'FontSize',fontsize)

hist(a,startage:histbin:endage)
% set the limits on the x-axis
xlim([startage endage])
hold on
% overlay probability
[AX,H1,H2] = plotyy(NaN,NaN,ageprob(:,1),ageprob(:,2));

% Sets the line style for the secondary axis
%set(colordef, AX(2), 'color','r')
set(H2,'LineWidth', 2)

% label the shared x-axis
xlabel('Age (Ma)')

% Label y-axis
ylabel(['Frequency (Bin = ',num2str(histbin),' Ma)'])
set(get(AX(2),'Ylabel'),'String','Probability','fontsize',fontsize)
set(AX(2),'fontsize',fontsize)

% Give it a title
set(gca, 'FontWeight', 'bold')
title([samplename, ' (n = ', num2str(length(a)), ')'])
set(gca, 'FontWeight', 'normal')

hold off

%%%%%%%%% Uncomment to save Illustrator file %%%%%%%%%
% If saveas command is uncommented the m-file will save an Adobe Illustrator file of the plot 
%saveas(gcf, [samplename, 'prob.ai'],'ai')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function newmat = removeNaN(oldmat)

% removes NaN values from a matrix
%
% syntax:
% newmat = removeNaN(oldmat)

newmat = oldmat;

% figure out which rows contain NaN and remove them
i = 1;
while i < 1 + length(newmat)
    if isnan(newmat(i,1)) == 1
        newmat = [newmat(1:i-1,:)
                  newmat(i+1:end,:)];
        i = 0;
    end
    i = i + 1;
end