function [X,Y,T]=Surff(C,GEO,mu,wt)
a=GEO(1);
b=GEO(2);
c=sqrt(a^2-b^2);
e=sqrt(1-(b/a)^2);
p=a*(1-e^2);
rf=[-GEO(1):500:15*GEO(1)];
ang=[0:1e-2: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)=(R+rf(i))*cos(ANG);
Y(i,j)=(R+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),0,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
contour(X,Y,T,[0 0],'LineWidth',2.2);
grid on
legend('Asteroide','Curva de T=0')
axis equal
str=sprintf('Curva de velocidad cero en el plano z=0, C=%g',C);
title(str);
xlabel('Eje X (Km)')
ylabel('Eje Y (Km)')
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
