Selasa, 03 Juli 2012

Matlab: Progam Penentuan Epicenter

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

0 komentar:

Posting Komentar

next page