function out=jacob_inputs(m,varargin)
%JACOB_INPUTS computes the jacobian of the model
%    
%    The result OUT is an array of size number of outputs *
%    number of inputs containing the derivatives of the
%    outputs with respect to the inputs.
%    
%    JACOB_INPUTS(M,X) computes for each output the jacobian
%    of the model with respect to the x (the input of the
%    linear) at point X.
%    
%    JACOB_INPUTS(M,X,J) computes for output J the jacobian of
%    the model with respect to the x (the input of the linear)
%    at point X.
%    
%    

%MAN_PAGE_BEGIN
%
%   @purpose    computes the jacobian of the model.
%   The result <code>out</code> is an array of size number of outputs * number of inputs
%	containing the derivatives of the outputs with respect to the inputs.
%
%   @synopsis   jacob_inputs(m,x)
%   @description computes for each output the jacobian of
%  the model with respect to the x (the input of the
%  linear) at point <CODE>x</CODE>.
%
%   @synopsis   jacob_inputs(m,x,j)
%   @description computes for output J the jacobian of
%  the model with respect to the x (the input of the
%  linear) at point <CODE>x</CODE>.
%
%MAN_PAGE_END

%   Copyright (c) 1998 by Antoine Duchateau			
%	$Revision: 0.1$  $Date: 10/12/1997
%
%	All rights reserved

error(checkargs([{m} varargin],[2 3],'series','numeric','numeric'));

check(m);

x = varargin{1};
if (nargin == 2)|isempty(varargin{2}),
	j = 1:get(m,'n_out');
else
	j = varargin{2};
end

lsy1 = length(m.links.Sy1);
lsy2 = length(m.links.Sy2);
lsy3 = length(m.links.Sy3);

out = zeros(get(m,'n_out'),get(m,'n_in'));

out(1:lsy1,1:lsy1) = ones(lsy1);

e1 = eval(m.mapping1,x(:,m.links.Sm1));
j1 = jacob_inputs(m.mapping1,x(:,m.links.Sm1));

out(lsy1+(1:lsy2),m.links.Sm1) = out(lsy1+(1:lsy2),m.links.Sm1) + j1(m.links.Sy2,:);

um2 = [x(:,m.links.Sm2) e1(:,m.links.Sm21)]
if isfield(m.links,'perm1'),
	um2 = um2(:,m.links.perm1);
end

[nouse, invperm1] = sort(m.links.perm1);

j2 = jacob_inputs(m.mapping2,um2);
j2 = j2(:,invperm1);

out(lsy1+lsy2+(1:lsy3),m.links.Sm2) = out(lsy1+lsy2+(1:lsy3),m.links.Sm2) + j2(:,length(m.links.Sm21)+1:end);

out(lsy1+lsy2+(1:lsy3),m.links.Sm1) = out(lsy1+lsy2+(1:lsy3),m.links.Sm1) + j2(:,1:length(m.links.Sm21)) * j1;

if isfield(m.links,'perm2'),
	out = out(m.links.perm2,:);
end

out = out(j,:);
