%
%   version 04.03.02
%
% calc_A_D_B calculates the matrices A=df/dz, D=d(x,t)/dx and B=df/dx
% Parameter:
%    daefile - subroutine describing the DAE (see daefile.m)
%       the dae is given by f((d(x,t)',x,t)
%    y - first argument of f
%    x - second argument of f
%    t - third argument of f
%    varargin - additional arguments of daefile
%
function [dfy,D,B]= calc_A_D_B(daefile,y,x,t,s_A,s_D,s_B,varargin)
% if ~(s_A & s_B)
%     nz=length(y);
%     nx=length(x);
%     zx=[y;x];
%     jthresh=0;
%     f0=feval('ff',t,zx,daefile,nz,nx,varargin{:});
%     tresh=ones(size(zx,1),1);
%     %disp('vor numjac A+B');
%     [dfdy,fac,g1,nF] = numjac('ff',t,zx,f0,tresh,[],0,[],[],daefile,nz,nx,varargin{:});
%     dfy=dfdy(1:nx,1:nz);
%     B=dfdy(1:nx,nz+1:nz+nx);
% end;
%
%   calculation of f(z,x,t) if numerical computation of A or/and B is required
%
if ~(s_A & s_B)
    [f0]=feval(daefile,y,x,t,'',varargin{:});
end;
if s_A
    dfy=feval(daefile,y,x,t,'dfy',varargin{:});
else
    tresh=ones(size(y,1),1);
    %tresh=(abs(y)+1)*sqrt(eps);
    [dfy,fac,g1,nF] = numjac('fnachz',t,y,f0,tresh,[],0,[],[],x,daefile,varargin{:}); 
end;
if s_B
    B=feval(daefile,y,x,t,'dfx',varargin{:});
else
    tresh=ones(size(x,1),1);
    %tresh=(abs(x)+1)*sqrt(eps);
    [B,fac,g1,nF] = numjac('fnachx',t,x,f0,tresh,[],0,[],[],y,daefile,varargin{:});
end;
if s_D
    D=feval(daefile,y,x,t,'dx',varargin{:});
else
    f0=feval(daefile,y,x,t,'d',varargin{:});
    tresh=ones(size(x,1),1);
    %tresh=(abs(x)+1)*sqrt(eps);
    [D,fac,g1,nF] = numjac('fd',t,x,f0,tresh,[],0,[],[],y,daefile,varargin{:}); 
end