function [U,c] = austretender(M,r,P,d,k)
% Example:
%
% r = 1; m = [0;0]; p = [2.3;2]; d = -[1;1]; refraction = 1;
% [pa,da] = austretender(m, r, p, d, 1)
% t = linspace(0,2*pi,100);
% hold off
% plot(r*sin(t),r*cos(t),'r-'); axis equal; hold on
% plot([p(1),pa(1)],[p(2),pa(2)],'-o')
% [pa,da] = austretender(m,r,p,d,refraction);
% plot([pa(1),pa(1)+da(1)],[pa(2),pa(2)+da(2)],'-cx'), hold off
%

P = P(:); M = M(:);
% normiere d
d = d(:)/norm(d); 
sig = (d(2,:)>=0)-(d(2,:)<0); % Vorzeichen von dy

% berechne Schnittpunkt
dd = r^2 - ( d(2) * (P(1)-M(1)) - d(1) * (P(2)-M(2)) )^2;


if dd >= 0  
  U = P + (sig*sqrt(dd) - d' * (P-M)) * d;
else
  error('Strahl schneidet Sphaere nicht.')
end

% Normale am Schnittpunkt
n = U-M; 
n = n(:)/norm(n);

% neue Austrittsrichtung
sin1 = sig*((n(2) * d(1) - n(1) * d(2)));   % sin(alpha1) = n x d
sin2 = sin1/k;
if sin2<=1
  cos2 = sqrt(1-sin2^2);
  c = -[cos2,sin2;-sin2,cos2]*n;
else
  error('Totalreflexion') 
end