function [sys,J] = jacobian2ss(J,model,options,Ts,iostruct)
%JACOBIAN2SS
%
% SYS = JACOBIAN2SS Obtains a linear model from a the Jacobian structure
% information from linearization and the IO pairings.

%  Author(s): John Glass
%  Copyright 1986-2004 The MathWorks, Inc.
% $Revision: 1.1.6.10 $ $Date: 2004/12/10 19:40:55 $

%% Extract the upper parts of the LFT
a = J.A; b = J.B; c = J.C; d = J.D;

%% Get the model sample rates and block dimensions
[ny,nu] = size(d);
nxz = size(a,1);
Tsx = J.Ts(1:nxz);
Tsy = J.Ts(nxz+1:end);

%% Extract lower parts of the lft
E = J.Mi.E;
F = J.Mi.F;
G = J.Mi.G;
H = J.Mi.H;

%% Get the state names
StateName = J.StateName;

%% Compute the input and output index vectors
InputIdx = [J.Mi.InputIdx+1;nu+1];
OutputIdx = [J.Mi.OutputIdx+1;ny+1];
BlockHandles = J.Mi.BlockHandles;

%% Construct the input to block list
InputList = [];
for ct = 1:(length(InputIdx)-1)
    if (InputIdx(ct) ~= InputIdx(ct+1))
        InputList = [InputList;BlockHandles(ct)*ones(InputIdx(ct+1)-InputIdx(ct),1)];
    end
end

%% Construct the output to block list
OutputList = [];
for ct = 1:(length(OutputIdx)-1)
    if (OutputIdx(ct) ~= OutputIdx(ct+1))
        OutputList = [OutputList;BlockHandles(ct)*ones(OutputIdx(ct+1)-OutputIdx(ct),1)];
    end
end

%% Perform block reduction if requested
if strcmpi(options.BlockReduction,'on')
    %% Removing blocks that are not part of the linearization
    BlocksInPath = minjacobian(J);
    %% Tack the BlocksInPath to the Jacobian data
    J.Mi.BlocksInPath = BlocksInPath;

    %% Get the blocks that are not in the linearization
    blocks = find(~BlocksInPath);
    nblocks = length(blocks);

    %% Extract index elements, tack on the number of states, inputs, and
    %% outputs to account for the last block in the list.  Convert from
    %% C-indexing to Matlab-indexing.
    InputIdx = [J.Mi.InputIdx+1;nu+1];
    OutputIdx = [J.Mi.OutputIdx+1;ny+1];

    indx = [];indu = [];indy = [];
    for ct = 1:nblocks

        %% Get the indices in to the states that should be eliminated for
        %% each block being removed.
        try
            state_ind = find(strcmp(getfullname(J.Mi.BlockHandles(blocks(ct))),StateName));
            if ~isempty(state_ind)
                indx = [indx;state_ind];
            end
        end
        %% Get the indices for the inputs to blocks that are being removed.
        if (InputIdx(blocks(ct)) ~= InputIdx(blocks(ct)+1))
            indu = [indu,InputIdx(blocks(ct)):(InputIdx(blocks(ct)+1)-1)];
        end
        %% Get the indices for the outputs to blocks that are being
        %% removed.
        if (OutputIdx(blocks(ct)) ~= OutputIdx(blocks(ct)+1))
            indy = [indy,OutputIdx(blocks(ct)):(OutputIdx(blocks(ct)+1)-1)];
        end
    end

    %% Eliminate states
    if ~isempty(indx)
        a(indx,:) = [];
        a(:,indx) = [];
        StateName(indx) = [];
        b(indx,:) = [];
        c(:,indx) = [];
        Tsx(indx) = [];
        nxz = length(a);
    end
    %% Eliminate block inputs
    if ~isempty(indu)
        b(:,indu) = [];
        d(:,indu) = [];
        E(indu,:) = [];
        F(indu,:) = [];
        InputList(indu,:) = [];
        nu = size(b,2);
    end
    %% Eliminate block outputs
    if ~isempty(indy)
        c(indy,:) = [];
        d(indy,:) = [];
        E(:,indy) = [];
        G(:,indy) = [];
        OutputList(indy,:) = [];
        Tsy(indy) = [];
        ny = size(c,1);
    end
else
    %% Tack the BlocksInPath to the Jacobian data
    J.Mi.BlocksInPath = ones(size(J.Mi.BlockHandles));
end

