# “Power series” in Cos[n*x] for X coordinate of planetary ellipse

Consider the functions below:


ellipseAreaFromFocus[a_,b_,t_] = a*b*t/2 - Sqrt[a^2-b^2]*b*Sin[t]/2

(* NOTE: origin is midpoint of two focii, x = semimajor axis *)

ellipseMA2XY[a_,b_,ma_] := Module[{t},
t = x /. FindRoot[ellipseAreaFromFocus[a,b,x]==a*b*ma/2,{x,0}];
{a*Cos[t], b*Sin[t]}
];


If we sample the first coordinate of ellipseMA2XY[1.022,1,t] at several points, we can approximate it with a linear combination of Cos[n*x]:


t = Table[{t,ellipseMA2XY[1.022,1,t][[1]]},{t,0,2*Pi,.001}];

f[x_] = Fit[t,Table[Cos[n*x],{n,0,6}],x];

-0.10545614758232982 + 1.0057258170784746*Cos[x] +
0.10248576787806324*Cos[2*x] + 0.01567938482407289*Cos[3*x] +
0.0028440630822004482*Cos[4*x] + 0.0005668729202178632*Cos[5*x] +
0.00011997552131595974*Cos[6*x]


My question: assuming I carried this series to infinity, what's the formula for the coefficient of Cos[n*x] in terms of 1.022 (or any other value I put in there).

The eccentricity of this ellipse is Sqrt[1-1/1.022^2] or 0.206372, and the constant term appears to be about -1/2 of that (it might become exactly -1/2 in the infinite expansion), but the others leave me semi-stumped.

It is known that there's no finite expression for ellipseMA2XY[1.022,1,t][[1]], but I'm convinced that an infinite expansion in Cos[n*x] exists (and is probably well-known)

(I realize I could ask this in math.stackexchange.com, but since I'm using Mathematica, I figured I'd try here first)

• Could maybe get a good approximation using a discrete Fourier transform method. – Daniel Lichtblau Nov 29 '14 at 18:23
• What you're asking for is the Fourier series of ellipseMA2XY. Try FourierSeries. – Rahul Nov 30 '14 at 7:42

I got nowhere useful with the DCT so I went instead with straight Fourier. The idea is take the discrete Fourier transform of the points, use the first as the dc component, sum the rest with the reversal thereof to get cosines from complex exponentials, and only use the first half since we don't want to confuse negative frequencies with high frequencies.

fft = Fourier[t[[All, 2]]]/Sqrt[Length[t]];
dc = First[fft];
fftrest = Rest[fft];
fftcos = Re[fftrest + Conjugate[Reverse[fftrest]]];
len = Floor[Length[fftcos]/2]


Let's have a look at some of these.

dc
fftcos[[1 ;; 10]]

(* -0.1053099823173539 + 0. I

{1.005831877937645, 0.102413921396564, 0.01562846836386967,
0.002813229515042125, 0.0005474867224487756, 0.0001069486252433002,
0.00001716571734559302, -9.589586401412791*10^-7,
-4.036307493988713*10^-6, -4.029977347108143*10^-6} *)


They look to be reasonably close to what Fit gave, and that is good news of course.

Plot[dc + Take[fftcos, len].Table[Cos[n*x], {n, len}], {x, 0, 2*Pi}]