function [data_oi,pos,u_alt,u_ang] =  ...
    oi_top_track(file,dataregion,lonscl,latscl,timescl,err_var,deci,date_anal)
% Alongtrack optimal interpolation of topex altimeter data 
%
% [data_oi,pos,u_alt,u_ang] =  ...
%    oi_top_track(file,dataregion,lonscl,latscl,timescl,err_var,deci,date_anal)
%
% Inputs:
%
% FILE 		 input netcdf file of track data
% DATAREGION 	 [lonW lonE latS latN] geographic region of input data
% LONSCL 	 N-S scale of correlation in degrees
% LATSCL 	 E-W scale of correlation in degrees
% TIMESCL 	 time scale of correlation in days
% ERR_VAR	 assumed error variance in the track data 
% DECI 		 factor by which to decimate along track data
% DATE_ANAL 	 required date at which to produce OI
%
% Outputs:
%
% DATA_OI	 OI smoothed data (length N)
% POS		 position (lat lon time) of data_oi (N X 3)
% U_ALT,U_ANG    magnitude and direction of crosstrack velocity (by
%                geostrophy) corresponding to OI smoothed alongtrack height 
%                data (for subsequent comparison with other velocities)
%
% EXAMPLE:
% [data_oi,pos] =  ...
%    oi_top_trk('/home/wilkin/altim_map/top_93_tas',[150 162 -40 -28],...
%       0.5,0.5,15,0.01,3,[1993 3 13 0 0 0]);
%
% Alison Walker 27 July 97
%
% John Wilkin 29 July 97 
%    mods to use direct time data (3-d) returned by get_topex_3d
%    (should probably use option in get_topex_3d to get 6 cycles)

%  GET INPUT DATA FROM TOPEX NETCDF FILE
disp('')
disp(['Reading from altimeter netcdf file ' file '.nc'])
numcyc = 2;
[lon,lat,time,height_diff] = get_topex_3d(file,date_anal,dataregion,numcyc);

bell;disp(['Getting ' int2str(numcyc) ' cycles of Topex data'])

% FIND THE TEMPORAL LAG OF THE TRACKS
day_diffs = julian([1985 1 1 0 0 0]) + time/86400 - julian(date_anal);

% SETUP PARAMETERS FOR ALONG TRACK OI
% it would be advisable to narrow the lengthscl
% as we are only interested in a very local smoothing
% preferably the fn should be reasonably sharp around zero lag

% REGION LIMITS TO BE PASSED TO OI ROUTINE
smregions = [dataregion 0];

% INITIALIZE THE OUTPUTS 
data_oi = [];
pos     = [];
u_alt   = [];
u_ang   = [];

num_cycles = 1;

% FOR EACH TRACK COMBINE ALL AVAILABLE CYCLES
disp([' Optimally interpolating alongtrack data'])
disp(['  Each track and cycle is smoothed separately, and decimated by a'])
disp(['   factor of ' int2str(deci)])
for k = 1:size(height_diff,3)
  disp([' Cycle ' num2str(k) ' ...'])
  for j = 1:size(height_diff,2)
    disp(['  Track ' num2str(j) ' ...'])
    h = height_diff(:,j,k);
    t = day_diffs(:,j,k);
    
    % find only those obs that are not NaN
    obs = find(isnan(h)~=1);
    h = h(obs);
    
    if length(h) > 10 % define some minimum number of data pts necessary for OI
      
      % duplicate columns of lon/lat coordinates
      llat = lat(:,j); 
      llat = llat(:,ones([1 num_cycles]));
      llon = lon(:,j); 
      llon = llon(:,ones([1 num_cycles]));
      tt   = t;
      
      % clip out the coordinates where the height data was bad
      llat = llat(obs);
      llon = llon(obs);
      tt   = tt(obs);
      posdata = [llat(:) llon(:) tt(:)];
      
      % The analysis positions are the input data positions decimated by a
      % factor deci.
      % The three coordinates are carried through to the output.
      % Combining different observation times to a single analysis date is
      % left to the 2-d gridding optimal interpolation step 
      deciposdata = posdata;
      smsub = 1:deci:length(h);
      deciposdata = deciposdata(smsub,:);
      
      % DO THE OI ON EACH TRACK
      % Note that some parts of the tracks lie outside the dataregion.
      % These points are still part of DECIPOSDATA and turn up as
      % zeros in smdata, but they are subsequently clipped out and the end
      % of this loop      
      smdata = multi_oi_err(h,deciposdata,posdata,...
	  smregions,lonscl,latscl,timescl,err_var);
      
      % prepare to clip data that are outside the requested region
      inreg = findinrange(deciposdata(:,[2 1]),dataregion);
      
      if nargout > 2
	
	% compute magnitude and heading of cross-track velocity component
	% for later comparison with gridded velocities
	[ualt,angle] = ...
	    alt_vel_from_hgt(deciposdata(:,1),deciposdata(:,2),smdata);
	u_alt = [u_alt; ualt(inreg)];
	u_ang = [u_ang; angle(inreg)];
      end
      data_oi = [data_oi; smdata(inreg)];
      pos     = [pos; deciposdata(inreg,:)];

    end
  end
end
