function [X,Y,T]=Surff3(C,GEO,mu,wt,z)
a=GEO(1);
b=GEO(2);
c=sqrt(a^2-b^2);
e=sqrt(1-(b/a)^2);
p=a*(1-e^2);

rf=[0:1000:1.5*GEO(1)];
ang=[0:1e-2:pi];


for i=1:length(rf)
    for j=1:length(ang)
%         r=a*(1-e*cos(ang(j)));
%         %CALCULO DE THETA
%         if ang(j)<=pi
%             theta=acos((cos(ang(j))-e)/(1-e*cos(ang(j))));
%         elseif ang(j)>pi
%             theta=2*pi-acos((cos(ang(j))-e)/(1-e*cos(ang(j))));
%         end
%         %
        r=a*(1-e^2)/(1+e*cos(ang(j)));
        rx=r*cos(ang(j))+c;
        ry=r*sin(ang(j));
        R=sqrt((rx)^2+(ry)^2);
        %CALCULO DE ANG
            if ry>0
                ANG=acos(rx/R);
            else
                ANG=2*pi-acos(rx/R);
            end
        %
        X(i,j)=(0+rf(i))*cos(ANG);
        Y(i,j)=(0+rf(i))*sin(ANG);
    end
end




for i=1:length(rf)
    for j=1:length(ang)
        T(i,j)=(1/2*wt^2*(X(i,j)^2+Y(i,j)^2)+pot(X(i,j),Y(i,j),z,GEO,mu)-C);
% T(i,j)=X(i,j)+Y(i,j);
    end
end
% figure(20)
% mesh(X,Y,T)
figure(10)
% t = linspace(0,2*pi);plot(GEO(1)*cos(t),GEO(2)*sin(t),'LineWidth',2.5,'Color',[0.1,0.1,1])
hold on
T=real(T);
contour3(X,Y,z+T,[z z],'LineWidth',2.2);
grid on
axis([-2*GEO(1) 2*GEO(1) -2*GEO(1) 2*GEO(1) 0 z+1]*1.1)
str=sprintf('Curvas de velocidad cero, C=%g',C);
title(str);