%PROGAM PENENTUAN EPICENTER



clear 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');