 %%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;