function [Ys, iy] = stretch(y, nsegmed, nsegmin, nsegmax, dyfrac)

% Given an unevenly spaced grid, y, decimate it so that
% successive y increments change as slowly and smoothly
% as possible.  Ys is the new (longer) vector of y grid
% values, including original points and additional samples.
% iy is a vector of indices: Ys(iy) = y;
%
% Feed 1:length(Ys); transform input yi to iy;
% Then feed contour Ys.


%y = [1:0.5:4 4.2 4.4 4.5:20];
%nsegmed = 5;
%nsegmin = 2;
%nsegmax = 6;
% dyfrac = 0.2

[nr,nc] = size(y);
if nc == 1,
   y = y.';  % Make it a row.
   wascolumn = 1;
else
   wascolumn = 0;
end
[nr,ny] = size(y);
if (nr ~= 1)
   error('Input to stretch must be a vector, not a matrix.')
end

dy = diff(y);
dymed = median(dy);

dymedseg = dymed/nsegmed;

ns = round(dy/dymedseg);
ii = find(ns > nsegmax);
ns(ii) = nsegmax;
ii = find(ns < nsegmin);
ns(ii) = nsegmin;


iterations = 0;
nchanged = 1;
while (iterations < ny & nchanged > 0)
   dyseg = dy./ns;
   mdyseg = 0.5*(dyseg(1:(end-1)) + dyseg(2:end));
   ddyseg = diff(dyseg);
   incr_mask = [0 (ddyseg./mdyseg > dyfrac)];
   decr_mask = [(ddyseg./mdyseg < -dyfrac) 0];
   ii = find(incr_mask | decr_mask);
   nchanged = length(ii);
   if nchanged > 0, ns(ii) = ns(ii) + 1; end
   iterations = iterations + 1;
end

iy = cumsum([1 ns]);  % indices into ys of the original points.
npts = iy(end);
iys = 1:npts;

A = diag(6*ones(npts,1));
v1 = ones(npts-1,1);
v2 = ones(npts-2,1);
A = A - 4*(diag(v1,1) + diag(v1,-1)) + diag(v2,2) + diag(v2,-2);
A(2,[1 2]) = [-2 5];
A(end-1,[end-1 end]) = [5 -2];

Yf = zeros(npts, 1);  % column
Yf(iy) = y(:);       % Fixed values.

fmask = zeros(npts,1);
fmask(iy) = ones(size(iy));   % Ones at the fixed points.
vmask = ~fmask;               % Ones at the variable points.
ivar = find(vmask);

D = A * Yf;

Ar = A(ivar,ivar);         % Reduced A.
Dr = D(ivar);

% Ar * Yv = -Dr;
Yv = Ar\(-Dr);

Ys = Yf;
Ys(ivar) = Yv;

if ~wascolumn,
   Ys = Ys(:).';
end

