function value = rpbrk(rp,part)
%RPBRK Part(s) of a rpform.
%
%   VALUE = RPBRK(RP,PART)  returns the particular part specified by the string 
%   PART, which may be (the beginning character(s) of) one of the following 
%   strings: 
%   'breaks', 'coefs', 'pieces' or 'l', 'order' or 'k', dim'ension, 'interval'.
%
%   RPBRK(PP)  returns nothing, but prints all parts.
%
%   PJ = RPBRK(RP,J)  returns the rpform of the J-th polynomial piece of the 
%   function in RP.
%
%   PC = RPBRK(RP,[A B])  returns the restriction to the interval  [A .. B]  of 
%   the function in RP.
%
%   If RP contains an m-variate spline and PART is not a string, then it
%   must be a cell-array, of length m .
%
%   For example, if RP contains a bivariate spline with at least 4 pieces
%   in the first variable, then
%
%      rpp = rpbrk(rp,{4,[-1 1]});
%
%   gives the bivariate spline that agrees with the given one on the
%   rectangle  [rp.breaks{1}(4) .. [rp.breaks{1}(5)] x [-1 1] .
%
%   See also RPMAK, RSBRK, PPBRK, SPBRK, FNBRK.

%   Carl de Boor 23may99
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.3 $  $Date: 2000/05/09 02:27:01 $

if ~isstruct(rp)|~strcmp(rp.form([1 2]),'rp')
   error('The input does not seem to describe a function in rpform.')
end 

if nargin>1
   if nargout>1, error('Too many output arguments for the given input.')
   end
   if isstr(part)
      switch part(1)
         case 'd',       value = rp.dim;
         case {'l','p'}, value = rp.pieces;
         case 'b',       value = rp.breaks;
         case {'k','o'}, value = rp.order;
         case 'c',       value = rp.coefs;
         case 'i'
            if iscell(rp.breaks)
               for i=length(rp.order):-1:1
                  value{i} = rp.breaks{i}([1 end]); end
            else
               value = rp.breaks([1 end]);
            end
         otherwise
            error(sprintf('''%s'' is not part of a rpform.',part))
      end
   else % we are to restrict RP to some interval or piece
      pp = fn2fm(rp,'pp');
      if iscell(part)  % we are dealing with a tensor-product spline

         [breaks,c,l,k,d] = ppbrk(pp); m = length(breaks);
         sizec = [d,l.*k]; %size(c);
         if length(sizec)~=m+1, error('Information in RP is inconsistent.'),
         end
         for i=m:-1:1
            dd = prod(sizec(1:m));
            ppi = ppbrk1(ppmak(breaks{i},reshape(c,dd*l(i),k(i)),dd),...
                        part{i}) ;
            breaks{i} = ppi.breaks; sizec(m+1) = ppi.pieces*k(i);
            c = reshape(ppi.coefs,sizec);
            if m>1
               c = permute(c,[1,m+1,2:m]);
               sizec(2:m+1) = sizec([m+1,2:m]);
            end
         end
         out1 = ppmak(breaks,c,sizec);

      else  % we are dealing with a univariate spline

         out1 = ppbrk1(pp,part);
      end
   value = fn2fm(out1,'rp');
   end
else
   if nargout==0
      if iscell(rp.breaks) % we have a multivariate spline and, at present,
                           % I can't think of anything clever to do; so...
        disp(rp)
      else
        disp('breaks(1:l+1)'),        disp(rp.breaks)
        disp('coefficients((d+1)*l,k)'),  disp(rp.coefs)
        disp('pieces number l'),      disp(rp.pieces)
        disp('order k'),              disp(rp.order)
        disp('dimension d of target'),disp(rp.dim)
      end
   else
      error('More output arguments than expected.')
   end
end

function pppart = ppbrk1(pp,part)
%PPBRK1 restriction of pp to some piece or interval

[ignored, lp] = size(part);
if lp > 1 , % extract the part relevant to the interval 
            % specified by  part =: [a b]  
   pppart = ppcut(pp,part(1,1:2));
else        % extract the part-th polynomial piece of pp (if any)
   pppart = pppce(pp,part);
end

function ppcut = ppcut(pp,interv)
%PPCUT returns the part of pp  specified by the interval interv =: [a b]  

xl = interv(1); xr = interv(2); if xl>xr, xl = xr; xr = interv(1);  end
 
%  the first pol. piece is  jl ,
% the one responsible for argument  xl
jl=pp.pieces; index=find(pp.breaks(2:jl)>xl); 
                                   % note that the resulting  index  ...
if (~isempty(index)), jl=index(1); % ... is shifted down by one  ...
end                                % ... because of  breaks(2: ...
%  if xl ~= breaks(jl), recenter the pol.coeffs.
x=xl-pp.breaks(jl);
di = pp.dim;
if x ~= 0
   a=pp.coefs(di*jl+[1-di:0],:);
   for ii=pp.order:-1:2
      for i=2:ii
         a(:,i)=x*a(:,i-1)+a(:,i);
      end
   end
   pp.coefs(di*jl+[1-di:0],:)=a;
end
 
%  the last pol. piece is  jr ,
% the one responsible for argument  xr .
jr=pp.pieces;index=find(pp.breaks(2:jr+1)>=xr); 
                                   % note that the resulting ...
if (~isempty(index)), jr=index(1); % index  is shifted down by
end                                % ... one because of  breaks(2: ...
 
%  put together the cut-down  pp
di = pp.dim;
ppcut = ppmak([xl pp.breaks(jl+1:jr) xr], ...
                        pp.coefs(di*(jl-1)+[1:di*(jr-jl+1)],:),di);

function pppce = pppce(pp,j)
%PPPCE returns the j-th polynomial piece of pp  (if any).

%  if  pp  has a  j-th  piece, ...
if (0<j)&(j<=pp.pieces)  %             ...  extract it
   di = pp.dim;
   pppce = ppmak([pp.breaks(j) pp.breaks(j+1)], ...
              pp.coefs(di*j+[1-di:0],:),di);
else
   error(['The given rp function does not have that many pieces.']);
end
