dataSeis=load('tunguST.txt');
[utm_x,utm_y,zoneUtm]=deg2utm(dataSeis(:,2),dataSeis(:,1));

A=[utm_x, utm_y];

writeIt(A,'utmCoorArray.txt');

refX=A(6,2); % reference point ARRAY x (station TBUG ) direction (using longitude)
refY=A(6,1); % reference point ARRAY y (station TBUG ) direction (using latitude)

% Calculate distances between the reference point (station TBUG) and the
% other stations.
distx=(refX-A(:,2))/1000;
disty=(refY-A(:,1))/1000;

% Create slowness matrix
Smin=-1/.5;
Smax=1/.5;
incr=.1;
Sx=Smin:incr:Smax;
Sy=Sx;

% Calculate power for each value of slowness (Sx,Sy)
power=zeros(length(Sx),length(Sy));
 % tmp_seistraceTLEF = seisV_TLEF-mean(seisV_TLEF);
 % tmp_seistraceTBAG = seisV_TBAG-mean(seisV_TBAG);
 % tmp_seistraceTGRS = -(seisV_TGRS-mean(seisV_TGRS));
 % tmp_seistraceTBUG = seisV_TBUG-mean(seisV_TBUG);
 % tmp_seistraceTREX = seisV_TREX-mean(seisV_TREX);
 

% 1Hz-5Hz for LP events 
 
sps = hdr_TLAV.sps; % get sample rate from header
npoles = 2; % specify order of filter
low = 1; % corner frequency at 1 Hz
high = 5; % corner frequency at 3 Hz
[B,A] = butter(npoles,[low high]/(sps/2));
 
 tmp_seistraceTLAV = filter(B,A,seisV_TLAV-mean(seisV_TLAV));
 tmp_seistraceTLEF = filter(B,A,seisV_TLEF-mean(seisV_TLEF));
 tmp_seistraceTMAR = filter(B,A,seisV_TMAR-mean(seisV_TMAR));
 tmp_seistraceTBAG = filter(B,A,seisV_TBAG-mean(seisV_TBAG));
 tmp_seistraceTGRS = filter(B,A,seisV_TGRS-mean(seisV_TGRS));
 tmp_seistraceTBUG = filter(B,A,seisV_TBUG-mean(seisV_TBUG));
 %tmp_seistraceTREX = filter(B,A,seisV_TREX-mean(seisV_TREX));

 
 %%
 
% for l=1:length(Sx)
%     for m=1:length(Sy)
%         shift(l,m)=tmp_seistraceTLAV-(Sx(l)*distx)-(Sy(m)*disty);
%     end
% end

window=10*60*100;
interval=length(tmp_seistraceTLAV)/window;

arrayTLAV=reshape(tmp_seistraceTLAV,interval,window);
arrayTLEF=reshape(tmp_seistraceTLEF,interval,window);
arrayTMAR=reshape(tmp_seistraceTMAR,interval,window);
arrayTBAG=reshape(tmp_seistraceTBAG,interval,window);
arrayTGRS=reshape(tmp_seistraceTGRS,interval,window);
arrayTBUG=reshape(tmp_seistraceTBUG,interval,window);

