
clear
close all

% Program for the earthquake hypocentre location problem
% Problem taken from Stein + Wysession, Introduction to seismology, Chapter
% 7.2


nr = 10;     % number of receivers
nd = 3;      % number of space dimensions

% Initialize Station coordinates [km]
xr=[ 35 9 0; 
    -44 10 0;
    -11 -25 0;
     23 -39 0;
     42 -27 0;
    -12 50 0;
    -45 16 0;
      5 -19 0;
     -1 -11 0;
     20 11 0];
 
 

% Velocity model [km/s]
v=5;
 
% true model vector m (x, y, z, origin time)
m0=[0 0 10 0]';

% Calculate true data using hypo_forward.m
d0=hypo_forward(m0,xr,v)';

% End of initialization of true noise-free data

% Definite initial model for inverse problem
m=[3 4 20 2]';

maxni = 10;  % maximum number of iterations

% Iterative inversion

for ii=1:maxni,
    
    % calculate data difference for current model m
    dd = d0 - hypo_forward(m,xr,v)'; 
    
    % Calculate misfit energy
    err(ii) = (dd'*dd);
    disp(' -----------------------------------------------')
    disp(sprintf(' Misfit energy at iteration %i : %g ',ii,err(ii)));
    disp(' True solution vs. estimate : ')
    disp(sprintf(' x_true = %g km, x_est = %g km',m0(1),m(1))) 
    disp(sprintf(' y_true = %g km, y_est = %g km',m0(2),m(2))) 
    disp(sprintf(' z_true = %g km, z_est = %g km',m0(3),m(3))) 
    disp(sprintf(' t_true = %g  s, t_est = %g  s',m0(4),m(4))) 
    
    % initialize the system matrix G with the partial derivatives
    G = hypo_gradient(m,xr,v);
    
    % Solve the inverse problem for dm
    dm = inv(G'*G)*G'*dd;
    % update model
    m=m+dm;      

    
    
end    

plot(log10(err))
ylabel(' Log Misfit energy (s^2) ')
xlabel(' Iteration number ')
title(' Convergence behaviour - noise-free case ')