1n := 4; 2 3on rational, rat; 4off allfac; 5 6array p(n/2+2); 7 8harmonic u,v,w,x,y,z; 9 10weight e=1, b=1, d=1, a=1; 11 12%% Step1: Solve Kepler equation 13bige := fourier 0; 14for k:=1:n do << 15 wtlevel k; 16 bige:=fourier e * hsub(fourier(sin u), u, u, bige, k); 17>>; 18write "Kepler Eqn solution:", bige$ 19 20%% Ensure we do not calculate things of too high an order 21wtlevel n; 22 23%% Step 2: Calculate r/a in terms of e and l 24dd:=-e*e; hh:=3/2; j:=1; cc := 1; 25for i:=1:n/2 do << 26 j:=i*j; hh:=hh-1; cc:=cc+hh*(dd^i)/j 27>>; 28 29bb:=hsub(fourier(1-e*cos u), u, u, bige, n); 30aa:=fourier 1+hdiff(bige,u); ff:=hint(aa*aa*fourier cc,u); 31 32 33%% Step 3: a/r and f 34uu := hsub(bb,u,v); uu:=hsub(uu,e,b); 35vv := hsub(aa,u,v); vv:=hsub(vv,e,b); 36ww := hsub(ff,u,v); ww:=hsub(ww,e,b); 37 38%% Step 4: Substitute f and f' into S 39yy:=ff-ww; zz:=ff+ww; 40xx:=hsub(fourier((1-d*d)*cos(u)),u,u-v+w-x-y+z,yy,n)+ 41 hsub(fourier(d*d*cos(v)),v,u+v+w+x+y-z,zz,n); 42 43%% Step 5: Calculate R 44zz:=bb*vv; yy:=zz*zz*vv; 45 46on fourier; 47 48p(0):= fourier 1; p(1) := xx; 49for i := 2:n/2+2 do << 50 wtlevel n+4-2i; 51 p(i) := fourier ((2*i-1)/i)*xx*p(i-1) - fourier ((i-1)/i)*p(i-2); 52>>; 53 54wtlevel n; 55for i:=n/2+2 step -1 until 3 do p(n/2+2):=fourier(a*a)*zz*p(n/2+2)+p(i-1); 56 57yy*p(n/2+2); 58 59showtime; 60 61end; 62