clear

cre_mw1 = 5.8;
cre_mw2 = 8.0;
jday1 = 1;
jday2 = 1;




fnamein = '2000.001.FNET.pick.dat';
[YEAR,jjday,ARRAY,mst,time,angle,xT,ifin,afin] = readpick(fnamein);
mstla = mst(1);
mstlo = mst(2);

NT = length(xT);
dT = time(2)-time(1);

[dumm1,dumm2,gV]=cal_CQU_228(1./xT);

% Selecting possible picks for detecting Rayleigh Wave.

mm_fin = size(ifin);
tfin = ifin*dT+time(1);
iref = 3;  % reference frequency. In this case, 240 sec

idet = 1;

for j = 1 : mm_fin(2) % loop for backazimuth
%for j = 36 : 36

imax_ref = sum(sign(ifin(:,j,iref)));
for i1 = 1 : imax_ref % loop for reference frequency.
%for i1 = 27 : 27 % loop for reference frequency.

i1_tref = ifin(i1,j,iref);

if (tfin(i1,j,iref) > 0) & (tfin(i1,j,iref) < 24*3600)
ixT=zeros(mm_fin(3),1);
for k = 1 : mm_fin(3) % loop for xT
itempxT = 1;

for i = 1 : mm_fin(1) % loop for time
    itemp = ifin(i,j,k);
    if (itemp > i1_tref-100) & (itemp < i1_tref+100) & (itemp ~= 0)
       ixT(k)=ixT(k)+1;
       itTEMP(k,itempxT) = tfin(i,j,k);
       atTEMP(k,itempxT) = afin(i,j,k);
       itempxT = itempxT+1;
    end

end % End loop for time
end % End loop for xT

if (sum(sign(ixT)) == NT)
itTEMP(iref,:) = 0;
atTEMP(iref,:) = 0;
itTEMP(iref,1) = tfin(i1,j,iref);
atTEMP(iref,1) = afin(i1,j,iref);
[dumm,iatmax]=max(atTEMP');
if (size(atTEMP,2) == 1)
iatmax = zeros(1,NT);
iatmax = iatmax + 1;
dumm = atTEMP';
end


for i = 1:NT
ttime(i) = itTEMP(i,iatmax(i));
end
[otime,dist,qual]=cal_GV_DIST1(ttime,gV,iref);
fintime(idet)=otime;
findist(idet)=dist;
finbaz(idet)=angle(j);
finerror(idet,:)=qual;
fin_amp(idet,:)=dumm;
fin_time(idet,:)=ttime;
idet=idet+1;
itTEMP = 0;
atTEMP = 0;
end

end % End if
end % End loop for reference freq.
end % End loop for back azimuth.


idetfin = 1;
err1 = 600;
for ii = 1 : idet-1
temp_err1 = finerror(ii,2:NT+1);
[dumm,isort1]=sort(abs(temp_err1));
temp_err2 = sqrt(sum(abs(dumm(2:NT-1)).^2)/(NT-1));
temp_err3 = sqrt(sum(abs(dumm(2:NT-2)).^2)/(NT-2));
if (fin_amp(ii,NT) >= 5e-10)
if (finerror(ii,1) <= err1)
fin_baz1(idetfin)=finbaz(ii);
fin_error1(idetfin,:)=finerror(ii,:);
fin_amp1(idetfin,:)=fin_amp(ii,:);
fin_time1(idetfin,:)=fin_time(ii,:);

idetfin = idetfin + 1;
end
end


end
