%% Interpolate Asteroid Data
% In a paper by Gauss, he describes an approach to estimating the orbit of the Pallas asteroid. He starts with the following twelve 2-D positional data points x and y.

x = 0:30:330;
y = [408 89 -66 10 338 807 1238 1511 1583 1462 1183 804];
plot(x,y,'ro')
xlim([0 360])

% Gauss models the asteroid's orbit with a trigonometric polynomial of the following form.
% \begin{array}{rcl}
% y = a_{0} &+& a_{1}\cos(2\pi(x/360))+b_{1}\sin(2\pi(x/360))\\
% &&{}a_{2}\cos(2\pi (2x/360))+b_{2}\sin(2\pi(2x/360))\\
% &&{}\cdots \\
% &&{}a_{5}\cos(2\pi(5x/360))+b_{5}\sin(2\pi(5x/360))\\
% &&{}a_{6}\cos(2\pi(6x/360))
% \end{array}
%
% Use fft to compute the coefficients of the polynomial.

m = length(y);
n = floor((m+1)/2);
z = fft(y)/m;

a0 = z(1); 
an = 2*real(z(2:n));
a6 = z(n+1);
bn = -2*imag(z(2:n));

%
% Plot the interpolating polynomial over the original data points.

hold on
px = 0:0.01:360;
k = 1:length(an);
py = a0 + an*cos(2*pi*k'*px/360) ...
        + bn*sin(2*pi*k'*px/360) ...
        + a6*cos(2*pi*6*px/360); 

plot(px,py)

% Copyright 2012 The MathWorks, Inc.