home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / VINTERP.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.1 KB  |  137 lines
  1. %function y = vinterp(u,stepsize,finaliv,order)
  2. %
  3. %   The output vector Y (VARYING) is an interpolated version
  4. %   of the input U.  The INDEPENDENT VARIABLEs of Y start at
  5. %   MIN(GETIV(U)) and finish at FINALIV (or MAX(GETIV(U)) if only
  6. %   2 arguments), in increments of STEPSIZE.  The choices of
  7. %   order are:
  8. %       0    Zero order hold (default if only 3 args)
  9. %       1    Linear interpolation
  10. %
  11. %   See also: DTRSP, TRSP, and VDCMATE.
  12.  
  13. function y = vinterp(u,stepsize,finaliv,order)
  14.  
  15. if nargin == 0,
  16.     disp('usage: y = vinterp(u,stepsize,finaliv,order)')
  17.     return
  18.     end
  19.  
  20. if nargin < 2 | nargin > 4,
  21.     error('usage: y = vinterp(u,stepsize,finaliv,order)')
  22.     return
  23.     end
  24.  
  25. [utype,nru,ncu,nptsu] = minfo(u);
  26. if utype == 'syst'
  27.     error('cannot interpolate SYSTEM matrices')
  28.     return
  29. elseif utype == 'vary'
  30.     [udat,uptr,ut] = vunpck(u);
  31. else
  32.     ut = 0;
  33.     udat = u;
  34.     end
  35.  
  36. if ncu ~= 1,
  37.     if nru == 1,
  38.         udat = reshape(udat.',ncu*nptsu,1);
  39.         nru = ncu;
  40.     else
  41.         error('input is not a vector')
  42.             return
  43.             end
  44.     end
  45.  
  46. if nargin < 4,
  47.     order = 0;
  48.     if nargin < 3
  49.         finaliv = max(ut);
  50.         end
  51.     end
  52.  
  53. if finaliv < ut(1);
  54.     error('final iv preceeds starting iv')
  55.     return
  56.     end
  57.  
  58. yt = [ut(1):stepsize:finaliv].';
  59. npts = length(yt);
  60. teps = stepsize*1e-8;
  61.  
  62. %   initial vector first 
  63.  
  64. y = zeros(nru*npts+1,2);
  65.  
  66. %    we index on the shorter of yt or ut.  Note that this relies
  67. %    on the command y([]) = []; leaving y unchanged.
  68.  
  69. if order == 0,
  70.     if npts <= length(ut),
  71.         for i = 1:npts,
  72.             index = max(find(ut <= yt(i)+teps));
  73.             y((i-1)*nru+1:i*nru,1) = udat((index-1)*nru+1:index*nru);
  74.             end
  75.     elseif length(ut) == 1,
  76.         y(1:nru*npts,1) = reshape(udat*ones(1,npts),nru*npts,1);
  77.     else
  78.         for i = 1:length(ut)-1,
  79.             index = find(yt >= ut(i) & yt < ut(i+1));
  80.             y((min(index)-1)*nru+1:max(index)*nru,1) = ...
  81.               reshape(udat((i-1)*nru+1:i*nru)*ones(1,length(index)), ...
  82.                 nru*length(index),1);
  83.             end
  84.         i = length(ut);
  85.         index = find(yt >= ut(length(ut)));
  86.         y((min(index)-1)*nru+1:max(index)*nru,1) = ...
  87.           reshape(udat((i-1)*nru+1:i*nru)*ones(1,length(index)), ...
  88.           nru*length(index),1);
  89.         end
  90.     y(1:length(yt),2) = yt;
  91.     y((nru*length(yt))+1,2) = inf;
  92.     y((nru*length(yt))+1,1) = length(yt);
  93. elseif order == 1,
  94.     if length(ut) == 1,
  95.         y(1:nru*npts,1) = reshape(udat*ones(1,npts),nru*npts,1);
  96.     elseif npts <= length(ut),
  97.         for i = 1:length(yt);
  98.             lix = max(find(ut <= yt(i)));
  99.             hix = min(find(ut > yt(i)));
  100.             if isempty(hix),
  101.                 j = [i:length(yt)];
  102.                 y((min(j)-1)*nru+1:max(j)*nru,1) = ...
  103.                   reshape(udat((lix-1)*nru+1:lix*nru)*ones(1,length(j)), ...
  104.                   nru*length(j),1);
  105.                 break
  106.             else
  107.                 tscl = (yt(i) - ut(lix))/(ut(hix)-ut(lix));
  108.                 udiff = udat((hix-1)*nru+1:hix*nru)-udat((lix-1)*nru+1:lix*nru);
  109.                 y(((i-1)*nru+1):i*nru,1) = ...
  110.                                  udat((lix-1)*nru+1:lix*nru)+ tscl*udiff;
  111.             end
  112.             end
  113.     else
  114.         for i = 1:length(ut)-1,
  115.             index = find(yt >= ut(i) & yt < ut(i+1));
  116.             diffvect = udat((i)*nru+1:(i+1)*nru)-udat((i-1)*nru+1:i*nru);
  117.             y((min(index)-1)*nru+1:max(index)*nru,1) = ...
  118.               reshape(udat((i-1)*nru+1:i*nru)*ones(1,length(index)) ...
  119.               + diffvect*((yt(index).'-ut(i))./(ut(i+1)-ut(i))), ...
  120.               nru*length(index),1);
  121.             end
  122.         i = length(ut);
  123.         index = find(yt >= ut(length(ut)));
  124.         y((min(index)-1)*nru+1:max(index)*nru,1) = ...
  125.           reshape(udat((i-1)*nru+1:i*nru)*ones(1,length(index)), ...
  126.             nru*length(index),1);
  127.         end
  128.     y(1:length(yt),2) = yt;
  129.     y((nru*length(yt))+1,2) = inf;
  130.     y((nru*length(yt))+1,1) = length(yt);
  131. else
  132.     error('order must be 0 or 1')
  133.     return
  134.     end
  135. %
  136. % Copyright MUSYN INC 1991,  All Rights Reserved
  137.