clear

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

fname_cmt_out = 'ref_cmt.dat';
fid1 = fopen(fname_cmt_out,'w');
fname_cmt_out1 = 'ref_cmt1.dat';
fid3 = fopen(fname_cmt_out1,'w');

cre_mw1 = 5.0;
cre_mw2 = 8.0;
cre_jday1 = 1;
cre_jday2 = 31;
err_cre = 8e+3; % Error is tied to location error

mm1=size(tmw);
j1=1;
for j=1:mm1(1)
    if (tmw(j) >= cre_mw1) & (tmw(j) <= cre_mw2)
    if (jday(j) >= cre_jday1) & (jday(j) <= cre_jday2)
       ttlat = -lat(j);
       ttlon = 180+lon(j);  
       fprintf(fid1,'%s %5.3f %5.3f\n',char(eventid(j,:)),lat(j),lon(j));
       fprintf(fid3,'%s %5.3f %5.3f\n',char(eventid(j,:)),ttlat,ttlon);
       event_cmt(j1,:)=eventid(j,:);
       time_cmt(j1)=((jday(j)).*24+hour(j)+tmin(j)./60+sec(j)./3600)/24;
       lat_cmt(j1)=lat(j);
       lon_cmt(j1)=lon(j);
       mw_cmt(j1)=log10(moment(j))/1.5-10.73;
       j1=j1+1;
    end
    end
end

fclose(fid1);
fclose(fid3);

% Reading picking file.

fnamein = 'BAND2_LOC.OUT';
fid = fopen(fnamein,'r');

npeve = fscanf(fid,'%i',1);

for i = 1 : npeve
    time(i) = fscanf(fid,'%f',1);
    elat(i) = fscanf(fid,'%f',1);;
    elon(i) = fscanf(fid,'%f',1);
    error(i) = fscanf(fid,'%e',1);

%    index(i) = fscanf(fid,'%i',1);
    time01(i) = fscanf(fid,'%f',1);
    time02(i) = fscanf(fid,'%f',1);
    time03(i) = fscanf(fid,'%f',1);
    baz01(i) = fscanf(fid,'%f',1);
    baz02(i) = fscanf(fid,'%f',1);
    baz03(i) = fscanf(fid,'%f',1);
    amp01(i) = fscanf(fid,'%e',1);
    amp02(i) = fscanf(fid,'%e',1);
    amp03(i) = fscanf(fid,'%e',1);
end
fclose(fid);


amp = (amp01+amp02+amp03)/3;
%amp = amp01;
fname_pick_out = 'pick.dat';
fid2 = fopen(fname_pick_out,'w');


% Remove detections with identical points.(??)

ttime0 = 0;
jj = 1;
for j = 1 : npeve-1
    ttime1 = time01(j);
    ttime2 = time01(j+1);
    if (ttime0 ~= ttime1)
        j1 = j;
        ttime0 = ttime1;
    end
    if (ttime1 ~= ttime2)
       j2 = j;
       [dumm,imin]=min(error(j1:j2));
       iimin = j1+imin-1;
       fin_index(jj)=iimin;
       fin_error(jj)=dumm;
       fin_lat(jj) = elat(iimin);
       fin_lon(jj) = elon(iimin);
       fin_time(jj) = time(iimin);
       fin_amp(jj) = amp(iimin);
       jj =jj + 1;
    end
end

fin_len = length(fin_error);

% We need option which can remove R2 and other multi-circling phases.

Erad = 6371;
refgvel = 3.8;
xdtime = 180*Erad*pi/180/refgvel/(24*3600);
%cre_r2;

% Index 0: Event 1: Multiple circling.

index = zeros(fin_len,1);

for i = 1 : fin_len-1
if (index(i) == 0)
time0 = fin_time(i);
amp0 = fin_amp(i);
lat0 = fin_lat(i);
lon0 = fin_lon(i);
lat1 = -lat0;
lon1 = sign(lon0)*(abs(lon0)-180);

k = i + 1;
    for j = k : fin_len
       time1 = fin_time(j);
       amp1 = fin_amp(j);
       if (time1 <= time0+xdtime*(1+1/5)) & (time1 >= time0+xdtime*(1-1/5)) 
       lat2 = fin_lat(j);
       lon2 = fin_lon(j);
       epi1 = cal_epi(lat0,lon0,lat2,lon2);
       epi2 = cal_epi(lat1,lon1,lat2,lon2);
%       [time0,time1,min(epi1,epi2),lat0,lon0,lat1,lon1,lat2,lon2]
%       pause
       if (min(epi1,epi2) <= 30) & (amp1 < amp0)
       time0 = time1;
       amp0 = amp1;
       index(j) = 1;
       end
       end
       end