%% Default is least common multiple of all sample times.  Compute the unique
%% sample times.  If there are no states also assume that the sample time
%% is zero so the d2d conversions are not performed.
Ts_all = [Tsx;Tsy];
Tuq = unique(Ts_all(Ts_all >= 0));
if (Ts == -1)
    if (length(Tsx) == 0) || (length(Tuq) == 0) || (max(Tuq) == 0)
        Ts = 0;
    else
        Ts = local_vlcm(Tuq);
        if isempty(Ts)
            %% If model has fixed step size specified use that
            modelSolver = get_param(model,'Solver');
            modelFixedStepSize = get_param(model,'FixedStep');
            if ( ~strcmpi(modelFixedStepSize,'auto') && ...
                    ( strcmp(modelSolver,'ode5') || ...
                    strcmp(modelSolver,'ode4') || ...
                    strcmp(modelSolver,'ode3') || ...
                    strcmp(modelSolver,'ode2') || ...
                    strcmp(modelSolver,'ode1') || ...
                    strcmp(modelSolver,'FixedStepDiscrete')))
                warning(['No sample time found in the model. ',...
                    'Defaulting to fixed step size']);
                Ts = str2num(modelFixedStepSize);
            else
                warning('No sample time found in the model.  Defaulting to 1');
                Ts = 1;
            end
        end
    end
end

%% Linearize without discrete states if the user desires, otherwise perform
%% the discrete linearization if there are discrete states.
if ((Ts == 0)) && (strcmp(options.IgnoreDiscreteStates,'on') || ...
        (all(Tsx==0)))

    % Remove discrete states
    if any(Tsx)
        warning(['Ignoring discrete states']);
        ix = find(Tsx ~= 0);
        a(ix,:) = [];
        a(:,ix) = [];
        b(ix,:) = [];
        c(:,ix) = [];
        StateName(ix) = [];
    end

    P = speye(size(d,1)) - d*E;
    Q = P \ c;
    R = P \ d;

    % Close the LFT
    a = a + b * E * Q;
    b = b * (F + E * R * F);
    c = G * Q;
    d = H + G * R * F;
