Mine

function r=r_gen(theta_sel)
   G = 50;
   beta = 3;
   K = 5;
   thetai = [-180:7.2:180];
   thetai(end) = [];
   f = G.*exp(beta.*(cos(deg2rad(theta_sel-thetai))-1)) + K;
   r = zeros(size(f));
   for i=1:50
      r(i) = poirv(1,f(i));
   end