function [sd1,sd2] = ellipse(NBe,delNBe,NAl,delNAl,plotSelect,plotString);

% Produces uncertainty ellipses for [Be-10] - [Al-26]/[Be-10] plots, that
% is, the so-called Lal-Klein-Nishiizumi two-isotope diagram. Sometimes 
% known as the 'banana diagram' by the uncultured. 
%
% Syntax:
%   [h1,h2] = ellipse(NBe,delNBe,NAl,delNAl,2,plotString)
%   [h1] = ellipse(NBe,delNBe,NAl,delNAl,1,plotString)
%   [x,y] = ellipse(NBe,delNBe,NAl,delNAl,0,plotString)
%
% Arguments are: NBe, delNBe, NAl, delNAl, plotSelect, plotString
%
% The N and delN arguments are Al-26 and Be-10 concentrations and 
% uncertainties therein. 
%   Note: in most cases one will want to plot ellipses with nuclide 
%   concentrations normalized to site production rates: this function 
%   doesn't do the normalization. Normalize nuclide concentrations to match
%   your desired plot axes upstream. 
%
% The argument plotSelect tells the function what, if anything, to plot. 
%   plotSelect = 0 causes function not to plot anything, and to return the x,y 
%       coordinates of the 68% confidence ellipse. In this case the syntax is
%       [x,y] = ellipse(...).
%   plotSelect = 1 plots the 68% confidence ellipse and returns its handle. 
%   plotSelect = 2 plots the 68% and 95% confidence ellipses and returns their
%       handles, in that order.
% 
% The argument plotString is an optional string argument containing a plot 
%   linetype, e.g., [h1 h2] = ellipse(1e5,0.1e5,6.1e5,0.61e5,1,'r');
%
% Here's what actually goes on in the function:
%
% Forms joint PDF by evaluating the expression
%
% x  * exp(-0.5 * ((((y * x) -NAl) / delNAl)^2 + ((x - NBe) / delNBe)^2))
%	
% on a grid +/- 4 SD in NBe and NAl/NBe directions, grid spacing 0.1 SD.
%
% Then, contours the resulting PDF at 68% and 95%.  
%
%
% Written by Greg Balco -- UW Cosmogenic Nuclide Lab
% balcs@u.washington.edu
% Modified from Mathematica code by John Stone dated Summer 1997.  
% First MATLAB version, August 2000
% Revised March 2006.
% Part of the CRONUS-Earth online calculators: 
%      http://hess.ess.washington.edu/math
%
% Copyright 2001-2007, University of Washington
% All rights reserved
% Developed in part with funding from the National Science Foundation.
%
% 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).


% Check args

if nargin >= 5 
	if plotSelect ~= 1 & plotSelect ~= 2 & plotSelect ~= 0;
		plotSelect = 1;
	end;
end;

if nargin < 5
	plotSelect = 1;
end;

% Estimate range and create mesh

R = (NAl/NBe);
delR = ((delNAl / NAl)^2 + (delNBe / NBe)^2)^0.5;

[x,y] = meshgrid((NBe - 4 * delNBe):(0.1 * delNBe):(NBe + 4 * delNBe),...
		(R * (1 - 4 * delR)):(0.1 * R * delR):(R * (1 + 4 * delR)));

% calculate PDF

Prob = x  .* exp(-0.5 .* ((((y .* x) - NAl) ./ delNAl).^2 + ((x - NBe) ./ delNBe).^2));
		
% Now find the 68% probability contour. 

% Normalize to volume = 1

normP = Prob ./ sum(sum(Prob));

% Now we need to figure out cumulative probabilities. 

% multiply by 10000 to achieve manageable number of values, round to get integers:

normP = normP * 10000;
intP = round(normP);

for a = 1:max(max(intP));
	cumprob(a) = sum(intP(find(intP >= a)))/10000;
end;

probs = 1:max(max(intP));

sigma1 = find(abs(cumprob - 0.68) == min(abs(cumprob - 0.68)));

sigma2 = find(abs(cumprob - 0.95) == min(abs(cumprob - 0.95)));

% weed out cases where adjacent probs are same -- rounding error

if length(sigma1) ~= 1;
	sigma1 = min(sigma1);
end;

if length(sigma2) ~= 1;
	sigma2 = min(sigma2);
end;

% Now draw the contours.

cmat = contourc(x(1,:),y(:,1),normP,[sigma1 sigma2]);

cl1 = cmat(2,1);
s = size(cmat);

% Sometimes contourc returns a bunch of contours for one level -- grid size issue?
% This is spurious, so plot only the major one. 

contourStarts = find(cmat(1,:) == sigma1);
contourSizes = cmat(2,contourStarts);
contourToPlot = find(contourSizes == max(contourSizes));

x1 = cmat(1,(contourStarts(contourToPlot)+1):(contourStarts(contourToPlot) + contourSizes(contourToPlot)));
y1 = cmat(2,(contourStarts(contourToPlot)+1):(contourStarts(contourToPlot) + contourSizes(contourToPlot)));

if nargin < 6;
	plotString = 'k';
end;

% return results

if plotSelect == 1;
    sd1 = plot(x1,y1,plotString);

elseif plotSelect == 2;
	
    sd1 = plot(x1,y1,plotString);
    
    hold on;
	
    contourStarts = find(cmat(1,:) == sigma2);
    contourSizes = cmat(2,contourStarts);
    contourToPlot = find(contourSizes == max(contourSizes));

    x2 = cmat(1,(contourStarts(contourToPlot)+1):(contourStarts(contourToPlot) + contourSizes(contourToPlot)));
    y2 = cmat(2,(contourStarts(contourToPlot)+1):(contourStarts(contourToPlot) + contourSizes(contourToPlot)));

    sd2 = plot(x2,y2,plotString);

elseif plotSelect == 0;
    sd1 = x1;
    sd2 = y1;    
end;










		

						