else
    [ny,nu] = size(d);
    nxz = size(a,1);
    Tslist = [ Tuq ; Ts ];
    Eslow  = E;
    %% Get the sampling methodology
    RateConvMethod = options.RateConversionMethod;

    for k=1:length(Tslist)-1
        % Start with the fastest rate
        ts_current = Tslist(k);
        ts_next    = min(Ts, Tslist(k+1));

        xix = find(Tsx == ts_current);
        nix = find(Tsx ~= ts_current);

        % Find the indecies of the blocks at the current sample rates these
        % indecies correspond the the outputs of the blocks
        ind_Tsy_current = find(Tsy==ts_current);

        % Find the indecies of the inputs of the blocks
        blks = unique(OutputList(ind_Tsy_current));
        ind_Tsu_current = [];
        for ct = 1:length(blks);
            ind_Tsu_current = [ind_Tsu_current;find(blks(ct) == InputList)];
        end

        % Close the fastest loops
        [ix,jx,px] = find(Eslow);
        ux1   = ismember(ix,(ind_Tsu_current));
        ux2   = ismember(jx, ind_Tsy_current);
        ux = ux1 & ux2;
        Efast = sparse(ix(ux),jx(ux),px(ux),nu,ny);

        % And mark them as closed (remove them from interconnection matrix)
        Eslow = Eslow - Efast;
        P = speye(ny) - d*Efast;

        % Update the block outputs to indicate the rate conversion
        Tsy(ind_Tsy_current) = ts_next;

        % But leave the matrices "full-sized" for later connection
        c = P \ c;
        d = P \ d;
        a = a + b * Efast * c;
        b = b * (speye(nu) + Efast*d);

        % if the target rate is the slowest rate we are done
        if ts_current ~= ts_next

            % Otherwise use c2d/d2d/d2c to reach the next rate
            nxfast = length(xix);
            atmp = full(a(xix,xix));

            if isinf(ts_current)
                if ts_next ~= 0
                    Phi = eye(size(atmp));
                    Gam = zeros(size(b(xix,:)));
                else
                    Phi = zeros(size(atmp));
                    Gam = zeros(size(b(xix,:)));
                end
            else
                %% Find the input and output channels that
                %% correspond the states that are being converted
                indu = find(sum(b(xix,:)));
                indy = find(sum(c(:,xix)'));
                btmp = full(b(xix,indu));
                ctmp = full(c(indy,xix));
                dtmp = full(d(indy,indu));
                sysold = ss(atmp,btmp,ctmp,dtmp,ts_current);

                %% Switch over rate conversion methods
                switch RateConvMethod
                    case 'zoh'
                        if ts_current ~= 0
                            lastwrn = lastwarn; warning('off');
                            if ts_next ~= 0
                                sysnew = d2d(sysold,ts_next);
                            else
                                sysnew = d2c(sysold);
                            end
                            warning('on'); lastwarn(lastwrn);

                            if length(sysold.a) ~= length(sysnew.a)
                                %% Expand the number of states to account for states that have been added
                                deltanstates = length(sysnew.a)-length(sysold.a);
                                nx_iter = size(a,1);
                                laststate = xix(end);
                                %% Insert rows into the a matrix
                                a = [a(1:laststate,:);zeros(deltanstates,nx_iter);a(laststate+1:nx_iter,:)];
                                %% Insert columns into the a matrix
                                a = [a(:,1:laststate),zeros(nx_iter+deltanstates,deltanstates),a(:,laststate+1:nx_iter)];
                                %% Insert rows into the b matrix
                                b = [b(1:laststate,:);zeros(deltanstates,nu);b(laststate+1:nx_iter,:)];
                                %% Insert columns into the c matrix
                                c = [c(:,1:laststate),zeros(ny,deltanstates),c(:,laststate+1:nx_iter)];
                                %% Insert rows into the StateName field
                                EmptyStates = cell(deltanstates,1);
                                for ct = 1:deltanstates
                                    EmptyStates{ct} = '?';
                                end
                                if laststate == nx_iter;
                                    StateName = {StateName{:},EmptyStates{:}};
                                else
                                    StateName = {StateName{1:laststate},EmptyStates{:},StateName{laststate+1:nx_iter}};
                                end
                                %% Insert rows into the Tsx vector
                                Tsx = [Tsx(1:laststate);zeros(deltanstates,1);Tsx(laststate+1:nx_iter)];
                                
                                %% Update the xix and nix vectors
                                ind = find(nix > laststate);
                                xix = [xix;(1:deltanstates)'+laststate];
                                nix(ind) = nix(ind)+1;
                            end
                        else
                            sysnew = c2d(sysold,ts_next);
                        end
                    case {'tustin','prewarp'}
                        %% Find the input and output channels that
                        %% correspond the states that are being converted
                        indu = find(sum(b(xix,:)));
                        indy = find(sum(c(:,xix)'));
                        btmp = full(b(xix,indu));
                        ctmp = full(c(indy,xix));
                        dtmp = full(d(indy,indu));
                        sysold = ss(atmp,btmp,ctmp,dtmp,ts_current);
                        if ts_current ~= 0
                            if ts_next ~= 0
                                switch RateConvMethod
                                    case 'tustin'
                                        sysmid = d2c(sysold,'tustin');
                                        sysnew = c2d(sysmid,ts_next,'tustin');
                                    case 'prewarp'
                                        PreWarpFreq = LocalComputeScalarVal(options.PreWarpFreq);
                                        sysmid = d2c(sysold,'prewarp',PreWarpFreq);
                                        sysnew = c2d(sysmid,ts_next,'prewarp',PreWarpFreq);
                                end
                            else
                                switch RateConvMethod
                                    case 'tustin'
                                        sysnew = d2c(sysold,'tustin');
                                    case 'prewarp'
                                        PreWarpFreq = LocalComputeScalarVal(options.PreWarpFreq);
                                        sysnew = d2c(sysold,'prewarp',PreWarpFreq);
                                end
                            end
                        else
                            switch RateConvMethod
                                case 'tustin'
                                    sysnew = c2d(sysold,ts_next,'tustin');
                                case 'prewarp'
                                    PreWarpFreq = LocalComputeScalarVal(options.PreWarpFreq);
                                    sysnew = c2d(sysold,ts_next,'prewarp',PreWarpFreq);
                            end
                        end
                        %% Remove any statenames that may have been used in the
                        %% conversion for Tustin routines
                        for ct = 1:length(xix)
                            StateName{xix(ct)} = '?';
                        end
                end
                a(xix,xix) = sysnew.a;
                b(xix,indu) = sysnew.b;
                c(indy,xix) = sysnew.c;
                d(indy,indu) = sysnew.d;
            end
        end
        Tsx(xix) = ts_next;
    end

    % Close the remaining loops if they are any still open
    if any(any(Eslow))
        % Close the final loop if needed
        P = speye(ny) - d*Eslow;
        c = P \ c;
        d = P \ d;
        a = a + b * Eslow * c;
        b = b * (speye(nu) + Eslow*d);
    end

    %% Apply analysis input and output selection
    b = b * F;
    c = G * c;
    d = H + G * d * F;
end

%% Reorganize the IO due to back propagation
Bnew = b(:,iostruct.InputInd);
Cnew = c(iostruct.OutputInd,:);
Dnew = d(iostruct.OutputInd,:);
Dnew = Dnew(:,iostruct.InputInd);

%% Create the linear model
sys = ss(full(a),full(Bnew),full(Cnew),full(Dnew),Ts);
%% Set the labels and loop over to find unique names
sys.InputName = iostruct.InputName;
sys.OutputName = iostruct.OutputName;
sys.StateName = StateName;

%---

function M = local_vlcm(x)
%% VLCM  find least common multiple of several sample times

%% Protect against a few edge cases, remove zeros before computing LCM
x(~x) = [];
x(isinf(x)) = [];
if isempty(x), M = []; return; end;

[a,b]=rat(x);
v = b(1);
for k = 2:length(b), v=lcm(v,b(k)); end
d = v;

y = round(d*x);         % integers
v = y(1);
for k = 2:length(y), v=lcm(v,y(k)); end
M = v/d;

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LocalComputeScalarVal
%% Computes a scalar variable in the base workspace
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function valout = LocalComputeScalarVal(val)

if ischar(val)
    try
        valout = evalin('base',val);
    catch
        error('slcontrol:InvalidSampleTimeExpression',...
            sprintf('Invalid expression %s for the sample time option.',...
            val))
    end
else
    valout = val;
end