kepler's equation not working. it solves the Kepler equation for E given e and M (your function should accept M and e as inputs, and return E). The actual solution should be performed using the MATLAB built in function fzero() to solve the following

it solves the Kepler equation for E given e and M (your function should accept M and e as inputs, and return E). The actual solution should be performed using the MATLAB built in function fzero() to solve the following form of Kepler's equation:

0 = E - e*sin[E] – M

function E = kepler(M, e)
f=@(E) E-e*sin(E)-M;
E=fzero(f,0);
%E = eccentric anomaly measured from perihelion about
%the center of the elliptical orbit
%M = mean anomaly = 2p t/P
M=(0:0.01:6);
N=length(M);

%e = eccentricity of the orbit
e=[0 0.25 0.5 0.75 1];

n=length(e);
for x=1:N
     for y=1:n
         E(x,y)=kepler(e(y),M(x));
     end
end

figure, plot(M,E), grid on,

I can't seem to make it work. What am I doing wrong? Any help is highly appreciated

[Rick Rosson]

Please format your code.

Answer 1

[James Tursa]

Issues:

1) Separate your code into two files, your kepler function file and your script file that will call the kepler function file.

2) The e and M argument order in your kepler call are backwards from what the function is expecting, so you need to fix that as well.

3) Learn to use spaces in your code so that it is more readable.

E.g., here is the first file kepler.m

% File kepler.m solves Kepler's equation, E and M are in radians
function E = kepler(M, e)
f = @(E) E - e * sin(E) - M;
E = fzero(f,0);  % <-- I would use M as the initial guess instead of 0
end

And here is the second file kepler_script.m

% E = eccentric anomaly measured from perihelion about
% the center of the elliptical orbit

% M = mean anomaly = 2p t/P
M = (0:0.01:6);
N = length(M);

% e = eccentricity of the orbit
e = [0 0.25 0.5 0.75 1];
n = length(e);

for x=1:N
     for y=1:n
%         E(x,y)=kepler(e(y),M(x));
         E(x,y) = kepler(M(x),e(y));  % <-- Fixed the order of the inputs
     end
end

figure, plot(M,E), grid on,