function prt=dfdp(x,f,p,dp,func)
% numerical partial derivatives (Jacobian) df/dp for use with leasqr
% --------INPUT VARIABLES---------
% x=vec or matrix of indep var(used as arg to func) x=[x0 x1 ....]
% f=func(x,p) vector initialsed by user before each call to dfdp
% p= vec of current parameter values
% dp= fractional increment of p for numerical derivatives
%      dp(j)>0 central differences calculated
%      dp(j)<0 one sided differences calculated
%      dp(j)=0 sets corresponding partials to zero; i.e. holds p(j) fixed
% func=string naming the function (.m) file
%      e.g. to calc Jacobian for funcion expsum
%   prt=dfdp(x,f,p,dp,'expsum')
%----------OUTPUT VARIABLES-------
% prt= Jacobian Matrix prt(i,j)=df(i)/dp(j)
%================================
m = length(x); n = length(p);                       %dimensions
ps = p; prt = zeros(m,n); del = zeros(n,1);         % initialise Jacobian to Zero
for j = 1:n,
  del(j) = dp(j) .* p(j);                    %cal delx=fract(dp)*param value(p)
  if p(j) == 0,
    del(j) = dp(j);                          %if param=0 delx=fraction
  end;
  p(j) = ps(j) + del(j);
  if del(j) ~= 0,
    f1 = feval( func, x, p );
    if dp(j) < 0,
      prt(:,j) = ( f1 - f ) ./ del(j);
    else
      p(j) = ps(j) - del(j);
      prt(:,j) = ( f1 - feval(func,x,p) ) ./ ( 2 .* del(j) );
    end;
  end;
  p(j) = ps(j);                              %restore p(j)
end;
return;
