%%Función que calcula el potencial gravitatorio de un asteroide elipsoidal
%%dada su geometría y su densidad
function [GGU U]=GGpot(x,y,z,GEO,mu)
%%%GEO=[a b c]%
%%%mu=constante gravitacional (volumen*densidad*G)
a=GEO(1);b=GEO(2);c=GEO(3);
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;
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));
Ur=-r^-1*(U)+mu/r*(-J2*(Rm^2/r^3)*(1-3*sin(phi)^2)+6*J22*(Rm^2/r^3)*cos(phi)^2*cos(2*lambda));
Ul=mu/r*(6*J22*(Rm/r)^2*cos(phi)^2*sin(2*lambda));
Up=mu/r*(-J2*(Rm/r)^2*3*sin(phi)*cos(phi)+6*J22*(Rm/r)^2*cos(phi)*sin(phi)*cos(2*lambda));
B=[cos(phi)*cos(lambda) cos(phi)*sin(lambda) sin(phi);-r*cos(phi)*sin(lambda) r*cos(phi)*cos(lambda) 0;-r*sin(phi)*cos(lambda) -r*sin(phi)*sin(lambda) r*cos(phi)];
GU=B\[Ur Ul Up]';
% er=[cos(phi)*cos(lambda) cos(phi)*sin(lambda) sin(phi)]';
% ep=[-sin(phi)*cos(lambda) -sin(phi)*sin(lambda) cos(phi)]';
% el=[-cos(phi)*sin(lambda) cos(phi)*cos(lambda) 0]';
%
% GU=Ur*er+1/r*Up*ep+1/r*cos(phi)*Ul*el;
%Uxx
inc=1e-10;
x=x+inc;
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));
Ur=-r^-1*(U)+mu/r*(-J2*(Rm^2/r^3)*(1-3*sin(phi)^2)+6*J22*(Rm^2/r^3)*cos(phi)^2*cos(2*lambda));
Ul=mu/r*(6*J22*(Rm/r)^2*cos(phi)^2*sin(2*lambda));
Up=mu/r*(-J2*(Rm/r)^2*3*sin(phi)*cos(phi)+6*J22*(Rm/r)^2*cos(phi)*sin(phi)*cos(2*lambda));
B=[cos(phi)*cos(lambda) cos(phi)*sin(lambda) sin(phi);-r*cos(phi)*sin(lambda) r*cos(phi)*cos(lambda) 0;-r*sin(phi)*cos(lambda) -r*sin(phi)*sin(lambda) r*cos(phi)];
GUinc=B\[Ur Ul Up]';
GGU(1)=(GUinc(1)-GU(1))/inc;
%Uyy
y=y+inc;
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));
Ur=-r^-1*(U)+mu/r*(-J2*(Rm^2/r^3)*(1-3*sin(phi)^2)+6*J22*(Rm^2/r^3)*cos(phi)^2*cos(2*lambda));
Ul=mu/r*(6*J22*(Rm/r)^2*cos(phi)^2*sin(2*lambda));
Up=mu/r*(-J2*(Rm/r)^2*3*sin(phi)*cos(phi)+6*J22*(Rm/r)^2*cos(phi)*sin(phi)*cos(2*lambda));
B=[cos(phi)*cos(lambda) cos(phi)*sin(lambda) sin(phi);-r*cos(phi)*sin(lambda) r*cos(phi)*cos(lambda) 0;-r*sin(phi)*cos(lambda) -r*sin(phi)*sin(lambda) r*cos(phi)];
GUinc=B\[Ur Ul Up]';
GGU(2)=(GUinc(2)-GU(2))/inc;
%Uzz
z=z+inc;
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));
Ur=-r^-1*(U)+mu/r*(-J2*(Rm^2/r^3)*(1-3*sin(phi)^2)+6*J22*(Rm^2/r^3)*cos(phi)^2*cos(2*lambda));
Ul=mu/r*(6*J22*(Rm/r)^2*cos(phi)^2*sin(2*lambda));
Up=mu/r*(-J2*(Rm/r)^2*3*sin(phi)*cos(phi)+6*J22*(Rm/r)^2*cos(phi)*sin(phi)*cos(2*lambda));
B=[cos(phi)*cos(lambda) cos(phi)*sin(lambda) sin(phi);-r*cos(phi)*sin(lambda) r*cos(phi)*cos(lambda) 0;-r*sin(phi)*cos(lambda) -r*sin(phi)*sin(lambda) r*cos(phi)];
GUinc=B\[Ur Ul Up]';
GGU(3)=(GUinc(3)-GU(3))/inc;
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
