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);
%  tmp_seistraceTLAV = seisV_TLAV-mean(seisV_TLAV);
%  tmp_seistraceTMAR = seisV_TMAR-mean(seisV_TMAR);
%  

% 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

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=shiftTLAV+shiftTLEF+shiftTMAR+shiftTBAG+shiftTGRS+tmp_seistraceTBUG;
        % Tstack=shiftTLAV+shiftTMAR+shiftTGRS+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
window = 3600*5/60/60; % 24 hours or 86400 s of data
dasTBUG = 'AA32';
%[seisV_TBUGlg,hdr_TBUGlg] = segycut(start_time,window,dasTBUG,'1');
%tslg = (1:length(seisV_TBUGlg))/hdr_TBUGlg.sps; % create time series vector in seconds
% tmp_seistraceTBUGlg = filter(B,A,seisV_TBUGlg-mean(seisV_TBUGlg));
%tmp_seistraceTBUGlg=seisV_TBUGlg-mean(seisV_TBUGlg);

subplot(2,1,1)
plot(ts,seisV_TBUG,'k')
ylabel('counts')
title('TBUG')
subplot(2,1,2)
norm_power = power/max(power(:));

%contour(Sx,Sy,norm_power,[.95:.003:.999])
contour(Sx,Sy,norm_power,25)
colorbar
%contour(Sx,Sy,norm_power)
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=