clear

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




fnamein = '2000.001.FNET.pick.dat';
fid = fopen(fnamein,'r');

YEAR = fscanf(fid,'%s',1);
jjday = fscanf(fid,'%d',1);
ARRAY = fscanf(fid,'%s\n',1);
mstla = fscanf(fid,'%f',1);
mstlo = fscanf(fid,'%f\n',1);
stime = fscanf(fid,'%f',1);
dtime = fscanf(fid,'%f\n',1);
sang = fscanf(fid,'%f',1);
dang = fscanf(fid,'%f\n',1);
dumm = fscanf(fid,'%d',3);
xT = fscanf(fid,'%f',dumm(3));
ifin = zeros(dumm(1),dumm(2),dumm(3));
afin = zeros(dumm(1),dumm(2),dumm(3));
for i = 1:dumm(1)
    for j = 1:dumm(2)
        for k = 1:dumm(3)
            ifin(i,j,k)=fscanf(fid,'%d',1);
            afin(i,j,k)=fscanf(fid,'%e',1);
        end
    end
end

angle = sang:dang:dang*(dumm(2)-1);
NT = length(xT);
dT = dtime;

disp('Reading... Done');

% Take reference frequency. (240 here because it is likely having a large 
% amplitude due to Airy phase

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

% Selecting possible picks for detecting Rayleigh Wave.

mm_fin = size(ifin);
tfin = ifin*dT-stime;

kk1 = 1;
for j = 1 : mm_fin(2) % loop for back azimuth.
%for j = 5 : 5
mm_ref_len = sum(sign(tfin(:,j,7))); % Number of picks for gV(7)

for i = 1 : mm_ref_len
%for i = 30 : 30
    if (tfin(i,j,NT) > 4000) & ( tfin(i,j,NT) < 96800) % ignore tapering.
       ref_time = tfin(i,j,NT);
           time_min = ref_time - 100;
           time_max = ref_time + 700;
           k1 = 0;
           for k = 1 : mm_fin(3)
               k2 = 0;
               for i1 = 1 : mm_fin(1)  
               if (tfin(i1,j,k) < time_max) & (tfin(i1,j,k) > time_min)
                  k1 = k1 + 1;
                  k2 = k2 + 1;
                  if (k2 >= 2)
                     k1 = k1 - 1;
                  else
                  ttime(k) = tfin(i1,j,k);                  
                  fin_array(kk1,k)=i1;
                  end
               end
               end
           end

           if (k1 == NT)
              fin_array(kk1,k+1)=j;
              [otime,dist,qual] = cal_GV_DIST(ttime,gV,3);
              dumm = sprintf('%5.3f %5.3f %5.3f',ref_time,dist,qual(1));
              disp(dumm);
              fin_array(kk1,k+2)=qual(1);
              kk1 = kk1 + 1;
           end
    end
end
end

%ttime = [51490,51490,51450,51430,51410,51380,51300];

mmf = size(fin_array);
err_cre = 300;
ii = 1;
for i = 1:mmf(1)-1
    if (fin_array(i,NT+2) <= err_cre)
        fin_array1(ii,:)=fin_array(i,:);   
        fintime(ii) = ifin(fin_array(i,3),fin_array(i,NT+1),3)*dT;
        finbaz(ii) = angle(fin_array(i,NT+1));
    for j = 1:NT
        tdumm(j)=afin(fin_array(i,j),fin_array(i,NT+1),j);
    end
        finamp(ii) = mean(tdumm);
        ii = ii + 1;     
    end
end

fintime = fintime/(24*3600)+jday1;

% READING CMT catalog

[eventid,jday,moment,hour,tmin,sec,lat,lon]=readCMT('TEMPCMT.dat');
tmw=log10(moment)/1.5-10.73;

mm1=size(tmw);
j1=1;
for j=1:mm1(1)
    if (tmw(j) >= cre_mw1) & (tmw(j) <= cre_mw2)
    if (jday(j) >= jday1) & (jday(j) <= jday2)
       ttlat = -lat(j);
       ttlon = 180+lon(j);  
       event_cmt(j1,:)=eventid(j,:);
       time_cmt(j1)=((jday(j)).*24+hour(j)+tmin(j)./60+sec(j)./3600);
       xepi = cal_epi(mstla,mstlo,lat(j),lon(j));
       xbaz = cal_baz(mstla,mstlo,lat(j),lon(j));
       addtime = xepi*6371*pi/180/gV(3);
       time_cmt1(j1) = time_cmt(j1)+addtime;
       lat_cmt(j1)=lat(j);
       lon_cmt(j1)=lon(j);
       baz_cmt(j1)=xbaz;
       mw_cmt(j1)=log10(moment(j))/1.5-10.73;
       j1=j1+1;
    end
    end
end
