function df=diff3p(f,h,o,d)
%DIFF3P calculate numerical (or partial) derivate of one (or two) variable function
%of f(x)(or f(x,y)) by the three points formula.
%   The f is the value of f(x) (or f(x,y)) at discrete data. The return value, df, is
%   the same dimenson as f, h is the spacing between points in each direction.
%   o - order, the number of the order of derivate. d - dimension, stands for derivation
%   with respect to (w.r.t.) the first argument or second. h,o and d may be in default,
%   in default of them, they are all equal 1.
%
%	 If f is a vector, df is the numerical derivate of function f(x) at the discrete point
%   x1,x2,...,xn. If f is a matrix ,when d=1,or default, the every row of df is the 
%   numerical derivate of function f(x) or partial derivate of function f(x,y) w.r.t. x, 
%   when d=2, the every column of df is the numerical partial derivate of function f(x,y)
%   w.r.t. y. The row of the matrix is the value of f(x,y) at x, the column is the value 
%   of f(x,y) at y.
%  
%   Author Dongyi Liu. Data 2000/9/20.

if nargin<2; h=1; end
if nargin<3; o=1; end
if nargin<4; d=1; end
dim=size(f);
if(dim(2)==1&d~=1)
   strf=setstr(inputname(1));
   strd=setstr(inputname(4));
   fprintf('the dimension of %s doesn''tmatch with %s\n',strf,strd);
   df=[];
   return;
end
if dim(2)==1;
   f=f.';
end
if all(dim>[1,1]);
   if d==1;
      dd=dim(2);
   else
      dd=dim(1);
   end
elseif dim(1)==1; 
   dd=dim(2);
else
   dd=dim(2);
end
if d==2
   f=f.';
end
df(:,1)=-3.*f(:,1)+4.*f(:,2)-f(:,3);
df(:,dd)=f(:,dd-2)-4.*f(:,dd-1)+3.*f(:,dd);  
for k=2:dd-1
   df(:,k)=f(:,k+1)-f(:,k-1);
end
df=df./(2*h);
if dim(2)==1; f=f.';end
if d==2; df=df.';f=f.';end
if o~=1
   for kk=2:o
      df=diff3p(df,h,o-1,d);
   end
end