function jumps = fnjmp(f,x)
%FNJMP Jumps, i.e., f(x+) - f(x-) .
%
%   JUMPS = FNJMP(F,X)
%
%   returns the matrix JUMPS containing the jump  f(X+) - f(X-)  at
%   the point(s) in X in the function  f  specified by F.
%
%   For example,  fnjmp( ppmak( 1:4 , 1:3), 1:4 )  returns the vector
%   [ 0, 1, 1, 0]  (since the pp function here is 1 on [1 .. 2], 2 on [2 .. 3]
%   and 3 on [3 .. 4], hence has zero jump at 1 and 4 and a jump of 1 across
%   both 2 and 3).

%   Carl de Boor 20 feb 94
%   cb : 9 dec 95 (correct indexing mistake in pp treatment)
%   29aug96 (include BB-form) 06sep95 (correct for vector-valued  f)
%   2nov97 (improve the help), 06may98 (accept older versions of the forms)
%   3sep99 (use the string 'left' as the third argument for ppual)
%   7jan00 (also handle rational splines)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.9 $

if ~isstruct(f), f = fn2fm(f); end

if length(f.order)>1
   error('FNJMP only works with univariate functions.'), end

%  If necessary, sort  x :
[mx,nx] = size(x); lx = mx*nx; xs = reshape(x,1,lx);
jumps = zeros(fnbrk(f,'dim'),lx);

tosort = 0;
if any(diff(xs)<0)
   tosort = 1; [xs,ix] = sort(xs);
end

form = f.form;
switch form
case {'B-','BB','rB'}
   if form(1)=='r', f = fn2fm(f,'B-'); end
   [t,a,n,k,d] = spbrk(f);

   indexr = sorted(t,xs); indexl = n+k - sorted(-t,-xs);
   indexl = indexl(lx:-1:1);
   index = find(indexr-indexl>=k);
   if ~isempty(index)
      if form(1)=='B'
         a = [zeros(d,1) a zeros(d,1)];
         jumps(:,index) = a(:,indexr(index)+2-k) - a(:,indexl(index)+1);
      else
         a = [zeros(d,1) a zeros(d,1)];
         d = d-1;
         jumps(:,index) = ratjump(a(:,indexr(index)+2-k),a(:,indexl(index)+1));
      end
   end
case {'pp','rp'}
   if form(1)=='r', f = fn2fm(f,'pp'); end
   [breaks,coefs,l,k,d] = ppbrk(f);

   if l==1, if form(1)=='r', d=d-1; end
            jumps = zeros(d*mx,nx); return
   end
   index = sorted(breaks(2:l+1),xs); cands = find(index>0&index<l);
   index1 = find(xs(cands)-breaks(index(cands)+1)==0);
   if ~isempty(index1)
      index = cands(index1);
      if form(1)=='p'
         jumps(:,index) = ppual(f,xs(index)) - ppual(f,xs(index),'left');
      else
         d = d-1;
         jumps(:,index) = ratjump(ppual(f,xs(index)),ppual(f,xs(index),'left'));
      end
   end
otherwise
   error('The form of F is not (yet) recognized.')
end

if tosort>0,  jumps(:,ix) = jumps; end
jumps = reshape(jumps,d*mx,nx);

function jumps = ratjump(vplus,vminus)
%RATJUMP Compute jump in rational function s/w from the values of [s;w] .

d = size(vplus,1)-1;
jumps = (vplus(1:d,:).*repmat(vminus(d+1,:),d,1) ...
         -vminus(1:d,:).*repmat(vplus(d+1,:),d,1))./ ...
                          repmat(vplus(d+1,:).*vminus(d+1,:),d,1);
