function phi = ellipse_phi (X, a, b, phi, myeps);%ELLIPSE_PHI    Compute nearest points to ellipse%% phi = ellipse_phi (X, a, b, phi{}, myeps{sqrt(myeps)});% compute angles for nearest points on ellipse to given points X% % X: given points <X(i,1),X(i,2)>% a, b: ellipse parameters with a-axis in x-coordinate% phi: starting values for iteration% myeps: convergence limit%% phi: phi(i) angle of point on ellipse (in parametric form)%      nearest to point <X(i,1), X(i,2)>  if (nargin < 5), myeps = sqrt(eps); end;  if (nargin < 4), phi = []; end;  if (phi == []),  phi = angle (X(:,1)/a + sqrt(-1)*X(:,2)/b); end;    m = size(X,1);  for i = 1:m,    par = [b^2 - a^2; a*X(i,1); b*X(i,2)];    ni  = phi(i);    step = 0;    while (1),      c = cos(ni);      s = sin(ni);      dni = - (par(1)*s*c + par(2)*s - par(3)*c) / ...              (par(1)*(c^2 - s^2) + par(2)*c + par(3)*s);      while ((X(i,1) - a*c)^2 + (X(i,2) - b*s)^2) < ...            ((X(i,1) - a*cos(ni+dni))^2 + (X(i,2) - b*sin(ni+dni))^2),        dni = -0.2*dni;      end      ni = ni + dni;      if (abs(dni) <= sqrt(eps)) break; end;      step = step + 1;      if (step > 40),         disp ('warning: no convergence');        str = sprintf ('angle, diff-angle = %g, %g', ni, dni);        disp (str);        break;       end    end    phi(i) = ni;  end % for all pointsend % ellipse_phi