function [x,y] = ellipse(a,e,perh,npts);
% [x,y] = ellip(a,e,perh,npts)
% generates an elliptical orbit with eccentricity = e, perihelion = perh
% npts points equally spaced in time
% uses iteration, and formula from Landau & Lifshitz Mechanics pg 38
% r = a(1-e*cos(xi));	
% t = sqrt(a^3/GM)*(xi - e*sin(xi));
% x = a*(cos(xi)-e);
% y = a*sqrt(1-e^2)*sin(xi);

% constants
	GM = 4*pi^2;
	sqrtGM = 2*pi;		% GM = 4*pi^2 in units with years, AU. 
	a3 = a^3;
% period T
	period = sqrt(a3);	% T = 2*pi*a^(3/2)/sqrt(GM);
	t = linspace(0,period,npts)';
% solve for xi
	fac = sqrt(a3/GM);
	xi = linspace(0,2*pi,npts)';	% starting value
		xiold = 0;	ni = 0;
	for k = 1:20;	%number of iterations.  Most results not very sensitive.
		xiold = xi;
		xi = e*sin(xi) + t/fac;		
	end
	x1 = a*(cos(xi)-e);
	y1 = a*sqrt(1-e^2)*sin(xi);
% rotate to perihelion perh;
	x = x1*cos(perh) - y1*sin(perh);
	y = y1*cos(perh) + x1*sin(perh);