/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 0.8.5 ] */

/* [wxMaxima: section start ]
Refutation of Einstein General Relativity orbits
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
define radius function r(theta) of a precessing ellipse
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
r = alpha/(1+epsilon*cos(x*theta));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
multiply by (epsilon*cos(x*theta)+1)
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
%*(epsilon*cos(x*theta)+1);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
divide by r
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
%/r;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
subtract 1
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
%-1;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
divide by epsilon
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
%/epsilon;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
simplify
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
expand(ratsimp(%));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
transform cos to sin by using theorem sin(x)^2 + cos(x)^2 = 1
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
sin(x*theta)^2 = 1 - cos(x*theta)^2;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
insert result above for cos(x*theta)
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
sin(x*theta)^2 = 1 - (alpha/(epsilon*r)-1/epsilon)^2;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
work out square of sum and take common denominator
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
radcan(%);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
simplify
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
expand(%);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
This expression has to be eqated with the result of EGR
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
E: (2*alpha)/(epsilon^2*r)-alpha^2/(epsilon^2*r^2)-1/epsilon^2+1 =
(alpha/(x*epsilon))^2*(1/b^2-1/a^2+r0/a/r-1/r^2+r0/r^3);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
multiply by epsilon^2*x^2
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
E*epsilon^2*x^2;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
It is seen that the left hand side is a polynom of maximum degree 1/r^2 
while the right hand side is one of maximum degree 1/r^3. Both sides can 
never be equal for a certain range of r. However the variable r must be
able to vary because it describes an ellipse - reduction ad absurdum.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Plot of precessing ellipses
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
r: alpha/(1+epsilon*cos(x*theta));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
r1: ev(r, [alpha=1, epsilon=0.2, x=1]);
r2: ev(r, [alpha=1, epsilon=0.2, x=1.1]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: section start ]
Orbital Plots
   [wxMaxima: section end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(user_preamble = "set grid polar; set size square;",
             nticks        = 1000,
             line_width    = 2,
             xrange        = [-1.5,1.5],
             yrange        = [-1.3,1.7],
             title         = "Precessing Ellipse, 2 orbits",
key           = "r(theta), x=1.0",
color         = red,
polar(r1, theta, 0, 4.*%pi),
key           = "r(theta), x=1.1",
color         = blue,
polar(r2, theta, 0, 4.*%pi))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(user_preamble = "set grid polar; set size square;",
             nticks        = 1000,
             line_width    = 2,
             xrange        = [-1.5,1.5],
             yrange        = [-1.3,1.7],
             title         = "Precessing Ellipse, 12 orbits",
key           = "r(theta), x=1.0",
color         = red,
polar(r1, theta, 0, 24.*%pi),
key           = "r(theta), x=1.1",
color         = blue,
polar(r2, theta, 0, 24.*%pi))$
/* [wxMaxima: input   end   ] */

/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$