function rangepf=orange(pf, cf, n) 
% orange(pf, cf, n)  generates a probability function for the range of n 
%                    samples of the given probability and cumulative 
%                    functions 
% 
%   If you'd like nice results, the functions had better be matched 
% 
%   It also won't work too well if the x samples on a continuous function 
%      are not evenly spaced (all generated functions are) 
% 
%   One last thing: this function can be REALLY slow. 
if ispf(pf)==0 
   error('invalid probability function') 
end 
if iscf(cf)==0 
   error('invalid cumulative function') 
end 
slope=1; 
p=pvals(pf); 
P=pvals(cf); 
x=xvals(pf); 
R=max(x)-min(x); 
newx=[0:R/2000:R]; 
newp=zeros(length(newx), 1); 
cornerp=zeros(1,4); 
rangeindex=1; 
while rangeindex<=length(newx), 
   TP=0; 
   maxindex=round(1+newx(rangeindex)*slope); 
   if maxindex>length(x) 
      maxindex=length(x); 
   end 
   if x(maxindex)-x(1)>newx(rangeindex) 
      while x(maxindex)-x(1)>newx(rangeindex) & maxindex>1 
         maxindex=maxindex-1; 
      end 
      if maxindex>1 & newx(rangeindex)>0 
         slope=(maxindex-1)/newx(rangeindex); 
      end 
   else 
      while x(maxindex)-x(1)<newx(rangeindex) & maxindex<length(x) 
         maxindex=maxindex+1; 
      end 
      if newx(rangeindex)>0 
         slope=(maxindex-1)/newx(rangeindex); 
      end 
   end 
   if maxindex>1 
      Pr=P(maxindex:length(P)); 
      pr=p(maxindex:length(p)); 
      Pnr=P(1:length(P)-maxindex+1); 
      pnr=p(1:length(p)-maxindex+1); 
      ip=perm(n, 2)*((Pr-Pnr).^(n-2)).*pr.*pnr; 
      ix=x(maxindex:length(x)); 
      NP=trapz(ix, ip); 
      if ~(NP==[]) 
         TP=TP+NP; 
      end 
   end 
   newp(rangeindex)=TP; 
   rangeindex=rangeindex+1; 
end 
rangepf=[pf(1,:); newx' newp]; 
if ispf(rangepf)==0 
   newp=newp./ttlprob(rangepf); 
   rangepf=[pf(1,:); newx' newp]; 
end