load 'Wadi_Qudaid_1.mat';  % stored in csg
% wq=csg;

dim = size(csg);
disp(dim);
nt  = dim(1);
ns  = dim(2);
ng  = dim(3);
disp (nt);
% wq = reshape(csg, nt, ng, ns);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% show original data set
wq = zeros([nt  ng  ns]);

data = zeros([nt ng]);
trace = zeros([nt 1]);
data(:, :) = csg(:, :, 70);

% figure;
% imagesc( data, [-15 15]);
% colormap( gray );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Extracting First arrivals
figure;
farr = readbin('processed_pick.bin',ns,ng);
imagesc(farr);
xlabel('Receivers');
ylabel('Sources');
set(gca, 'xtick',[0 20 40 60 80 100 117]);
set(gca, 'xticklabel',[0 20 40 60 80 100 117]);
set(gca, 'ytick',[0 20 40 60 80 100 117]);
set(gca, 'yticklabel',[0 20 40 60 80 100 117]);
title('Picked first arrival matrix');

dwin = zeros([ns 2*ng]);
for i = 1:ns
    srcarr = floor(farr(i,i)/0.001)-5;
    for j = 1:ng                 
        t1  = floor(farr(i,j)/0.001)+1;        
        t2  = nt;
        t11 = t1 - srcarr;
        if( t11 < 1 )
            t11 = 1;
        end
        t22 = t11 + (t2-t1);        
        wq(t11:t22, j, i) = csg(t1:t2, j, i);
        
%         dwin(j+(i-1)*ns, 1) = i;
%         dwin(j+(i-1)*ns, 2) = j;
        dwin(i, 2*(j-1)+1) = 1 + (t11-1)*8;
        dwin(i, 2*(j-1)+2) = 1 + (t22-1)*8;        
    end
    disp('i ='); disp(i);
    disp(farr(i,i));
end
% writebin('wadi_data_window.bin',dwin);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% one CSG 
figure;
for it = 1:nt
    for ig = 1:ng
        data( it, ig ) = wq(it, ig, 70);
    end
end

% subplot( 311  );
imagesc( data );
colormap( gray);
title('csg 70');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% OBS windowed early arrival dataset
ext_trace = zeros([nt*8-7 1]);
ext_data = zeros([nt*8-7 ng]);

for is = 30:30 %1:ns
%     path = strcat('./OBS/wadi',num2str(is),'.bin');
    for ig=1:ng
        for it=1:nt
            trace(it) = wq(it,ig,is);        % for each trace
        end
        trace = trace/max(abs(trace)); % trace by trace normalization
        for it = 1:nt
            ext_trace(1+(it-1)*8) = trace(it);
            if(it<nt)
                for k = 1:1:7
                    lam = 0.125*k;
                    ext_trace(1+(it-1)*8+k) = (1-lam)*trace(it)+(lam)*trace(it+1);
                end
            end
        end
        ext_data(:,ig) = ext_trace(:);
    end
%     writebin( path, ext_data(1:2000, :) );
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%show one CSG
figure;
mycsg = ext_data(1:2000, :);
for ig = 1:ng
    mycsg(:,ig) = mycsg(:,ig)/max(abs(mycsg(:,ig)));
end
imagesc(mycsg);
colormap(gray);
hold on;

for ig = 1:ng
    x = ig;
    y = dwin(30, 2*(ig-1)+1);
    plot(x,y,'+r');
    hold on;
end
xlabel( 'Distance X(m)' );
ylabel( 'Time(s)' );
set( gca, 'xtick',[0.0 20 50 80 110] );
set( gca, 'xticklabel',[0.0 40 100 160 220] );
set( gca, 'ytick',[0 200 800 1400 2000] );
set( gca, 'yticklabel', [0 0.025 0.1 0.175 0.25] );
title( 'Common shot gather #30 with picked first arrivals' );
writebin( 'example_data.bin', mycsg );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Spectrum of a trace
L = 2000;
fs = 8000;
nfft = 3*L;
% yy = zeros([nt 1]);
yy = ext_data(1:2000, 30);
disp(size(yy));
fyy = fft(yy,nfft);
disp('fft size');
disp(size(fyy));
w = (0:length(fyy)-1)*fs/length(fyy);

figure;
plot(w,abs(fyy));
% plot(abs(fftshift(fftshift(fyy))));
axis([0 500 0 max(abs(fyy))]);
grid on;
xlabel('Hz');
ylabel('Amplitude');
title('Data spectrum (is = 30, ig = 30)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%