%function y = vinterp(u,stepsize,finaliv,order)
%
%   The output vector Y (VARYING) is an interpolated version
%   of the input U.  The INDEPENDENT VARIABLEs of Y start at
%   MIN(GETIV(U)) and finish at FINALIV (or MAX(GETIV(U)) if only
%   2 arguments), in increments of STEPSIZE.  The choices of
%   order are:
%       0	Zero order hold (default if only 3 args)
%       1	Linear interpolation
%
%   See also: DTRSP, TRSP, and VDCMATE.

function y = vinterp(u,stepsize,finaliv,order)

if nargin == 0,
    disp('usage: y = vinterp(u,stepsize,finaliv,order)')
    return
    end

if nargin < 2 | nargin > 4,
    error('usage: y = vinterp(u,stepsize,finaliv,order)')
    return
    end

[utype,nru,ncu,nptsu] = minfo(u);
if utype == 'syst'
    error('cannot interpolate SYSTEM matrices')
    return
elseif utype == 'vary'
    [udat,uptr,ut] = vunpck(u);
else
    ut = 0;
    udat = u;
    end

if ncu ~= 1,
    if nru == 1,
        udat = reshape(udat.',ncu*nptsu,1);
        nru = ncu;
    else
        error('input is not a vector')
            return
            end
    end

if nargin < 4,
    order = 0;
    if nargin < 3
        finaliv = max(ut);
        end
    end

if finaliv < ut(1);
    error('final iv preceeds starting iv')
    return
    end

yt = [ut(1):stepsize:finaliv].';
npts = length(yt);
teps = stepsize*1e-8;

%   initial vector first 

y = zeros(nru*npts+1,2);

%    we index on the shorter of yt or ut.  Note that this relies
%    on the command y([]) = []; leaving y unchanged.

if order == 0,
    if npts <= length(ut),
        for i = 1:npts,
            index = max(find(ut <= yt(i)+teps));
            y((i-1)*nru+1:i*nru,1) = udat((index-1)*nru+1:index*nru);
            end
    elseif length(ut) == 1,
        y(1:nru*npts,1) = reshape(udat*ones(1,npts),nru*npts,1);
    else
        for i = 1:length(ut)-1,
            index = find(yt >= ut(i) & yt < ut(i+1));
            y((min(index)-1)*nru+1:max(index)*nru,1) = ...
              reshape(udat((i-1)*nru+1:i*nru)*ones(1,length(index)), ...
                nru*length(index),1);
            end
        i = length(ut);
        index = find(yt >= ut(length(ut)));
        y((min(index)-1)*nru+1:max(index)*nru,1) = ...
          reshape(udat((i-1)*nru+1:i*nru)*ones(1,length(index)), ...
          nru*length(index),1);
        end
    y(1:length(yt),2) = yt;
    y((nru*length(yt))+1,2) = inf;
    y((nru*length(yt))+1,1) = length(yt);
elseif order == 1,
    if length(ut) == 1,
        y(1:nru*npts,1) = reshape(udat*ones(1,npts),nru*npts,1);
    elseif npts <= length(ut),
        for i = 1:length(yt);
            lix = max(find(ut <= yt(i)));
            hix = min(find(ut > yt(i)));
            if isempty(hix),
                j = [i:length(yt)];
                y((min(j)-1)*nru+1:max(j)*nru,1) = ...
                  reshape(udat((lix-1)*nru+1:lix*nru)*ones(1,length(j)), ...
                  nru*length(j),1);
                break
            else
                tscl = (yt(i) - ut(lix))/(ut(hix)-ut(lix));
                udiff = udat((hix-1)*nru+1:hix*nru)-udat((lix-1)*nru+1:lix*nru);
                y(((i-1)*nru+1):i*nru,1) = ...
                                 udat((lix-1)*nru+1:lix*nru)+ tscl*udiff;
	        end
            end
    else
        for i = 1:length(ut)-1,
            index = find(yt >= ut(i) & yt < ut(i+1));
            diffvect = udat((i)*nru+1:(i+1)*nru)-udat((i-1)*nru+1:i*nru);
            y((min(index)-1)*nru+1:max(index)*nru,1) = ...
              reshape(udat((i-1)*nru+1:i*nru)*ones(1,length(index)) ...
              + diffvect*((yt(index).'-ut(i))./(ut(i+1)-ut(i))), ...
              nru*length(index),1);
            end
        i = length(ut);
        index = find(yt >= ut(length(ut)));
        y((min(index)-1)*nru+1:max(index)*nru,1) = ...
          reshape(udat((i-1)*nru+1:i*nru)*ones(1,length(index)), ...
            nru*length(index),1);
        end
    y(1:length(yt),2) = yt;
    y((nru*length(yt))+1,2) = inf;
    y((nru*length(yt))+1,1) = length(yt);
else
    error('order must be 0 or 1')
    return
    end
%
% Copyright MUSYN INC 1991,  All Rights Reserved
