PROGRAM PENENTUAN HIPOSENTER - MATLAB

berikut adalah script m.file matlab untuk program penentuan Hiposenter 3 dimensi.

%PROGRAM PENENTUAN HIPOSENTER

clear all;clc;

%Membuka file data
fid=fopen('wes_dadiII.dat');
data=textscan(fid,'%s %f %f %f','delimited');
fclose(fid);

code=data{1};
x=data{2};
y=data{3};
t=data{4};
N=length(x);
xp=[];yp=[];zp=[];

%Model Awal*
x0 = input ('masukkan longitude awal : ');
y0 = input ('masukkan latitude awal : ');
z0 = input ('masukkan kedalaman awal : ');
vp=4;Erms=1;

%Model Awal**
%x0=25;y0=25;z0=10;vp=4;Erms=0.1;

%MODEL AWAL BISA DIPILIH MENGGUNAKAN Model Awal* atau Model Awal**


%Pemodelan Inversi Menggunakan Iterasi
k=0;
while Erms>0.1
tp=sqrt((x0-x).^2+(y0-y).^2+(z0).^2)/vp;
t0_i=abs(t-tp);t0=mean(t0_i);
tcal=t0+sqrt((x0-x).^2+(y0-y).^2+(z0).^2)/vp;

dt=t-tcal;

%Membuat Matrix Jacobi
J1=(x0-x)./(vp*sqrt((x0-x).^2+(y0-y).^2+(z0).^2));
J2=(y0-y)./(vp*sqrt((x0-x).^2+(y0-y).^2+(z0).^2));
J3=z0./(vp*sqrt((x0-x).^2+(y0-y).^2+(z0).^2));
J=[J1 J2 J3];

%Menghitung Error/Misfit
Erms=sqrt(sum((tcal-t).^2)./N);
dm=inv(J'*J)*J'*dt;
xp=[xp;x0];
yp=[yp;y0];
zp=[zp;z0];

%Model Perturbation
x0=x0+dm(1,1);
y0=y0+dm(2,1);
z0=z0+dm(3,1);

k=k+1;

if Erms <=0.001
break
end
end
initial_x=xp(1,1);
initial_y=yp(1,1);
initial_z=zp(1,1);
z=zeros(N,1);

pusat_x=x0;
pusat_y=y0;
pusat_z=z0;
figure(1)
plot3(xp,yp,zp,'--or','markerfacecolor','b')
text(initial_x,initial_y,initial_z,'Model Awal','verticalalignment','top')
hold on
plot3(x0,y0,z0,'kp','markerfacecolor','y','markersize',14)
text(pusat_x,pusat_y,pusat_z,'Pusat Gempa','verticalalignment','top')
hold on
plot3(x,y,z,'^b','markerfacecolor','b')
text(x,y,z,code,'verticalalignment','top')
grid on
set(gca,'ZDir','rev')
xlabel('Bujur');ylabel('Lintang');zlabel('Kedalaman');
fprintf(1,'Origin Time=%6.2f\n',t0_i)
fprintf(1,'Longitude=%6.2f\n',x0)
fprintf(1,'Latitude=%6.2f\n',y0)
fprintf(1,'Depth=%6.2f\n',z0)
fprintf(1,'Iterasi=%6d\n',k)


-> ganti text miring dengan data anda,,dari notepad dengan extension .dat