function    [DD,W] = d3echotool2_pt(recdir,prefix,cue,CL,EXTENT,DD)
%
%     DD=d3echotool(recdir,prefix,cue,CL)
%     or
%     DD=d3echotool(recdir,prefix,cue,CL,EXTENT)
%     or
%     DD=d3echotool(recdir,prefix,cue,CL,EXTENT,DD)
%
%     Inputs:
%     recdir is the deployment directory e.g., 'e:/eg15/eg15_207a'.
%     prefix is the base part of the name of the files to analyse e.g., 
%        if the files have names like 'eg207a001.wav', put prefix='eg207a'.
%     cue is the time in seconds-since-tag-on to start working from.
%        If cue = [start_cue end_cue], the length of that interval will
%        be used as the frame length instead of the default 20s.
%     CL is a vector of click cues in seconds.
%     EXTENT = [left right] are the times in seconds to display with
%        respect to each click. Default is [-0.0005 0.02].
%     DD is the data structure defined below. DD can be passed as an
%        input argument to allow multiple work sessions.
%
%  	Each matrix in the cell array DD contains the following columns:
%  	1  click_cue
%  	2  two_way_travel_time to selected echo, s
%  	3  echo envelope level, dB
%
%	  Valid commands in the figure window are:
%		left-button		select an echo
%     +  extend the two-way travel time extent
%     -  reduce the two-way travel time extent
%     >  move the two-way travel time extent to the right (1/3 of window length)
%     <  move the two-way travel time extent to the left (1/3 of window length)
%     a  accept the current echo sequence
%     e  edit the echo sequence closest to the crosshairs
%     f  move forward to the next block of clicks
%     b  move backward to the previous block of clicks
%     k decrease low caxis
%     l increase low caxis
%     u decrease high caxis
%     U increase high caxis
%     2 move twtt to between 90 and 200 m
%     3 move twtt to between 190 and 300 m
%     i move twtt 50 m to the right
%     o move twwt 50 m to the left
%	  r  remove the current echo sequence
%	  q  quit

%  markjohnson@st-andrews.ac.uk
%  Sept 2018
%

if nargin<4,
    help d3echotool2
    return
end

nshow = 0.0015 ;           % length of waveform to display in detail figure
figs = [1 2] ;             % which figure windows to use
MAXICI = 0.4 ;               % draw a horizontal black line when ICI is more than this
SRCH = 0.0025 ;            % time window over which to search for z command
TH = 0.7 ;
DF = 8 ;
nCL = 70 ;                 % default is to read in 70 clicks at a time
CH = 1 ;
W = [] ;
maxENR= 50;

% the parameters below are:
%   cax = display colour limits in dB
%   f_env_ls = envelope detector signal bandpass filter cut-off frequencies in Hz
%   f_env_hs = same as f_env_ls but these values are used if the sampling
%      rate is > 100kHz.

switch prefix(1:2),
   case 'zc',      % for ziphius use:
      cax = [-93 -35] ;         
      f_env_ls = [25000 45000] ;
      f_env_hs = [30000 65000] ;

   case 'md',      % for mesoplodon use:
      cax = [-91 -35] ;       
      f_env_ls = [25000 45000] ;
      f_env_hs = [25000 65000] ;

   case {'pw','gm'}      % for pilot whale use:
      cax = [-95 -30] ;      
      f_env_ls = [20000 40000] ;
      f_env_hs = [30000 65000] ;

   case {'hp','pp'}      % for hartbour porpoise use:
      cax = [-103 -30] ;      
      f_env_ls = 40000 ;
      f_env_hs = [100000 230000] ;

   case 'by',
      cax = [-95 -10] ;
      f_env_ls = 20000 ;
      f_env_hs = [20000 55000] ; 
      
   case {'sw','pm'},      % for sperm whale use:
      cax = [-93 -15] ;
      f_env_vls = [2000 15000] ;
      f_env_ls = [2000 30000] ;
      f_env_hs = [5000 40000] ;

   otherwise,      % for others use:
      cax = [-95 -10] ;
      f_env_ls = [2000 35000] ;
      f_env_hs = [2000 35000] ;
end

cax = [3 maxENR];

% handle variable arguments
if nargin<5 || isempty(EXTENT),
   EXTENT = [-0.0005 0.02] ;
end
if length(EXTENT)==1,
   EXTENT = [-0.0005 EXTENT] ;
end

if nargin<6,
   DD = {} ;       % initialize cell array of echo sequences
end

if isempty(cue),
   cue = CL(1)-0.1 ;
end

if length(cue)==1
    cue = max(CL(1,1),cue) ;
elseif  length(cue)==2,
    len = diff(cue);
    nCL = [] ;
%     cue = cue(1) ;
end

CL = CL(:,1) ;

% find sampling rate
[x fs] = d3wavread(cue+[0 0.1],recdir,prefix,'wav') ;
if isempty(x), return, end