end
end

k1 = 1;
for k = 1 : fin_len
    if (index(k) == 0) & (fin_error(k) <= err_cre)
        fin_lat1(k1)= fin_lat(k);
        fin_time1(k1) = fin_time(k);
        fin_lon1(k1) = fin_lon(k);
        fin_error1(k1) = fin_error(k);
        fin_index1(k1) = fin_index(k);
        fin_amp0(k1) = fin_amp(k);
        k1 = k1 + 1;
    end
end

for k = 1 : fin_len
    if (fin_error(k) <= err_cre) & (index(k) == 0)      
       fprintf(fid2,'%5.3f %5.3f %5.3f\n',fin_lat(k),fin_lon(k),fin_error(k));
    end
end

fclose(fid2);

% Final check

mm_fin = size(fin_lat1);
mm_cmt = size(lat_cmt);

ttempout = 'loc_eve.dat';
fid5=fopen(ttempout,'w');
ttempout1 = 'loc_non_eve.dat';
fid6 = fopen(ttempout1,'w');

k = 1;
kk = 1;
for i = 1 : mm_fin(2)
ii = 0;
    for j = 1 : mm_cmt(2)
        timediff = abs(fin_time1(i)-time_cmt(j));
        difepi = cal_epi(fin_lat1(i),fin_lon1(i),lat_cmt(j),lon_cmt(j));
        if (timediff <= 0.03) & (difepi <= 40)
        lat00(k) = fin_lat1(i);
        lon00(k) = fin_lon1(i);
        lat01(k) = lat_cmt(j);
        lon01(k) = lon_cmt(j);
        mw_cmt3(k) = mw_cmt(j);
        iik = fin_index1(i);
        fin_index2(k) = iik;
        fin_amp3(k) = fin_amp0(i); % Detected amplitude
        fin_error2(k) = fin_error1(i);
        fin_loc(k) = cal_epi(lat00(k),lon00(k),lat01(k),lon01(k));

fprintf(fid5,'%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f %8.3f %8.3f %8.3f\n', ...
time01(iik),time02(iik),time03(iik),baz01(iik),baz02(iik),baz03(iik), ...
time_cmt(j),lat_cmt(j),lon_cmt(j));
        k = k + 1;
        ii = 1;
        break;
        end
    end
if (ii == 0) % no related events
   iik1 = fin_index1(i);
   ttime02(kk) = fin_time1(i);
   tamp02(kk) = fin_amp0(i);
   lat02(kk) = fin_lat1(i);
   lon02(kk) = fin_lon1(i);
   error02(kk) = fin_error1(i);

fprintf(fid6,'%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f %8.3f %8.3f %8.3f\n', ...
time01(iik1),time02(iik1),time03(iik1),baz01(iik1),baz02(iik1),baz03(iik1), ...
ttime02(kk),lat02(kk),lon02(kk));

kk = kk + 1;
end
end

fclose(fid5);
fclose(fid6);



fin_amp1 = log(fin_amp0)/1.5+15.5;

clf
subplot(6,1,1),plot(fin_time1,fin_amp1,'r.')
hold on
plot(time_cmt,mw_cmt,'bo')
plot(ttime02,log(tamp02)/1.5+15.5,'k*')
axis([1 6 3 9])
subplot(6,1,2),plot(fin_time1,fin_amp1,'r.')
hold on
plot(time_cmt,mw_cmt,'bo')
plot(ttime02,log(tamp02)/1.5+15.5,'k*')
axis([6 11 3 9])
subplot(6,1,3),plot(fin_time1,fin_amp1,'r.')
hold on
plot(time_cmt,mw_cmt,'bo')
plot(ttime02,log(tamp02)/1.5+15.5,'k*')
axis([11 16 3 9])
subplot(6,1,4),plot(fin_time1,fin_amp1,'r.')
hold on
plot(time_cmt,mw_cmt,'bo')
plot(ttime02,log(tamp02)/1.5+15.5,'k*')
axis([16 21 3 9])
subplot(6,1,5),plot(fin_time1,fin_amp1,'r.')
hold on
plot(time_cmt,mw_cmt,'bo')
plot(ttime02,log(tamp02)/1.5+15.5,'k*')
axis([21 26 3 9]);
subplot(6,1,6),plot(fin_time1,fin_amp1,'r.')
hold on
plot(time_cmt,mw_cmt,'bo')
plot(ttime02,log(tamp02)/1.5+15.5,'k*')
axis([26 32 3 9])
