% Source Wavelet Inversion
%
% Written by C. Boonyasiriwat, UTAM, 2007

clear;clc;

% Initialize problem parameters
nt = 1500;
nz = 301;
dx = 5;
dt = 0.001;
t = (0:nt-1)'*dt;
z = (0:nz-1)'*dx;

% Setup velocity profile
c = 2000*ones(nz,1);
c(51:80) = 1600;
c(121:150) = 2500;
c0 = c;

p0 = zeros(nz,1);
p1 = p0;
p = p0;

% True source
src_true = ricker(nt, 5, dt);
src_true = [src_true'; zeros(nt-length(src_true),1)];

% Generate data
zs = 5;
p0(zs) = src_true(1);
p1(zs) = src_true(2);
[seis_true,p_all] = modeling(nz,nt,src_true,c,zs,dx,dt,p0,p1);

% Trial source
src_trial = ricker(nt, 10, dt);
src_trial = [src_trial'; zeros(nt-length(src_trial),1)];

n_iter = 10;
residual = zeros(n_iter);
set(gcf,'DoubleBuffer','on');
for i=1:n_iter
    
    % Compute synthetic data using trial source
    [seis_trial,p_all] = modeling(nz,nt,src_trial,c,zs,dx,dt,p0,p1);
    res = norm(seis_true-seis_trial);
    residual(i) = res;
    subplot(312);plot(t,seis_true,'b',t,seis_trial,'r');
    legend('Observed Data','Synthetic Data');
    xlabel('Time (s)');ylabel('Amplitude');
    title('Observed and Synthetic Data');

    % Back propagate the residual
    dp = seis_true - seis_trial;
    dp = flipud(dp)/max(abs(dp));
    p0(zs) = dp(1);
    p1(zs) = dp(2);
    [dseis,dp_all] = modeling(nz,nt,dp,c,zs,dx,dt,p0,p1);
    
    % Compute velocity update
    ds = flipud(dseis);
    
    % Determine a step length that reduces the residual.
    alpha = 1e-2;
    for j=1:20
        src = src_trial + alpha*ds;
        p0(zs) = src(1);
        p1(zs) = src(2);
        [seis_trial,p_all] = modeling(nz,nt,src,c,zs,dx,dt,p0,p1);
        res_new = norm(seis_true-seis_trial);
        if res_new < res
            break;
        else
            alpha = alpha/2;
        end;
    end;
    
    % Update the source
    src_trial = src;

    subplot(311);plot(t,src_true,'b',t,src_trial,'r');
    legend('True Source','Inverted Source');
    xlabel('Depth (m)');ylabel('Velocity (m/s)');
    title('True and Inverted Sources');
    
    subplot(313);plot(residual);title('Residual over Iteration');
    xlabel('Iteration');ylabel('Residual');
    drawnow;
end;
