%
%   02.12.02
%
%   nl_consi delivers the nonlinear function and its Jacobian
%   to calculate consistent initial values for DAEs
%
%   The nonlinear function is given by
%
%              f(y,x,t0) = 0
%   W1(BD^-(y-D'Ux)+f_t) = 0
%               (I - R)y = 0
%                   z(x) = 0
%          DP_1(x-alpha) = 0
%
%   yx - vector of [y;x]
%   newJac - 0 return of right only
%            1 return of right and Jac
%   n - dimension of f
%   daefile - description of the DAE
%   t - timepoint
%   z - function with additional constraints
%   alpha - "initial" value
%   s_A -   1  analytic representation of df/dy=A,
%           0  numerical calculation of A
%   s_D -   1  analytic representation of d d(x,t)/dx=D in daefile
%           0  numerical calculation of D
%   s_B  -  1  analytic representation of df/dx=B in daefile
%           0  numerical calculation of B
%   s_Bt -  0  numerical calculation of B_t
%           -1 not time dependent
%   s_Dt -  0  numerical calculation of D_t
%           -1 not time dependent
%   s_ft -  1  analytic representation of d f(z,x,t)/dt=f_t in daefile
%           0  numerical calculation of f_t
%           -1 not time dependent
%   given_index - 
%   threshold -
%   fid -
%   varargin -
%
function [ok_f,right,JacQR]=nl_consi(yx,newJac,print_level,n,daefile,t,z,alpha,...
    s_A,s_D,s_B,s_Bt,s_Dt,s_ft,given_index,threshold,varargin);
%
persistent W1 A B Dm U DP1 ImR ips index Dprime rankDP1 rankD rankW1 BDm
%
%   output file
%
fid=file_open_date(print_level,'nl_consi');
fprintf(fid,'%s \n',headline(daefile));
%

ok_f=1;
%
fprintf(fid,'%s \n','in nl_consi');
%
%   extract y and x from yx
%
dimyx=size(yx,1);
m=dimyx-n;
y=yx(1:m);
x=yx(m+1:end);
%
%----------------- handling of additional conditions
zx=feval(z,x);
nz=size(zx,1);
%-----------------
if newJac
    %
    [ok_s,index,A,R,D,rankD,B,Dm,Q0,P0,W0,DP1,rankDP1,W1,rankW1,U]=...
        matrix_sequence(daefile,'xintpol',t,s_A,s_D,s_B,given_index,threshold,fid,2,t,x,varargin{:});
    %disp(['nl_consi: index = ',num2str(index),' x(2) = ',num2str(x(2),'%25.17e')]);
    if ~ok_s
        ok_f=0;
        %
        %   only to assign something
        %
        right=[];
        JacQR=[];
        disp('nl_consi: Index determination not possible!');
        fclose(fid);
        return;
    end;
    % abtesten rankD und size(A,2)
    if rankD==m
        ImR_and_zero=[];
    else
        ImR=speye(m)-R;
        ImR_and_zero=[ImR,sparse(m,n)];
    end;
    switch index
    case 2
        print_mat(fid,' U =',U);
        print_mat(fid,' Dm =',Dm);
        print_mat(fid,' DP1 =',DP1);
        print_mat(fid,' R =',R);
        switch s_Dt
        case -1
            Dprime=sparse(m,n);
        case 0
            Dprime=cal_Dt(s_D,daefile,y,x,t,varargin{:});
            print_mat(fid,' D''t =',Dprime);
        otherwise
            error('wrong value for s_Dt');
        end;
        switch s_Bt
        case -1
            Bprime=sparse(n,n);
        case 0
            Bprime=cal_Bt(s_B,daefile,y,x,t,varargin{:});
            print_mat(fid,' B''t =',Bprime);
        otherwise
            error('wrong value for s_Bt');
        end;
        BDm=B*Dm;
        Jac=[A,B; W1*BDm, W1*(BDm*Dprime*U+Bprime); ImR_and_zero];
        if nz > 0
            zpr=z_x(z,x,zx);
            Jac=[Jac;sparse(nz,m),zpr];
            rankJac=n+rankW1+(size(A,2)-rankD)+nz;
            if nz==rankDP1
                BB=[];
            else
                BB=[sparse(m,m),DP1];
            end;
        else
            Jac=[Jac;sparse(m,m),DP1];
            rankJac=n+rankW1+(m-rankD)+rankDP1;
            BB=[];
        end;
    case 1
        print_mat(fid,' R =',R);
        Jac=[A,B; ImR_and_zero];
        if nz > 0
            zpr=z_x(z,x,zx);
            Jac=[Jac; sparse(nz,m),zpr];
            rankJac=n+(m-rankD)+nz;
            if nz==rankDP1
                BB=[];
            else
                BB=[sparse(m,m),DP1];
            end;
        else
            Jac=[Jac;sparse(m,m),DP1];
            rankJac=n+(m-rankD)+rankD;
            BB=[];
        end;
    case 0
        Jac=[A,B ;sparse(m,m),speye(n)];
        rankJac=n+n;
        BB=[];
    otherwise
        ok_f=0;
        disp(['nl_consi: The index is probably ',num2str(index),' > 2 !']);
        %
        %   only to assign something
        %
        right=[];
        JacQR=[];
        fclose(fid);
        return
    end;
    print_mat(fid,' Jac =',Jac);
    [JacQ,JacR,ips,ok]=reduction(Jac',rankJac,BB',threshold);
    if ~ok 
        disp('nl_consi: Jacobian hasn''t full rank');
        ok_f=0;
        %
        %   only to assign something
        %
        right=[];
        JacQR=[];
        fclose(fid);
        return
    end;
    JacQR=[JacQ,JacR];
end;
%
fx=feval(daefile,y,x,t,'',varargin{:});
%
%   calculation of f(x) for index = 2,1,0
%
switch index
case 2
    %
    %   calclation of right side vector
    %
    right=fx;
    text=' is computed ';
    num='numerically';
    ana='analytically';
    switch s_ft
    case -1
        fprime=sparse(n,1);
    case 0
        fprime = cal_ft(daefile,y,x,t,varargin{:});
        text_ft=[' f_t',text,ana];
    case 1
        fprime=feval(daefile,y,x,t,'dft',varargin{:});
        text_ft=[' f_t',text,num];
    end;
    right=[right;W1*(BDm*(y-Dprime*U*x)+fprime)];
    if rankD~=m
        right=[right;ImR*y];
    end;
    if nz > 0
        if nz==rankDP1
            right=[right;zx];
        else
            right=[right;zx;DP1*(x-alpha)];
        end;
    else
        right=[right;DP1*(x-alpha)];
    end;
    %
    %   deleting of the unused components
    %
    right(ips)=[];
case 1 
    %
    %   calclation of right side vector
    %
    right=fx;
    if rankD~=m
        right=[right;ImR*y];
    end;
    if nz > 0
        if nz==rankDP1
            right=[right;zx];
        else
            right=[right;zx;DP1*(x-alpha)];
        end;
    else
        right=[right;DP1*(x-alpha)];
    end;
    %
    %   deleting of the unused components
    %
    right(ips)=[];
case 0
    %
    %   calclation of right side vector
    %
    right=fx;
    right=[right;x-alpha];
end;
print_mat(fid,'right = ',right');
fprintf(fid,'%s \n','out nl_consi');
if print_level==2
    fclose(fid);
end;