Forward Gravity

grav2d.m

close all
clear all
clc

Gg=6.67e-11; %G
faktor=1; %faktor pengali nilai 				
z=10; %kedalaman model
x=20; %lebar model 
dx=1; %jarak spasial x
dz=1; %jarak spasial z
gridz=z/dz;
gridx=x/dx;
xi=1*dx:dx:x*dx;
zi=1*dx:dz:z*dz;

dens=zeros(z,x);
%dens(gridz/2-2:gridz/2-2,gridx/2+2:gridx/2+2)=10;

%BENTUK MODEL
aa=1;
% dens=dlmread('dens.txt');
for i=1:x
    for j=1:z
        dens(j,i)=9.0;
    end
end

for i=x/2:x
    for j=2:5
        dens(j,i)=1.0;
    end
end

for i=6:13
    for j=1:aa
        dens(j,i)=3.0;
    end
    aa=aa+1;
end

%BENTUK MODEL


for i=1:length(xi)
    for ii=1:length(xi)
    for jj=1:length(zi)
        A=xi(i)-xi(ii)+dx/2;
        B=xi(i)-(xi(ii)+dx/2);
        C=zi(1)-zi(jj)+dz/2;
        D=zi(1)-(zi(jj)+dz/2);
        a(jj,ii)=Gg*(A*log((A^2+D^2)/(A^2+C^2))-B*log((B^2+D^2)/(B^2+C^2))+2*D*(atan(A/D)-atan(B/D))-2*C*(atan(A/C)-atan(B/C)));
    end
    end
    tot=0;
    for jm=1:length(zi)
        for im=1:length(xi)
        tot=tot+a(jm,im)*dens(jm,im);
        end
    end
    anom(i)=tot;
end
anom=anom*1e5;
subplot(211)
plot(xi,anom)
ylabel('anomali (mgal)')
title('Plot Data Anomali Gravitasi')
subplot(212)
imagesc(xi,zi,dens);
%leocd91@gmail.com
xlabel('jarak (m)');
ylabel('kedalaman (m)');
str = sprintf('kontras densitas = %1.1f g/cm^3',max(max(dens)));
title(str)