%PROGAM PENENTUAN EPICENTERclear all; clc;
%Membuka file data
fid=fopen('stasiun_data.dat');
data=textscan(fid, '%s %f %f %f');
fclose(fid);
N=length(data);
code = data{1};
x = data{2};
y = data{3};
t = data{4};
xp=[]; yp=[];
%Model awal
x0 = input('Masukkan Longitude Awal : ');
y0 = input('Masukkan Latitude Awal : ');
vp=4; Erms=1;
%Pemodelan Inversi Menggunakan Iterasi
k=0;
while Erms > 0.00001
tp = sqrt((x0-x).^2+(y0-y).^2)/vp;
t0=abs(t-tp); t0=mean(t0);
tcal=t0+tp;
dt=t-tcal;
%Membuat matrix Jacobi
J1=(x0-x)./(vp*sqrt((x0-x).^2+(y0-y).^2));
J2=(y0-y)./(vp*sqrt((x0-x).^2+(y0-y).^2));
J=[J1 J2];
%Menghitung error/misfit
Erms = sqrt(sum((tcal-t).^2)./N);
%Permodelan Inversi
dm=inv(J'*J)*J'*dt;
xp=[xp;x0];
yp=[yp;y0];
%Plotting hasil
figure (1)
plot(x0,y0,'or','markerfacecolor','r')
hold on
plot(xp,yp,'--b')
hold on
plot(x,y,'^b','markerfacecolor','b')
text(x,y,code,'verticalalignment','top')
daspect([1 1 1])%memberikan gambaran bujur sangkar dengan ratio yang sama
axis([0 70 0 70])
%Model Perturbation
x0=x0+dm(1,1);
y0=y0+dm(2,1);
k=k+1;
if Erms <= 0.1
break
end
end
fprintf(1,'Epicenter => ') %1 maksudnya memunculkan text di cmd window
fprintf(1,'Lon ;%6.2f ',x0)
fprintf(1,'Lat ;%6.2f\n ',y0)
fprintf(1,'Jumlah Iterasi ;%2.0f\n ',k)
init_x = xp(1,1); % init maksudnya inisial (model awal)
init_y = yp(1,1);
figure (1)
text(init_x,init_y,'Model Awal','verticalalignment','top');
text(x0,y0,'Epicenter','verticalalignment','top');
%Membuka file data
fid=fopen('stasiun_data.dat');
data=textscan(fid, '%s %f %f %f');
fclose(fid);
N=length(data);
code = data{1};
x = data{2};
y = data{3};
t = data{4};
xp=[]; yp=[];
%Model awal
x0 = input('Masukkan Longitude Awal : ');
y0 = input('Masukkan Latitude Awal : ');
vp=4; Erms=1;
%Pemodelan Inversi Menggunakan Iterasi
k=0;
while Erms > 0.00001
tp = sqrt((x0-x).^2+(y0-y).^2)/vp;
t0=abs(t-tp); t0=mean(t0);
tcal=t0+tp;
dt=t-tcal;
%Membuat matrix Jacobi
J1=(x0-x)./(vp*sqrt((x0-x).^2+(y0-y).^2));
J2=(y0-y)./(vp*sqrt((x0-x).^2+(y0-y).^2));
J=[J1 J2];
%Menghitung error/misfit
Erms = sqrt(sum((tcal-t).^2)./N);
%Permodelan Inversi
dm=inv(J'*J)*J'*dt;
xp=[xp;x0];
yp=[yp;y0];
%Plotting hasil
figure (1)
plot(x0,y0,'or','markerfacecolor','r')
hold on
plot(xp,yp,'--b')
hold on
plot(x,y,'^b','markerfacecolor','b')
text(x,y,code,'verticalalignment','top')
daspect([1 1 1])%memberikan gambaran bujur sangkar dengan ratio yang sama
axis([0 70 0 70])
%Model Perturbation
x0=x0+dm(1,1);
y0=y0+dm(2,1);
k=k+1;
if Erms <= 0.1
break
end
end
fprintf(1,'Epicenter => ') %1 maksudnya memunculkan text di cmd window
fprintf(1,'Lon ;%6.2f ',x0)
fprintf(1,'Lat ;%6.2f\n ',y0)
fprintf(1,'Jumlah Iterasi ;%2.0f\n ',k)
init_x = xp(1,1); % init maksudnya inisial (model awal)
init_y = yp(1,1);
figure (1)
text(init_x,init_y,'Model Awal','verticalalignment','top');
text(x0,y0,'Epicenter','verticalalignment','top');
0 komentar:
Posting Komentar