

clear
close all

% Inversion for Earth model from reflection travel times

nr = 2;     % number of receivers


% Initialize Station coordinates [km]
xr=[ 10 20 ];
 
% true model vector m (x, y, z, origin time)
m0=[5 10]';

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

% End of initialization of true noise-free data

% Definite initial model for inverse problem
m=[4 8]';

maxni = 10;  % maximum number of iterations

% Iterative inversion

for ii=1:maxni,
    
    % calculate data difference for current model m
    dd = d0 - refl_forward(m,xr)'; 
    

    % initialize the system matrix G with the partial derivatives
    G = refl_gradient(m,xr);
    
    % Solve the inverse problem for dm
    dm = inv(G'*G)*G'*dd;
 %   dm = inv(G)*dd;
    % update model
    m=m+dm;      

    % 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(' v_true = %g km, v_est = %g km',m0(1),m(1))) 
    disp(sprintf(' h_true = %g km, h_est = %g km',m0(2),m(2))) 
    
    
end    

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