%%Función que calcula el potencial gravitatorio de un asteroide elipsoidal
%%dada su geometría y su densidad
function U=pot(x,y,z,GEO,mu)
%GEO=[a b c]
%mu=constante gravitacional (volumen*densidad*G)
a=GEO(1);b=GEO(2);c=GEO(3);
%tic
Rm=(a+b+c)/2;
ti=[a-b;b-c];
A=[(b/2*(Rm/a)^2-a/2*(Rm/b)^2) -(3*b*(Rm/a)^2+3*a*(Rm/b)^2);(c/2*(Rm/b)^2+b*(Rm/c)^2) (3*c*(Rm/b)^2)];
Cost=A\ti;
%toc
J2=Cost(1);
J22=Cost(2);
r=sqrt(x^2+y^2+z^2);
phi=asin(z/r);
if y>=0
lambda=acos(x/r/cos(phi));
else
lambda=2*pi-acos(x/r/cos(phi));
end
U=mu/r*(1+J2/2*(Rm/r)^2*(1-3*sin(phi)^2)-3*J22*(Rm/r)^2*cos(phi)^2*cos(2*lambda));
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
