% Test m-code for looking at characteristics of events.

clear

YEAR = '2000';
ARRAY = 'GRSN';
Erad = 6371;
jday1 = 1;
jday2 = 31;

cre_mw1 = 5.5;
cre_mw2 = 8.0;

idet=1;
for jjday = jday1:jday2

if (jjday < 10)
JDAY = sprintf('00%i',jjday);
elseif (jjday < 100)
JDAY = sprintf('0%i',jjday);
else
JDAY = jjday;
end
disp(JDAY);
fnamein = sprintf('%s.%s.%s.pick1.dat',YEAR,JDAY,ARRAY);
fid = fopen(fnamein,'r');

YEAR = fscanf(fid,'%s',1);
jtemp = 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;
tfin = ifin*dT+stime;

[a,b,gV]=cal_CQU_228(1./xT);

[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) >= jjday) & (jday(j) < jjday+1)
       otime = hour(j)*3600+tmin(j)*60+sec(j);
       xepi = cal_epi(mstla,mstlo,lat(j),lon(j));
       xbaz = cal_baz(mstla,mstlo,lat(j),lon(j));
       ibaz = round(xbaz/dang)+1;
       addtime = xepi*Erad*pi/180./gV;
       fin_time_cmt(idet,:)=otime+addtime';
       fin_eveid(idet,:)=eventid(j,:);
       fin_mw(idet) = tmw(j);
       for jj = 1 : NT
           ipick = sum(sign(ifin(:,ibaz,jj)));           
           for k = 1:ipick-1
               tmp1 = tfin(k,ibaz,jj);
               tmp2 = tfin(k+1,ibaz,jj);
               tmp0 = fin_time_cmt(idet,jj);
               if (tmp1 < tmp0) & (tmp2 >= tmp0)
                  tmp = [tmp0-tmp1,tmp0-tmp2];
                  [dumm,imin]=min(abs(tmp));
                  if (imin ==1)
                  fin_time(idet,jj)=tmp1;
                  fin_amp(idet,jj)=afin(k,ibaz,jj);                  
                  else
                  fin_time(idet,jj)=tmp2;
                  fin_amp(idet,jj)=afin(k+1,ibaz,jj);   
                  fin_jday(idet) = jday(j);               
                  end
               
               end
           end
       end
          
    idet = idet + 1;
    end
    end
end




end % loop for jjday


idet1 = idet - 1;


for i = 1 : idet1
    temp = fin_time_cmt(i,:)-fin_time(i,:);
    temp2 = mean(abs(temp.^2));
    fin_error(i) = sqrt(temp2);
    [otime,dist,qual]=cal_GV_DIST1(fin_time(i,:),gV,3);
    fin_error1(i) = qual(1);
end
