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*c…