function [E,number_of_terms]=Kepler(e,M)

%%Kepler.m will take in the values e and M in the expression M=E-esinE and
%%solve for E exactly by using an expansion in Bessel functions of the
%%first kind of order m

sum_of_bessels=0;           %initialize sum to zero                 
i=1;                        %index of summation set to 1


while abs((1/i)*besselj(i,i*e)*sin(i*M)) > eps       %Term in expansion is negligible
    sum_of_bessels=sum_of_bessels+(1/i)*besselj(i,i*e)*sin(i*M);
    number_of_terms=i;
    i=i+1;
end

E=M+2*sum_of_bessels;