function [bt_arr]=abx_track_pick(hdf_file,trk_root,setno,clims,nwindow,nterp,iiseq,iisin,seqlim)
% Call [bt_arr]=abx_track_pick(hdf_file,setno,clims,nwindow,nterp,iiseq,iisin);
% Plots bearing vs. time image and allows user to pick points off 
% multiple traces.  Outputs time-brg track file for input into abxm.
% hdf_file = input
% trk_root for '_trk' output file;  default = hdf_file(1:end-4);
% setno     = dataset number in HDF file;
% clims     = [cmin,cmax] = min,max correlation or dB level ([]=min,max in data);
% nwindow	= # pixels in bearing on each side of chosen points to look for 
%			peak (0=no quadratic interpolation);
% nterp		= # samples in time to interpolate between chosen points;
% iiseq:	1=seq #, 2=MOY, 3=HOD on vertical axis of HDF file;
% iisin:	0=degrees, 1=sin(phi) on horizontal axis of HDF file;
%
% COPYRIGHT (c) 2008
% 
% Applied Research Laboratories, The University of Texas at Austin
% (ARL:UT)
% All rights reserved.
% 
% This material may be reproduced by or for the U.S. Government
% only. This notice must appear in all copies of this file and its
% derivatives. Further distribution outside the U.S. Government is
% prohibited without the explicit written authorization from the
% Director of the Environmental Sciences Laboratory, ARL:UT.
% 
% THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND,
% EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE
% IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
% PURPOSE, OR NONINFRINGEMENT.
% 
% IN NO EVENT WILL ARL:UT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING
% OUT OF THE USE OR INABILITY TO USE THIS SOFTWARE, EVEN IF ARL:UT HAS
% BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.

if isempty(trk_root)
    trk_root = hdf_file(1:end-4);
end

%tms:  use matlab hdfsd interface
[SD,status]=hdfsd_open(hdf_file);
dim_lims={[],seqlim};
[DATA,status]=hdfsd_read(SD,dim_lims,setno);
time=SD.scales{1}{3};
brg=SD.scales{1}{1};
data=squeeze(DATA.data)';

%keyboard

%time=timex(seqlim(1):seqlim(2));
%data=datax(seqlim(1):seqlim(2),:);
[ntime,nbrg]=size(data);
tfac=(time(ntime)-time(1))/(ntime-1);
bfac=(brg(nbrg)-brg(1))/(nbrg-1);

%
if isempty(clims), 
	dmin=min(min(data)); dmax=max(max(data)); clims=[dmin dmax]; 
end
figure; imagesc(brg,time,data,clims); 
colormap(1-gray);
%colormap('default');
grid;
xlabel('Bearing'); ylabel('Time');
plot_set3(gca,1,'Times',12,0);
%
disp('Click on (time,brg) points on traces.  <cr> to complete the trace')

[b_in,t_in] = ginput;
npt_in=length(t_in);
if npt_in==0, break; end;
%
% Interpolate in time:
if nterp>0,
	npt=(npt_in-1)*(nterp+1) + 1;
	t=interp1([1:npt_in]',t_in,linspace(1,npt_in,npt)','linear');
   b=interp1(t_in,b_in,t,'linear');
else
	t=t_in; b=b_in; npt=npt_in;
end

bt_arr=[]; t_uni=[];
for j=1:npt,
	tj=t(j); bj=b(j);
	iit=round((tj - time(1))/tfac + 1);
	iib=round((bj - brg(1))/bfac + 1);
%
	iibb=[max(1,iib-nwindow):min(nbrg,iib+nwindow)]; nbb=length(iibb);
	bwindow=brg(iibb); twindow=time(iit);
	dwindow=data(iit,iibb);
%	bd=[bwindow;dwindow],
	[d_max,iibrg]=max(dwindow);
%
% Note: use quad_max_3pt here:
	if nwindow==0 | iibrg==1 | iibrg==nbb, 
% Don't try to fit quadratic to max at endpoint or else maximum will be outside window:
		brg_fit_max=bwindow(iibrg);	
	else,
		iicubfit=iibrg-1:iibrg+1;
		abc=polyfit(bwindow(iicubfit),dwindow(iicubfit),2);
		brg_fit_max=-abc(2)/(2*abc(1));
	end
	bt_arr=[bt_arr; time(iit), brg_fit_max];
	t_uni=[t_uni; time(1) + tfac*(iit-1)];
end
%
% Overlay points on image and allow user to delete any:
hold on; ok=0;
while ok==0;
	nptx=size(bt_arr,1);
	h1=plot_cfill2(bt_arr(:,2),t_uni,'o','r',4);
	h2=text(bt_arr(:,2),t_uni,int2str([1:nptx]'));
	iibad=input('Enter points to delete, <CR> for none: ');
	if any(iibad>nptx) | any(iibad<0), 
		disp('Illegal point, try again'); 
		delete(h1); delete(h2);
	else,
		if isempty(iibad), 
			ok=1; 
		else
			iiok=1:nptx; iiok(iibad)=zeros(1,length(iibad));
			iiok=find(iiok~=0);
			bt_arr=bt_arr(iiok,:); t_uni=t_uni(iiok);
			delete(h1); delete(h2);
		end
	end
end	

%
eol=[char(13) char(10)];
fid = fopen([trk_root '_trk'],'wt');
fprintf(fid,'$$$ TEMPLATE FOR STEERING ANGLE FILE INPUT TO ABX WITH IIGRAM=2 $$$\n');
fprintf(fid,'$$$ NOTE: DATA WILL BE READ AFTER LINES BEGINNING WITH "*" $$$\n');
fprintf(fid,'------------------------------------------------------------------\n');
fprintf(fid,'NTRACK   ==> Number of tracks to form\n');
fprintf(fid,'IISEQ: 1 ==> Time axis given by sequence number of FFT file\n');
fprintf(fid,'       2 ==> Time axis given by MOY relative to first MOY in FFT file\n');
fprintf(fid,'       3 ==> Time axis given by Hour of Day\n');
fprintf(fid,'IISIN: 0 ==> Azimuth angles given in deg E of N\n');
fprintf(fid,'       1 ==> Azimuth angles given in sin(phi) rel LOB\n');
fprintf(fid, '*(1) # TRACKS; IISEQ;  IISIN\n');
fprintf(fid,'            1      %1i       %1i\n',iiseq,iisin);
fprintf(fid,'*(2) SEQ#,MOY, or HOD; NTRACK AZIMUTHS\n');
fprintf(fid,'%10.3f %10.3f\n',bt_arr');
fclose(fid)

disp(['Output in file ' trk_root '_trk'])
