% computes the vertical gravity anomaly from a buried cylinder of radius R,
% located at depth Z, with axis along the y-direction, x>0 to the right and
% z>0 downward
% = a benchmark for gravity routine gpoly.m
% J.A. Olive, MIT/WHOI Joint Program
% Dec. 2010
clear all

% CYLINDER
R=30; % cylinder radius (m)
Z=100; % cylinder depth (m) 
Drho=1000; % density contrast between cylinder and surrounding (kg/m^3)

% OBSERVATION PROFILE
no=100; % number of discrete points of observation
xo=linspace(-6*Z,6*Z,no);
zo=zeros(1,no);



N=100; % number of discrete points making up cylinder section
i=sqrt(-1);
idx=1:1:N;
cplx=R*exp(2*i*idx*pi./N);
xs=real(cplx);
zs=imag(cplx)+Z;

gz=zeros(1,no);
for u=1:no % sweep observation line
    gz(u)=gpoly(xo(u)/1000,zo(u)/1000,xs/1000,zs/1000,N,Drho);
end
  
% analytical solution for a buried cylinder
gamma=6.673e-11;
an=1e5*2*pi*(R)^2*Drho*gamma*Z.*(xo.^2+Z^2).^(-1);



subplot(211)
plot(xo,an,'r')
hold on
plot(xo,gz,'.')
legend('analytic','calculated')
ylabel('\Delta g_z (mgals)','fontsize',15)
xlabel('x (m)','fontsize',15)
grid on
subplot(212)
plot(xs,zs,'.')
axis equal
hold on
plot(xo,zo,'k.')
set(gca,'YDir','Reverse')
xlabel('x (m)','fontsize',15)
ylabel('z (m)','fontsize',15)
grid on