% Draw gravity well
constants
close all
offset = 0.05;
x = ((-2:0.1:2)-offset).*10^9;
y = ((-2:0.1:2)-offset).*10^9;

[X,Y] = meshgrid(x,y);
r = sqrt(X.^2+Y.^2);
z = -G*M_E./r;
surf(X,Y, z) %% Plot grav well.

% zlim([-600e6 0])

rsoi = AU*(M_E/M_S)^(2/5)
hold on
 theta = linspace(0, 2*pi, 1000);
 plot3(rsoi*cos(theta)+offset*10^9, rsoi*sin(theta)+offset*10^9, zeros(1000,1)) % Plot sphere (radius) of influence.
 
 plot3(R_E*cos(theta)+offset*10^9, R_E*sin(theta)+offset*10^9, zeros(1000,1))
axis square
camproj('perspective')