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

% [sd1,sd2] = ellipse(NBe,delNBe,NAl,delNAl,plotSelect,plotString)
%
% Produces uncertainty ellipses for "Banana" cosmogenic-isotope plots.
%
% Greg Balco
%
% UW Cosmogenic Isotope Lab
% August 2000
%
% Translation of John's ellipse-drawing Mathematica code (which is dated summer 1997 ). 
% 
% Basically a plotting function. Plots one sample and returns handles. 
% Plots 1-sigma and 2-sigma ellipses, returns handles to 1- and 2- sigma 
%   ellipses (actually 68% and 95% confidence intervals) 
%
% Arguments are NBe, delNBe, NAl, delNAl, plotSelect, plotString
%
% plotSelect argument: optional -- 1, only 1-sigma ellipse; 2, both 1- and 2-sigma
% plotString argument: optional string argument giving desired plot linetype 
%
% e.g. [h1 h2] = ellipse(1e5,0.1e5,6.1e5,0.61e5,1,'r');
%
% Note: remember to normalize N's to sea level/high latitude before feeding them to this function. 
%
% Here's what it's doing: 
%
% 1: forms PDF calculation -- evaluates 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.
% This gives so-called "multinormal" PDF, I think.             
%
% Then, contours the resulting PDF at 68% and 95%. Not mathematically pretty. 	
%

% Check args

if nargin >= 5 
	if plotSelect ~= 1 & plotSelect ~= 2;
		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));
		
% So we now know the probability (multinormal distribution?) 
% for each cell in the grid. 

% 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;

sd1 = plot(x1,y1,plotString);

if plotSelect == 2;
	
	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);
end;










		

						
