clear;clc;

f1 = 5;
f2 = 10;
nt = 1500;
nz = 301;
c = 2000*ones(nz,1);
c(51:80) = 1600;
c(121:150) = 2500;
c0 = c;
p0 = zeros(nz,1);
p1 = p0;
p = p0;
dx = 5;
dt = 0.001;
rick_f5 = ricker(nt, f1, dt);
rick_f5 = [rick_f5'; zeros(nt-length(rick_f5),1)];
t = (0:nt-1)'*dt;
z = (0:nz-1)'*dx;
set(gcf,'DoubleBuffer','on');

zs = 5;
p0(zs) = rick_f5(1);
p1(zs) = rick_f5(2);
[seis_f5,p_all_f5] = modeling(nz,nt,rick_f5,c,zs,dx,dt,p0,p1);

rick_f10 = ricker(nt, f2, dt);
rick_f10 = [rick_f10'; zeros(nt-length(rick_f10),1)];
p0(zs) = rick_f10(1);
p1(zs) = rick_f10(2);
[seis_f10,p_all_f10] = modeling(nz,nt,rick_f10,c,zs,dx,dt,p0,p1);

c = 2000*ones(nz,1);
n_iter = 50;
residual = zeros(n_iter);

% Low frequency inversion
for i=1:n_iter
    
    % Compute synthetic data using trial velocity
    [seis0,p0_all] = modeling(nz,nt,rick_f5,c,zs,dx,dt,p0,p1);
    res = norm(seis_f5-seis0);
    residual(i) = res;
    subplot(323);plot(t,seis_f5,'b',t,seis0,'r');
    legend('Input Trace','Synthetic Trace');
    xlabel('Time (s)');ylabel('Amplitude');

    % Back propagate the residual
    dp = seis_f5 - seis0;
    dp = flipud(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
    dc = sum(p_all_f5.*fliplr(dp_all),2);
    
    % Determine a step length that reduces the residual.
    alpha = 1;
    for j=1:5
        c1 = c - alpha*dc/100;
        p0(zs) = rick_f5(1);
        p1(zs) = rick_f5(2);
        [seis0,p0_all] = modeling(nz,nt,rick_f5,c1,zs,dx,dt,p0,p1);
        res_new = norm(seis_f5-seis0);
        if res_new < res
            break;
        else
            alpha = alpha/2;
        end;
    end;
    
    % Update the velocity
    c = c - alpha*dc/100;

    subplot(321);plot(z,c0,'b',z,c,'r');legend('True Velocity','Inverted Velocity');
    xlabel('Depth (m)');ylabel('Velocity (m/s)');
    title('Low frequency inversion');
    axis([0 1500 1000 3000]);
    drawnow;
end;
subplot(325);plot(residual);title('Residual over iteration');
xlabel('Iteration');ylabel('Residual');

% High frequency inversion
residual = zeros(n_iter);
for i=1:n_iter
    
    % Compute synthetic data using trial velocity
    [seis0,p0_all] = modeling(nz,nt,rick_f10,c,zs,dx,dt,p0,p1);
    res = norm(seis_f10-seis0);
    residual(i) = res;
    subplot(324);plot(t,seis_f10,'b',t,seis0,'r');
    legend('Input Trace','Synthetic Trace');
    xlabel('Time (s)');ylabel('Amplitude');

    % Back propagate the residual
    dp = seis_f10 - seis0;
    dp = flipud(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
    dc = sum(p_all_f10.*fliplr(dp_all),2);
    
    % Determine a step length that reduces the residual.
    alpha = 1;
    for j=1:5
        c1 = c - alpha*dc/100;
        p0(zs) = rick_f10(1);
        p1(zs) = rick_f10(2);
        [seis0,p0_all] = modeling(nz,nt,rick_f10,c1,zs,dx,dt,p0,p1);
        res_new = norm(seis_f10-seis0);
        if res_new < res
            break;
        else
            alpha = alpha/2;
        end;
    end;
    
    % Update the velocity
    c = c - alpha*dc/100;

    subplot(322);plot(z,c0,'b',z,c,'r');legend('True Velocity','Inverted Velocity');
    axis([0 1500 1000 3000]);
    xlabel('Depth (m)');ylabel('Velocity (m/s)');
    title('High frequency inversion');
    drawnow;
end;
subplot(326);plot(residual);title('Residual over iteration');
xlabel('Iteration');ylabel('Residual');
