% Input:
% A: Matrix A (of the continuous dynamics)
% B: Matrix B (of the continuous dynamics)
% t: double (stepsize)
% t: scalar (stepsize)
% L: Matrix (Direction vectors - rowwise)
% N: Integer (number of sets to generate)
% U: function handle (Support function of free input)
% X0: function handle (Support function of initial set)
% normvalue: 1/2/Inf (norm to use)
% options: cell array (parameters)

% Output:
% omegas: cell array of function handles (reachable sets represented by their support function)
% results: matrix (containing in each row the evaluation values for the
% correspondent direction in L)

% Description:
% Computes reachable sets Omega0 to OmegaN for stepwise constant u.
% Determines the initial overapproximation and delegates the flowpipe computation to subroutines.

function [omegas,results] = algoConstU(A,B,t,L,N,U,X0,normvalue,options)

starttime = tic;

% initial values
C = [ A B; zeros(size(B,2),size(A,2)) zeros(size(B,2),size(B,2))];
eC = expm(t*C);
Ad = eC(1:size(A,1),1:size(A,2));
Bd = eC(1:size(A,1),size(A,2)+1:size(C,2));

if options{23} == 1
    Omega0 = calcPinitConstU(Ad,X0,normvalue,t,L,options{32});
else
    V = @(l) U(B'*l);
    [~,phi2,~] = calcphis(abs(A),t);
    [~,~,edA] = calcphis(A,t); 
    [~,~,epsilonpsi,e_p,e_n] =calcepsilon(phi2,A,edA,X0,V);

    Omega0 = calcomega0delta(edA,t,V,X0,epsilonpsi,e_p, e_n,0,optimset('Display','off','Algorithm','active-set'));
end;
fprintf('initial values calculated after: %f \n\n',toc(starttime));

omegas = [];
results = [];

% reachability analysis
starttime = tic;
if strcmp(options{1} , 'Evaluation')
    
    fprintf('Using direct evaluation \n');
    results = evaluateConstU(L,Ad,Bd,U,Omega0,N);
    fprintf('reachable sets evaluated after: %f \n\n',toc(starttime));

end;

starttime = tic;
if strcmp(options{1} , 'Construction')
    
    fprintf('Using construction of the reachable sets \n');
    omegas = computeReachableSetsConstU(Ad,Bd,N,Omega0,U);
    fprintf('reachable sets computed after: %f \n\n',toc(starttime));
    set(0,'RecursionLimit',N+100) ;
end;
