function [U] = GridDef_add(lbmub,varargin)

% This functions augments the already chosen sampling domain for the
% following sampling methods:
%  - Latin Hypercube sampling with the function lhsdesign.m
%  - Orthogonal sampling with the function orthogonal_sample.m
%  - Sobol sampling with the function sobolset.m
%
% The input of this function is given by
%       lbmub:	3 x ntot matrix, with first row determining the minimum
%               value, the second row the reference value and the third row
%               the maximum value. This matrix includes all independent
%               variables, also the unvaried.
%
%       indin:  Indices ofthe lbmub matrix, which should be varied
%
%       opt:	Structure containing the options of the system, includes
%               .type:  String containing CG, LHS or OD
%               .Np:    LHS:    Total number of points
%                       OD:     Number of points per subdivision
%                       RND:    Random numbers generated
%                       SOB:    Sobol sequence sampling
%               .Ns     OD:     Number of subdivision per variable
%
%       prev:   Previous sampling points for augmenting LHS

% Initialize the global variable SOB giving the sobol sequence defined in
% the first section
global SOB

% Safety function to use all available input points and a LHS with 1000 points in
% the case no inputs or grids are defined.
if nargin < 2
    indin = 1:size(lbmub,2);
    opt.type = 'LHS';
    opt.Np = 1e3;
    opt.val = 'y';
elseif nargin < 3
    indin = varargin{1};
    opt.type = 'LHS';
    opt.Np = 1e3;
    opt.val = 'y';
elseif nargin < 4
    indin = varargin{1};
    opt = varargin{2};
else
    indin = varargin{1};
    opt = varargin{2};
    prev = varargin{3};
end

% Definition of the number of varied independent variables, the total
% number of variables, and the number of previous samples
nvar = length(indin);
nsamp = opt.Np;
nprev = size(prev,1);

switch opt.type                     
    case 'RND'
        % Performing the random sampling
        U.grid = rand(nsamp,nvar);
        U.var  = U.grid.*(ones(nsamp,1)*diff(lbmub([1 3],indin)))+ones(nsamp,1)*lbmub(1,indin);
        
        % Definition of the vectors for the feed composition without any change
        U.out = ones(nsamp,1)*lbmub(2,:);

        % Definition of the adjusted matrix
        U.out(:,indin) = U.var;
        
        % Rearragnement of the variables which need to be changed
        if isfield(opt,'dep')
            for k = 1:size(opt.dep,1)
                U.out(:,opt.dep(k,2)) = U.out(:,opt.dep(k,1))*opt.dep(k,3).*...
                                       (1+(U.grid(:,opt.dep(k,2))-.5)*2*opt.dep(k,4));%+opt.dep(k,4);
                
            end
        end
                
                
    case 'LHS'
        % Performing the LHS
        U.grid = lhsaugment(prev,nsamp);
        U.grid = U.grid(end-nsamp+1:end,:);
        U.var  = U.grid.*(ones(nsamp,1)*diff(lbmub([1 3],indin)))+ones(nsamp,1)*lbmub(1,indin);
        
        % Definition of the vectors for the feed composition without any change
        U.out = ones(nsamp,1)*lbmub(2,:);

        % Definition of the adjusted matrix
        U.out(:,indin) = U.var;
        
        % Rearragnement of the variables which need to be changed
        if isfield(opt,'dep') 
            for k = 1:size(opt.dep,1)
                if any(indin == opt.dep(k,2))
                    U.out(:,opt.dep(k,2)) = U.out(:,opt.dep(k,1))*opt.dep(k,3).*...
                                           (1+(U.grid(:,opt.dep(k,2))-.5)*2*opt.dep(k,4));%+opt.dep(k,4);
                end
            end
        end
        
    case 'SOB'
        % Calculating the new points of the sobols equence
        U.grid = SOB(nprev+1:nprev+nsamp,:);
        U.var  = U.grid.*(ones(nsamp,1)*diff(lbmub([1 3],indin)))+ones(nsamp,1)*lbmub(1,indin);
        
        % Definition of the vectors for the feed composition without any change
        U.out = ones(nsamp,1)*lbmub(2,:);

        % Definition of the adjusted matrix
        U.out(:,indin) = U.var;
        
        % Rearragnement of the variables which need to be changed
        if isfield(opt,'dep') 
            for k = 1:size(opt.dep,1)
                if any(indin == opt.dep(k,2))
                    U.out(:,opt.dep(k,2)) = U.out(:,opt.dep(k,1))*opt.dep(k,3).*...
                                           (1+(U.grid(:,opt.dep(k,2))-.5)*2*opt.dep(k,4));%+opt.dep(k,4);
                end
            end
        end
end
