% m-code for picking potential signals.
% Different picking algorithm will be tested.


clear
YEAR = '2000';
jday1 = 25;
jday2 = 25;
ARRAY1 = 'BDSN';
BAND = 'BAND1';
EXT = 'AMPEQ';
Erad = 6371;
refgvel = 3.8;
%amp_cre = 1.3e-8;
amp_fac = 2;

% Reading amp data.


jjday = jday1:jday2;
mmjday = length(jjday);

for i = 1:mmjday
jday = jjday(i);

if (jday < 10)
julday = sprintf('00%i',jday);
elseif (jday < 100)
julday = sprintf('0%i',jday);
else
julday = sprintf('%i',jday);
end

% Picking for first ARRAY
tempname01 = sprintf('%s.%s.%s.%s.%s.amp',YEAR,julday,ARRAY1,BAND,EXT);
[dumm,dumm1,tmstla,tmstlo,stime,sdtime,sang,sdang,tempout01] = ...
readamp(tempname01);

mm01 = size(tempout01);
time01 = 0 : sdtime : sdtime*(mm01(1)-1);
time = (time01 + stime)/(3600*24) + jday;

baz01 = 0 : sdang : sdang*(mm01(2)-1);
baz01 = baz01 + sang;
[tmaxamp01,imaxamp01]=max(tempout01');
maxamp = tmaxamp01';
maxbaz = baz01(imaxamp01);
amp_fac1 = 1.5;
amp_fac2 = 2;
idur4lta = 30;


 
mm0 = size(maxamp);
fac1 = 2;

ii=1;
for i =1 : mm0(1)-2;
    amp1 = maxamp(i);
    amp2 = maxamp(i+1);
    amp3 = maxamp(i+2);
    amp12 = amp1 - amp2;
    amp23 = amp2 - amp3;
    amp4 = amp12 * amp23;
    if (amp4 < 0) 
       finamp(ii)=maxamp(i+1);
       fintime(ii)=time(i+1);
       finbaz(ii)=maxbaz(i+1);
       iidetect(ii)=i+1;
       ii = ii + 1;
    end  
end


mm = size(finamp);

isign = 0;
ii = 1;
tmax0=max(finamp(1),finamp(2));
tmin0=min(finamp(1),finamp(2));
tmax3=tmax0-0.1*(tmax0-tmin0);
tmin3=tmin0+0.1*(tmax0-tmin0);

isign = 0;
i=1;
while (i < mm(2)-1)

% initialization.

amp1 = finamp(i);
amp2 = finamp(i+1);
amp12 = amp1 - amp2;


tmax1 = max(amp1,amp2);
tmin1 = min(amp1,amp2);
tmax2 = max(tmax0,tmax1);
tmin2 = min(tmin0,tmin1);

if isign == 0;

if (tmax2 > amp_fac1*tmin2) 
   n=i+1;
   [dumm,imax]=max(finamp(1:n));
   [dumm,imin]=min(finamp(1:n));
   if (imax > imin) % increase
       isign = 1; % pick minimum.
       icand = imin;
       tmax3 = finamp(imin);
       tmax0 = 0;
       tmin0 = 1;
   else
       isign = -1; % pick maximum.
       icand = imax;
       tmin3 = finamp(imax);
       tmax0 = 0;
       tmin0 = 1;
   end
i=i-1;
end

tmax0 = tmax2;
tmin0 = tmin0;

else
% Start main picking.

if (amp12 < 0) % increasing. 
   if (tmax2 > tmax3) & (isign == -1) % find another max.
      icand = i+1;
      tmax3 = finamp(icand);
      tmin2 = 1;
   end

   if (tmax2 > amp_fac1*tmin2) & (isign == 1) % find next maximum.
      finamp0(ii) = finamp(icand);
      fintime0(ii) = fintime(icand);
      finbaz0(ii) = finbaz(icand);
      idetect(ii) = iidetect(icand);
      factor(ii) = min(fac1,(tmax2/tmin2));
      fac1 = tmax2/tmin2;
      tmax3 = finamp(i+1);
      tmin3 = 0;
      tmax0 = 0;
      tmin0 = 1;
      tmax2 = 0;
      ii = ii + 1;
      isign = -1;
      icand = i + 1;
   else
      tmax0 = tmax2;
      tmin0 = tmin2;
   end

else % decreasing.

   if (tmin2 < tmin3) & (isign == 1) % find another min.
      icand = i + 1;
      tmin3 = finamp(icand);
      tmax2 = 0;
   end

   if (tmax2 > amp_fac1*tmin2) & (isign == -1) % find next minimum.
      finamp0(ii) = finamp(icand);
      fintime0(ii) = fintime(icand);
      finbaz0(ii) = finbaz(icand);
      idetect(ii) = iidetect(icand);
      factor(ii) = min(fac1,(tmax2/tmin2));
      fac1 = tmax2/tmin2;
      tmin3 = finamp(i+1);
      tmax3 = 0;
      tmax0 = 0;
      tmin0 = 1;
      tmin2 = 1;
      ii = ii + 1;
      isign = 1;
      icand = i+1;
   else
      tmax0 = tmax2;
      tmin0 = tmin2;
   end
end
end
i=i+1;
end

% End of first step.
% We picked maximum which have more than amp_fac1 times larger value than 
% local minimum.


% Second step
% Define duration for long term average.
iidur = floor(idur4lta/2);
ndur1 = -iidur:-4;
ndur2 = 4:iidur;
amp_cre_mean = median(maxamp);

mm1 = size(finamp0);
jj = 1;
for j = 1 : mm1(2)-2;
    tamp1 = finamp0(j);
    tamp2 = finamp0(j+1);
    tamp3 = finamp0(j+2);
    tamp12 = tamp1 - tamp2;
    tamp23 = tamp2 - tamp3;
    if (tamp12 < 0) & (tamp23 > 0)
       
% Need a step which can define amp_cre (lta over dur4time)
    if (idetect(j)-iidur+1 > 0) & (idetect(j)+iidur+1 < mm0(1))
        amp_cre1 = mean(maxamp(idetect(j+1)+ndur1));
        amp_cre2 = mean(maxamp(idetect(j+1)+ndur2));
        amp_cre = (amp_cre1+amp_cre2)/2*amp_fac2;

       if (tamp2 > amp_cre) | (tamp2 >= amp_cre_mean*2)
          fintime1(jj)=fintime0(j+1);
          finamp1(jj)=finamp0(j+1);
          finbaz1(jj)=finbaz0(j+1);
          jj = jj + 1;
       end
    end
    end
end

%time(idetect)


% End of second step
% Finally, maximum which have more than amp_fac2 times larger value than 
% average valude over dur4lta are selected. 

% Select picks for given julian day.
mm11 = size(fintime1);
jj = 1;
for j = 1 : mm11(2)
    if (fintime1(j) < jday+1) & (fintime1(j) >= jday)
        time21(jj)=fintime1(j);
        maxamp21(jj)=finamp1(j);
        maxbaz21(jj)=finbaz1(j);
        jj = jj + 1;
    end
end


end

clf
plot(time,maxamp)
hold on
plot(time21,maxamp21,'r.')
xlim([jday,jday+1])
