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);
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
