  function [xM,yM] = Mars(t)
% (xM,yM) is the location of Mars at time t (days).
P = 206.6; A = 249.22; phi = 336*pi/180; Y = 686.93;
x = (P-A)/2 + ((P+A)/2)*cos((2*pi/Y)*t);
y = sqrt(P*A)*sin((2*pi/Y)*t);
c = cos(phi); s = sin(phi);
xM = c*x - s*y;
yM = s*x + c*y;

alpha1 = 2*pi/Y;
alpha2 = -alpha1;
r1 = ((P+A)/2 + sqrt(A*P))/2;
r2 = ((P+A)/2 - sqrt(A*P))/2;
a = (P-A)/2;


xM1 = a*c + r1*cos(alpha1*t+phi) + r2*cos(alpha2*t+phi);
yM1 = a*s + r1*sin(alpha1*t+phi) + r2*sin(alpha2*t+phi);


