function th_max = orbit(e)

N_e = 1000;
th_e = linspace(0,pi,N_e);
result = zeros(1,N_e);

N_x = 4000;

for k = 1:N_e
    % Generate the th_x grid
    th_x = linspace(pi,2*pi,N_x);
    
    wt    = th_e(k) - 2*e*sin(th_e(k));
    error = abs(th_x - 2*e*sin(th_x) - wt - pi);
    
    min_index = 1;
    for l = 2:N_x
        if( error(l) < error(min_index) )
            min_index = l;
        end
    end
    result(k) = th_x(min_index);
end

delta_theta = result-th_e; % result holds th_x;

ratio = sin(delta_theta)./((1+e*cos(result))./(1+e*cos(th_e)) - cos(delta_theta));
alpha = atan(ratio);

plot(th_e,alpha);
grid on;
xlabel('\epsilon, eccentricity');
ylabel('\alpha, <XES angle');
max_index = 1;
for k = 2:N_e
    if( alpha(k) > alpha(max_index) )
        max_index = k;
    end
end

th_max = th_e(max_index);