function v=ppual(pp,xx,left) % % v = ppual(pp,xx[,left]) % % returns the value at x of the pp function described by pp . The % function is made left-continuous if the optional third argument is present. % copyright 1987-1993 Carl de Boor % C. de Boor / latest change: June 19, 1989 % C. de Boor / latest change: December 1, 1990 (correct misprint in last line) % C. de Boor / latest change: 19 se 94 (temporarily programmed around the % breaks ([1 1 1 1 ... ]) discontinuity ) % C. de Boor / latest change: 19 oc 94 (optional argument for left continuity) % if necessary, sort xx xs=xx(:)';ix=length(xs);tosort=0; if (length(find(diff(xs)<0))>0), tosort=1;[xs,indx]=sort(xs); end % take apart pp [x,c,l,k,d]=ppbrk(pp); % for each data point, compute its breakpoint interval if nargin<3 index=max([sorted(x(1:l),xs);ones(1,ix)]); else index=max([l+1-sorted(-x(l+1:-1:2),-xs);ones(1,ix)]); end % now go to local coordinates % xs=xs-x(index); if all(index==1), xs = xs - x(1)*index; else xs = xs - x(index); end % ... and apply nested multiplication: if (d>1), ones(d,1)*xs; xs=ans(:)'; 1+d*ones(d,1)*index+[-d:-1]'*ones(1,ix); index=ans(:); end if all(index==1), vv=c(1,1)*index; for i=2:k, vv = xs.*vv + c(1,i)*index; end else vv=c(index,1)'; for i=2:k, vv = xs.*vv + c(index,i)'; end end v=zeros(d,length(xx));v(:)=vv; if tosort>0,[junk,indx]=sort(indx); v=v(:,indx);end