function grav_r0z0linear(r0,z0)

%r = (0.1:0.5:10) ;
%z = (0.1:0.5:19.75) ;
r = (0.1:0.5:31) ;
z = (0.1:0.5:31) ;
k = 2 * pi * 6.67e-11;
drho = -2500 ;

input = load('gravity.dat');
x  = input(:,1);
dg = input(:,2);
n  = length(x) ;

dgp = zeros(n,1);

for i = 1:n
   dgp(i) = k * drho * r0^2 * z0/(x(i)^2 + z0^2);
end;
for i = 1:n
   G(i,1) = 2*k*drho*r0*z0/(x(i)^2 + z0^2);
   G(i,2) = k*drho*r0^2*(x(i)^2 - z0^2)/(x(i)^2 + z0^2)^2;
end

m = ones(2,1) ;
d = dg - dgp ;

for i = 1:length(r)
    for j = 1:length(z)
        m(1) = r(i) - r0 ;
        m(2) = z(j) - z0 ;
         E(j,i) = (d - G*m)'*(d - G*m) ;
    end
end

v = logspace(-14,14,150);
contour(r,z,E,v); 
title('Contours of Log(E)');  axis ij
xlabel('Radius r'); ylabel('Depth z');