power=cell(interval,1);
%power(1:interval)={[]};
for i=1:1
    pow=zeros(length(Sx),length(Sy));
    for l=1:length(Sx)
        for m=1:length(Sy)
            %time_shift  = round(100*(-(Sx(l)*distx)-(Sy(m)*disty)));
            shiftTLAV=circshift(arrayTLAV(i,:), round(100*(-(Sx(l)*distx(1))-(Sy(m)*disty(1)))));
            shiftTLEF=circshift(arrayTLEF(i,:), round(100*(-(Sx(l)*distx(2))-(Sy(m)*disty(2)))));
            shiftTMAR=circshift(arrayTMAR(i,:), round(100*(-(Sx(l)*distx(3))-(Sy(m)*disty(3)))));
            shiftTBAG=circshift(arrayTBAG(i,:), round(100*(-(Sx(l)*distx(4))-(Sy(m)*disty(4)))));
            shiftTGRS=circshift(arrayTGRS(i,:), round(100*(-(Sx(l)*distx(5))-(Sy(m)*disty(5)))));
            %shiftTREX=circshift(tmp_seistraceTREX, round(100*(-(Sx(l)*distx(7))-(Sy(m)*disty(7)))));
            Tstack=shiftTLAV+shiftTLEF+shiftTMAR+shiftTBAG+shiftTGRS+ arrayTBUG(i,:); % lahaar event
            %Tstack=shiftTLAV+shiftTLEF+shiftTMAR+shiftTBAG+shiftTGRS+shiftTREX+tmp_seistraceTBUG;
            %Tstack = shiftTLEF + shiftTBAG + shiftTGRS + shiftTREX + tmp_seistraceTBUG;
            pow(m,l)=sum(Tstack.^2);
    
     %   maxp=[Sxpower
%        shiftTLEF
            %tmp_seistraceTLAV-(Sx(l)*distx)-(Sy(m)*disty);
        end
    end
    power{i}={pow};
    C=max(max(pow));
    %C=max(max(pow))
    [X,Y,C]=find(pow==C);
    Sxmax=Sx(X);
    Symax=Sy(Y);

    azimuth(i)=atand(Symax(i)/Sxmax(i));
    velocity(i)=1/(sqrt((Symax(i)^2)+(Sxmax(i))^2));

    subplot(i,1,i)
    norm_power = power/max(power(:));
    contour(Sx,Sy,norm_power,[.95:.003:.999])
    xlabel('S_x(s/km)');
    ylabel('S_y(s/km)');
    axis equal
    axis square
    xlim([-2,2])
    ylim([-2,2])

    xtick=-2:2;
    ytick=-2:2;
    set(gca,'XTick',xtick,'YTick',ytick)
    grid on
    gtext(['Velocity: ' num2str(velocity(i)) 'km/s'])
    gtext(['Azimuth: ' num2str(azimuth(i)) '^o'])



end
title(['2D slowness vector histogram - ' ' start time: ' num2str(start_time)])
C=max(max(power));
[X,Y,C]=find(power==C);
Sxmax=Sx(X);
Symax=Sy(Y);

azimuth=atand(Symax/Sxmax);
velocity=1/(sqrt((Symax^2)+(Sxmax)^2));



close all

%%

% for l=1:length(Sx)
%     for m=1:length(Sy)
%         shift(l,m)=tmp_seistraceTLAV-(Sx(l)*distx)-(Sy(m)*disty);
%     end
% end

for l=1:length(Sx)
    for m=1:length(Sy)
        %time_shift  = round(100*(-(Sx(l)*distx)-(Sy(m)*disty)));
        shiftTLAV=circshift(tmp_seistraceTLAV, round(100*(-(Sx(l)*distx(1))-(Sy(m)*disty(1)))));
        shiftTLEF=circshift(tmp_seistraceTLEF, round(100*(-(Sx(l)*distx(2))-(Sy(m)*disty(2)))));
        shiftTMAR=circshift(tmp_seistraceTMAR, round(100*(-(Sx(l)*distx(3))-(Sy(m)*disty(3)))));
        shiftTBAG=circshift(tmp_seistraceTBAG, round(100*(-(Sx(l)*distx(4))-(Sy(m)*disty(4)))));
        shiftTGRS=circshift(tmp_seistraceTGRS, round(100*(-(Sx(l)*distx(5))-(Sy(m)*disty(5)))));
        shiftTREX=circshift(tmp_seistraceTREX, round(100*(-(Sx(l)*distx(7))-(Sy(m)*disty(7)))));
        Tstack=shiftTLAV+shiftTLEF+shiftTMAR+shiftTBAG+shiftTGRS+shiftTREX+tmp_seistraceTBUG;
        %Tstack = shiftTLEF + shiftTBAG + shiftTGRS + shiftTREX + tmp_seistraceTBUG;
        power(m,l)=sum(Tstack.^2);
     %   maxp=[Sxpower
%        shiftTLEF
        %tmp_seistraceTLAV-(Sx(l)*distx)-(Sy(m)*disty);
    end
end

C=max(max(power));
[X,Y,C]=find(power==C);
Sxmax=Sx(X);
Symax=Sy(Y);

azimuth=atand(Symax/Sxmax);
velocity=1/(sqrt((Symax^2)+(Sxmax)^2));



close all





%%
figure
norm_power = power/max(power(:));
contour(Sx,Sy,norm_power,[.95:.003:.999])
xlabel('S_x(s/km)');
ylabel('S_y(s/km)');
title(['2D slowness vector histogram - ' ' start time: ' num2str(start_time)])
axis equal
axis square
xlim([-2,2])
ylim([-2,2])

xtick=-2:2;
ytick=-2:2;
set(gca,'XTick',xtick,'YTick',ytick)
grid on
gtext(['Velocity: ' num2str(velocity) 'km/s'])
gtext(['Azimuth: ' num2str(azimuth) '^o'])

%contour(Sx,Sy,norm_power)

%power=