clear; close all; r = 1; V_v = [1 1.1 1.2 1.3 1.31 1.4]; R = r./(2*(1./V_v).^2-1); es = (R-r)./(R+r); as = (R+r)/2; bs = as.*sqrt(1-es.^2); theta = linspace(0,2*pi,100); for idx=1:length(bs)-1 b = bs(idx); a = as(idx); e = es(idx); x = a*cos(theta); y = b*sin(theta); plot(x-a,y,'DisplayName',sprintf('V/v=%.2g,e=%.2g',V_v(idx),es(idx))); hold on; end b = bs(end); a = as(end); x = a*cos(theta(1:10)); y = b*sin(theta(1:10)); plot(x-a,y,'DisplayName',sprintf('V/v=%.2g,e=%.2g',V_v(end),es(end))); axis equal; legend(); big;