% set this range
%irange = 1:5:501;
%irange = 1:2:502;  % 4-m intervals over top 1000 m

irange = 1:10:(0.5*floor(depth(2*ndepth,y(1))-10)); % full depth range
use_depth = 1;     %flag

% Set this flag to 1 for an east-west section, 0 for north-south:
x_is_longitude = 0;


% The following finds the indices for the pressures corresponding
% to the even depth grid set above by irange.  This new grid then
% replaces the evenly spaced irange.
if use_depth == 1,
   d = (irange -1) *2;
   dd = d(ones(nsta,1),:)'; 
   irange = 1 + fix(press(d)/2);
end

nd = length(irange);
pp = p(ones(nsta,1),irange)';
s = S(irange,:);
ibad = find(isnan(s));
if ~isempty(ibad),
   s(ibad) = 1.e38*ones(size(ibad));
end
t = T(irange,:);
ibad = find(isnan(t));
if ~isempty(ibad),
   t(ibad) = 1.e38*ones(size(ibad));
end


% change to y for latitude, x for longitude
if x_is_longitude == 1,
   xx = x(:,ones(1,nd))';
else
   xx = y(:,ones(1,nd))';
end
[sva, sig] = svan(s(:),t(:),pp(:));
if use_depth == 1,
   xptscon = [xx(:) -dd(:) t(:) s(:) theta(s(:),t(:),pp(:),0) sva];
else
   xptscon = [xx(:) -pp(:) t(:) s(:) theta(s(:),t(:),pp(:),0) sva];
end
if ~isempty(ibad),
   xptscon(ibad,5:6) = 1.e38*ones(length(ibad), 2);
end

blankpol = mkpol(xx(1,:), -pp(:,1), t);


if use_depth == 1
   eval(['save6 ' sprintf('xdts%.0f_%.0f.con', jsta(1), jsta(nsta))...
          ' xptscon -ascii']);
   eval(['save6 ' sprintf('xdts%.0f_%.0f.pol', jsta(1), jsta(nsta))...
          ' blankpol -ascii']);

else
   eval(['save6 ' sprintf('xpts%.0f_%.0f.con', jsta(1), jsta(nsta))...
          ' xptscon -ascii']);
   eval(['save6 ' sprintf('xpts%.0f_%.0f.pol', jsta(1), jsta(nsta))...
          ' blankpol -ascii']);

end


