function y = P(A,B,C,Z,TE,U,x1,y1,al)
%
ga=atan2(y1,x1); s=sqrt(x1^2+y1^2);
te=TE-al;
t=U-s*cos(ga-TE);

[NK,M]=size(TE);
z=Z(:)*ones(1,M);

ttt=1-(z/C).^2; ttt=sqrt(ttt.*(ttt > 0));
a=A*ttt; b=B*ttt;

%
at=(a.*cos(te)).^2+(b.*sin(te)).^2;

t2=t.^2;
ttt=at-t2; ttt=sqrt(ttt.*(ttt > 0));
y=(2*(a.*b).*ttt)./(at+1e20*(at == 0));