% make analysis filter
if fs<50e3,
   if exist(f_env_vls,'var'),
      [b a] = butter(4,f_env_vls/(fs/2)) ;
   else
      [b a] = butter(4,f_env_ls(1)/(fs/2),'high') ;
   end
elseif fs<100e3,
   if length(f_env_ls)==1,
      [b a] = butter(4,f_env_ls/(fs/2),'high') ;
   else
      [b a] = butter(4,f_env_ls/(fs/2)) ;
   end
else
   [b a] = butter(4,f_env_hs/(fs/2)) ;
end

% initialize
fs = fs/DF ;
currentseq = length(DD)+1 ;
seqplot = NaN*ones(length(DD),1) ;
figure(figs(1))
clf
next = 1 ; 
off_frame = [] ;
MODES = {'BOTTOM','TARGET'} ;
MODE = 1 ;
cleanh = [] ;

while next<2,
   while next==1,
      fprintf('reading at %d  ', cue) ;
      if ~isempty(nCL),
          k = find(CL>=cue,nCL) ;
          len = min(CL(k(end))-CL(k(1))+1,120)
      elseif isempty(nCL),
          k = find(CL>=cue(1) & CL<=cue(2)) ;
          cue = cue(1);
      end
          
      x = d3wavread(cue+[0 len+EXTENT(2)+0.01],recdir,prefix,'wav') ;
      if isempty(x), return; end

      overlap = 0.1*len ;  % 10% overlap between frames
      
      if ~isempty(k),
         cl = CL(k) ;
         fprintf('%d clicks\n', length(cl)) ;       
      else
         cl = [] ;
      end
      
      if length(cl)>3,
         x = hilbenv(filter(b,a,x(:,CH)));
         xe = mean(buffer(x.^2,2*DF,DF,'nodelay'));
         R = extract_cues(xe,fs,cl(:,1)-cue+EXTENT(1),diff(EXTENT));
         hold off, zoom off
         tk = (0:size(R,1)-1)/fs+EXTENT(1) ;
         NL = 10*log10(prctile(R(:),5));
         imagesc(tk,1:length(cl),adjust2Axis(10*log10(R)'-NL),cax); grid
         hold on
         colormap(jet) ;
         kk = find(diff(cl)>MAXICI) ;     % plot horizontal bars where ICI is large
         for kkk=kk',
            plot(EXTENT,(kkk+0.5)*[1;1],'k')
         end
         next = 0 ; 
      else
         cue = cue + len ;
         clf, plot(0,0); title('No clicks in block') ;
         drawnow ; 
         fprintf('Goto next block or quit (q=quit)? ')
         [gx gy button] = ginput(1) ;
         fprintf('\n') ;
         if button=='q',
            return ;
         end
      end
   end

   % plot old sequences in this block
   % find points already selected in this block
   for k=1:length(DD),
      ddd = DD{k} ;
      if ~isempty(ddd) & any(ddd(:,1)>=cue & ddd(:,1)<(cue+len+1)),
         kbot = nearest(cl,ddd(:,1),0.001) ;
         kk = find(~isnan(kbot)) ;
         if ~isempty(kk)
         seqplot(k) = plot(ddd(kk,2),kbot(kk),'k.-') ;
         end
      end
   end
   
   dd = NaN*ones(length(cl),2) ;
   hd = plot(dd(:,1),1:size(dd,1),'k*') ;
%    title(sprintf('%d to %d, caxis %d to %d dB',cue,cue+len,cax(1),cax(2))) ;    % Changed by Pernille 
   title(sprintf('caxis %d to %d dB',cax(1),cax(2))) ;

   while next==0,
      accept = 0 ;
      [xr yr button] = ginput(1) ;
      if isempty(button),continue,end
      switch button
      case 'a',               % type 'a' to accept the current sequence
         accept = 1 ;
      case 'k',
         cax(1) = cax(1)-3 ;
         caxis(cax) ;
         title(sprintf('caxis %d to %d dB',cax(1),cax(2))) ;    % Chenged by Pernille
      case 'l',
         cax(1) = cax(1)+3 ;
         caxis(cax) ;
         title(sprintf('caxis %d to %d dB',cax(1),cax(2))) ;    % Chenged by Pernille
      case 'u',
         cax(2) = cax(2)-3 ;
         caxis(cax) ;
         title(sprintf('caxis %d to %d dB',cax(1),cax(2))) ;    % Chenged by Pernille
      case 'U',
         cax(2) = cax(2)+3 ;
         caxis(cax) ;
         title(sprintf('caxis %d to %d dB',cax(1),cax(2))) ;    % Chenged by Pernille      
      case 'q',           % type 'q' to end session
         accept = 1 ;
         next = 2 ;
         if nargout==2,
            W.R = R ;
            W.cl = cl ;
            W.t = tk ;
         end
      case 'b',           % type 'b' to go to previous frame
         accept = 1 ;
         next = 1 ;
         cue = cue-len+overlap ;    % advance time cursor
      case 'f',           % type 'f' to go to next frame
         accept = 1 ;
         next = 1 ;
         cue = cue+len-overlap ;    % advance time cursor
      case 'r',           % type 'r' to clear all currently selected points
         dd = NaN*dd ;
         set(hd,'XData',dd(:,1)) ;
      case '+',
         EXTENT(2) = EXTENT(1)+min(400/750,diff(EXTENT)*2) ;
         accept = 1 ;
         next = 1 ;
      case '-',
         EXTENT(2) = EXTENT(1)+max(2/750,diff(EXTENT)/2) ;
         accept = 1 ;
         next = 1 ;
      case '>',
         EXTENT = EXTENT+diff(EXTENT)/3 ;
         accept = 1 ;
         next = 1 ;
      case '<',
         EXTENT = EXTENT-diff(EXTENT)/3 ;
         accept = 1 ;
         next = 1 ;
      case '2',
         EXTENT = [90/750 200/750] ;
         accept = 1 ;
         next = 1 ;
      case '3',
         EXTENT = [190/750 300/750] ;
         accept = 1 ;
         next = 1 ;
      case '4',
         EXTENT = [290/750 400/750] ;
         accept = 1 ;
         next = 1 ;
      case '5',
         EXTENT = [390/750 500/750] ;
         accept = 1 ;
         next = 1 ;
      case '6',
         EXTENT = [490/750 600/750] ;
         accept = 1 ;
         next = 1 ;
      case 'i',
         EXTENT = EXTENT-50/750 ;
         accept = 1 ;
         next = 1 ;
      case 'o',
         EXTENT = EXTENT+50/750 ;
         accept = 1 ;
         next = 1 ;
      case 'e',           % type 'e' to edit a sequence
         % look for nearest sequence to crosshairs
         yr = round(max([1 min([length(dd) yr])])) ;
         arat = abs(diff(EXTENT)/(cl(end)-cl(1))) ;
         ddist = NaN*ones(length(DD),1) ;
         for k=1:length(DD),
            ddd = DD{k} ;
            if ~isempty(ddd) & any(ddd(:,1)>=cue & ddd(:,1)<(cue+len+1)),
               ddist(k) = min(arat*abs(ddd(:,1)-cl(yr))+abs(ddd(:,2)-xr)) ;
            end
         end
         [m currentseq] = min(ddist) ;
         ddd = DD{currentseq} ;
         kbot = nearest(cl,ddd(:,1),1e-3) ;
         kk = find(~isnan(kbot)) ;
         dd(kbot(kk),:) = ddd(kk,2:end) ;
         set(seqplot(currentseq),'XData',NaN*kk) ;
         set(hd,'XData',dd(:,1)) ;
         kk = find(ddd(:,1)<cl(1)-1e-3 | ddd(:,1)>cl(end)+1e-3) ;
         off_frame = ddd(kk,:) ;

      case 1 
         if yr>0.5 & yr<length(cl)+0.5 & xr>EXTENT(1) & xr<EXTENT(2),
            kcl = round(yr) ;
            figure(figs(2)), clf
            dd(kcl,:) = showpiece(R(:,kcl),tk,xr,nshow*fs,fs) ;
            figure(figs(1))
            set(hd,'XData',dd(:,1)) ;
            newseq = 0 ;
         end
      end
      
      if accept,           % finish sequence and store in DD
         kcl = find(~isnan(dd(:,1))) ;
         if ~isempty(kcl),
            seqplot(currentseq) = plot(dd(kcl,1),kcl,'k.-') ;
            V = [off_frame;cl(kcl) dd(kcl,:)] ;
            [mm I] = sort(V(:,1)) ;
            DD{currentseq} = V(I,:) ;
            off_frame = [] ;
         end
         save echotool_Recover DD
         dd = NaN*dd ;
         currentseq = length(DD)+1 ;
         set(hd,'XData',dd(:,1)) ;
         newseq = 0 ;
      end
   end
end


function   X = showpiece(r,tk,xr,nshow,fs)
%
%
kns = round(nshow) ;
kkx = find(tk>xr,1)+(-kns:kns) ;
kkx = kkx(find(kkx>0 & kkx<length(tk))) ;
tk = tk(kkx) ;
r = 10*log10(r(kkx)) ;
plot(tk,r),grid
hold on

% plot default peak selection
[xpk m] = max(r) ;
hb = plot([1;1]*tk(m),get(gca,'YLim')'*[1 1],'k-') ;
[xr yr button] = ginput(1) ;

if button=='x'         
   X = NaN*[1 1] ;
   return
end

if button=='s' & xr>=tk(1) & xr<=tk(end)
   m = round(interp1(tk,(1:length(tk))',xr)) ;
   set(hb(1),'XData',[1;1]*tk(m)) ;
   xpk = r(m) ;
   pause(0.5)
end
   
X = [tk(m) xpk] ;
return

