camal.tst 1.3 KB

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