function Results = POISS(INPUT,Results_old,EstimOpt,OptimOpt)

global B_backup

% save tmp_MNL
% return

tic 

Results.bhat = [];
Results.R = [];
%Results.R_out = {};
Results.stats = [];


%% check input:

if nargin<3
    error('Too few input arguments for POISS(INPUT,EstimOpt,OptimOpt)')
end

warning off MATLAB:mir_warning_maybe_uninitialized_temporary

format shortG;
format compact;

if isfield(EstimOpt, 'NB') == 0
    EstimOpt.NB = 0;
end
EstimOpt.NVarA = size(INPUT.Xa,2);

if isfield(INPUT, 'Xnb') == 0 || numel(INPUT.Xnb) == 0
    if EstimOpt.NB > 0
        INPUT.Xnb = ones(size(INPUT.Y,1),1);
    else
        INPUT.Xnb = zeros(size(INPUT.Y,1),0);
    end
else
    INPUT.Xnb = [ones(size(INPUT.Y,1),1), INPUT.Xnb];
end
EstimOpt.NVarNB = size(INPUT.Xnb,2); % Covariates of alpha (overdispersion)

if EstimOpt.NVarNB > 0 && EstimOpt.NB == 0
   EstimOpt.NB = 1;
end

if EstimOpt.Display == 1
    disp(' ');
    disp('__________________________________________________________________________________________________________________');
    disp(' ');
    if EstimOpt.NB == 0
        disp('Estimating Poisson regression model ...')
    elseif EstimOpt.NB == 1
        disp('Estimating Negative Binomial regression model ...')
    elseif EstimOpt.NB == 2
        disp('Estimating NBP regression model ...')
    end
    if isfield(EstimOpt,'Censored') == 0 || EstimOpt.Censored < 1
        disp('without censoring.')
    else
        disp(num2str(EstimOpt.Censored, 'with censoring at %2.0f.'))
    end
end

if isfield(EstimOpt,'Censored') == 0 || EstimOpt.Censored < 1
    EstimOpt.Censored = 0;
end
if EstimOpt.NB == 2 && EstimOpt.Censored == 1
   error('Censoring not possible for NBP model') 
end
if isfield(INPUT,'W') == 0 || length(INPUT.W) ~= length(INPUT.Y)
    INPUT.W = ones(length(INPUT.Y),1);
else
    INPUT.W = INPUT.W(:);
end

if isfield(EstimOpt,'NamesA') == 0 || isempty(EstimOpt.NamesA) || length(EstimOpt.NamesA) ~= EstimOpt.NVarA
    EstimOpt.NamesA = (1:EstimOpt.NVarA)';
    EstimOpt.NamesA = cellstr(num2str(EstimOpt.NamesA));
elseif size(EstimOpt.NamesA,1) ~= EstimOpt.NVarA
    EstimOpt.NamesA = EstimOpt.NamesA';
end
if EstimOpt.NVarNB > 1
    if isfield(EstimOpt,'NamesNB') == 0 || isempty(EstimOpt.NamesNB)|| length(EstimOpt.NamesNB) ~= (EstimOpt.NVarNB-1)
        EstimOpt.NamesNB = (1:EstimOpt.NVarNB-1)';
        EstimOpt.NamesNB = cellstr(num2str(EstimOpt.NamesNB));
    elseif size(EstimOpt.NamesNB,1) ~= EstimOpt.NVarNB - 1
        EstimOpt.NamesNB = EstimOpt.NamesNB';
    end
    EstimOpt.NamesNB = [{'Cons'};EstimOpt.NamesNB];
else
    EstimOpt.NamesNB = {'Cons'};
end
if EstimOpt.NB == 2
    EstimOpt.NamesNB = [EstimOpt.NamesNB; {'P'}];
end

if isfield(INPUT,'W') == 0 || length(INPUT.W) ~= length(INPUT.Y)
    INPUT.W = ones(length(INPUT.Y),1);
% else
%     INPUT.W = INPUT.W(:);
%     INPUT.W = reshape(INPUT.W(:, ones(1, EstimOpt.NCT))', EstimOpt.NCT*EstimOpt.NP,1);
end

%% starting values 
if exist('B_backup','var') && ~isempty(B_backup) && size(B_backup,1) == EstimOpt.NVarA + EstimOpt.NVarNB + (EstimOpt.NB == 2)
    b0 = B_backup(:);
    if EstimOpt.Display == 1
        disp('Using the starting values from Backup')
    end
elseif isfield(Results_old,'POISS') && isfield(Results_old.POISS,'b0') % starting values provided
    Results_old.POISS.b0_old = Results_old.POISS.b0;
    Results_old.POISS = rmfield(Results_old.POISS,'b0');
    if length(Results_old.POISS.b0_old) ~=  EstimOpt.NVarA + EstimOpt.NVarNB + (EstimOpt.NB == 2)
       if EstimOpt.Display == 1
            disp('WARNING: Incorrect no. of starting values or model specification')
       end
       Results_old.POISS = rmfield(Results_old.POISS,'b0_old');
    else
        b0 = Results_old.POISS.b0_old(:);
    end
end
if  ~exist('b0','var')
    if EstimOpt.Display == 1
        disp('Using linear regression estimates as starting values')
    end
    b0 = regress(log(INPUT.Y+0.1), INPUT.Xa);
    if EstimOpt.NVarNB == 1
        b0 = [b0; zeros(EstimOpt.NVarNB,1)];
    elseif EstimOpt.NVarNB == 2
        b0 = [b0; zeros(EstimOpt.NVarNB,1); 2];
    end
end

%% Optimization Options
if EstimOpt.NB == 2
   EstimOpt.NumGrad = 1; 
end
if isfield(EstimOpt,'BActive')
    EstimOpt.BActive = EstimOpt.BActive(:)';
end

if EstimOpt.ConstVarActive == 1
    if ~isfield(EstimOpt,'BActive') || isempty(EstimOpt.BActive) || sum(EstimOpt.BActive == 0) == 0
        error ('Are there any constraints on model parameters (EstimOpt.ConstVarActive)? Constraints not provided (EstimOpt.BActive).')
    elseif length(b0) ~= length(EstimOpt.BActive)
        error('Check no. of constraints')
    end
    if EstimOpt.Display == 1
        disp(['Initial values: ' mat2str(b0',2)])
        disp(['Parameters with zeros are constrained to their initial values: ' mat2str(EstimOpt.BActive')])
    end
else
    if ~isfield(EstimOpt,'BActive') || isempty(EstimOpt.BActive) || sum(EstimOpt.BActive == 0) == 0
        EstimOpt.BActive = ones(1,length(b0));
        if EstimOpt.Display == 1
            disp(['Initial values: ' mat2str(b0',2)])
        end
    else
        if length(b0) ~= length(EstimOpt.BActive)
            error('Check no. of constraints')
        else
            if EstimOpt.Display == 1
                disp(['Initial values: ' mat2str(b0',2)])
                disp(['Parameters with zeros are constrained to their initial values: ' mat2str(EstimOpt.BActive')])
            end
        end
    end
end

if ((isfield(EstimOpt, 'ConstVarActive') == 1 && EstimOpt.ConstVarActive == 1) || sum(EstimOpt.BActive == 0) > 0) && ~isequal(OptimOpt.GradObj,'on')
    if EstimOpt.Display == 1
        cprintf(rgb('DarkOrange'), 'WARNING: Setting user-supplied gradient on - otherwise parameters'' constraints will be ignored - switch to constrained optimization instead (EstimOpt.ConstVarActive = 1) \n')
    end
    OptimOpt.GradObj = 'on';
end


if (isfield(EstimOpt, 'ConstVarActive') == 0 || EstimOpt.ConstVarActive == 0) && isequal(OptimOpt.Algorithm,'quasi-newton') && isequal(OptimOpt.Hessian,'user-supplied')
    if EstimOpt.Display == 1
        cprintf(rgb('DarkOrange'), 'WARNING: Setting user-supplied Hessian off - quasi-newton algorithm does not use it anyway \n')
    end
    OptimOpt.Hessian = 'off';
end
if EstimOpt.Display == 1
    fprintf('\n')
    cprintf('Opmization algorithm: '); cprintf('*Black',[OptimOpt.Algorithm '\n'])
    if strcmp(OptimOpt.GradObj,'on')
        if EstimOpt.NumGrad == 0
            cprintf('Gradient: '); cprintf('*Black','user-supplied, analytical \n')
        else
            cprintf('Gradient: '); cprintf('*Black',['user-supplied, numerical, ' OptimOpt.FinDiffType '\n'])
        end
    else
        cprintf('Gradient: '); cprintf('*Black',['built-in, ' OptimOpt.FinDiffType '\n'])
    end

    if isequal(OptimOpt.Algorithm,'quasi-newton')
        cprintf('Hessian: '); cprintf('*Black','off, ')
        switch EstimOpt.HessEstFix
            case 0
                cprintf('*Black','retained from optimization \n')
            case 1
                cprintf('*Black','ex-post calculated using BHHH \n')
            case 2
                cprintf('*Black','ex-post calculated using high-precision BHHH \n')
            case 3
                cprintf('*Black','ex-post calculated numerically \n')
            case 4
                cprintf('*Black','ex-post calculated analytically \n')
        end
    else
        if strcmp(OptimOpt.Hessian,'user-supplied')
            cprintf('Hessian: '); cprintf('*Black','user-supplied, BHHH, ')
        else
            cprintf('Hessian: '); cprintf('*Black',['built-in, ' OptimOpt.HessUpdate ', '])
        end
        switch EstimOpt.HessEstFix
            case 0
                cprintf('*Black','retained from optimization \n')
            case 1
                cprintf('*Black','ex-post calculated using BHHH \n')
            case 2
                cprintf('*Black','ex-post calculated using high-precision BHHH \n')
            case 3
                cprintf('*Black','ex-post calculated numerically \n')
            case 4
                cprintf('*Black','ex-post calculated analytically \n')
        end
    end
    fprintf('\n')
end

%% Restucturing Data - no need for now
INPUT.Y = INPUT.Y(INPUT.MissingInd == 0);
INPUT.Xa = INPUT.Xa(INPUT.MissingInd == 0,:);
INPUT.Xnb = INPUT.Xnb(INPUT.MissingInd == 0,:);

%% estimation
LLfun = @(B) LL_poiss_MATlike(INPUT.Y, INPUT.Xa,INPUT.Xnb,INPUT.W, EstimOpt,OptimOpt,B);
if EstimOpt.ConstVarActive == 0 
    
    if EstimOpt.HessEstFix == 0
        [Results.bhat, LL, Results.exitf, Results.output, Results.g, Results.hess] = fminunc(LLfun, b0, OptimOpt);
    else
        [Results.bhat, LL, Results.exitf, Results.output, Results.g] = fminunc(LLfun, b0, OptimOpt);
    end  
    
elseif EstimOpt.ConstVarActive == 1 % equality constraints
        
    EstimOpt.CONS1 = diag(1 - EstimOpt.BActive);
    EstimOpt.CONS1(sum(EstimOpt.CONS1,1)==0,:)=[];
    EstimOpt.CONS2 = zeros(size(EstimOpt.CONS1,1),1);
%     EstimOpt.CONS1 = sparse(EstimOpt.CONS1);
%     EstimOpt.CONS2 = sparse(EstimOpt.CONS2);
    if EstimOpt.HessEstFix == 0
        [Results.bhat, LL, Results.exitf, Results.output, Results.lambda, Results.g, Results.hess] = fmincon(LLfun,b0,[],[],EstimOpt.CONS1,EstimOpt.CONS2,[],[],[],OptimOpt);
    else
        [Results.bhat, LL, Results.exitf, Results.output, Results.lambda, Results.g] = fmincon(LLfun,b0,[],[],EstimOpt.CONS1,EstimOpt.CONS2,[],[],[],OptimOpt);
    end

end

%% Saving results

Results.LL = -LL;
Results.b0_old = b0;


if EstimOpt.HessEstFix == 1
	f = LL_poiss(INPUT.Y,INPUT.Xa,INPUT.Xnb,EstimOpt,Results.bhat);
    Results.jacobian = numdiff(@(B) LL_poiss(INPUT.Y,INPUT.Xa,INPUT.Xnb, EstimOpt,B),f,Results.bhat,isequal(OptimOpt.FinDiffType, 'central'),EstimOpt.BActive);
elseif EstimOpt.HessEstFix == 2
    Results.jacobian = jacobianest(@(B) LL_poiss(INPUT.Y,INPUT.Xa,INPUT.Xnb,EstimOpt,B),Results.bhat);
elseif EstimOpt.HessEstFix == 3
    Results.hess = hessian(@(B) sum(LL_poiss(INPUT.Y,INPUT.Xa,INPUT.Xnb,EstimOpt,B),1), Results.bhat);
end

if sum(EstimOpt.BActive == 0) > 0
    if EstimOpt.HessEstFix == 1 || EstimOpt.HessEstFix == 2
        Results.jacobian = Results.jacobian(:, EstimOpt.BActive == 1);
        Results.hess = Results.jacobian'*Results.jacobian;
    elseif EstimOpt.HessEstFix == 0 || EstimOpt.HessEstFix == 3
        Results.hess = Results.hess(EstimOpt.BActive == 1,EstimOpt.BActive == 1);
    end
    Results.ihess = inv(Results.hess);
    Results.ihess = direcXpnd(Results.ihess,EstimOpt.BActive);
    Results.ihess = direcXpnd(Results.ihess',EstimOpt.BActive);
    
else
    if EstimOpt.HessEstFix == 1 || EstimOpt.HessEstFix == 2
        Results.hess = Results.jacobian'*Results.jacobian;
    end
    Results.ihess = inv(Results.hess);
end

if EstimOpt.RobustStd == 1
    if EstimOpt.NumGrad == 0
        [~, Results.jacobian] = LL_poiss(INPUT.Y,INPUT.Xa,INPUT.Xnb,EstimOpt,Results.bhat);
        Results.jacobian = Results.jacobian.*INPUT.W(:, ones(1,size(Results.jacobian,2)));
    else
        Results.LLdetailed = LL_poiss(INPUT.Y,INPUT.Xa,INPUT.Xnb,EstimOpt,Results.bhat);
        Results.jacobian = numdiff(@(B) INPUT.W.*LL_poiss(INPUT.Y,INPUT.Xa,INPUT.Xnb,EstimOpt,B),Results.LLdetailed,Results.bhat,isequal(OptimOpt.FinDiffType, 'central'),EstimOpt.BActive);
    end
    RobJacob = zeros(EstimOpt.NP, size(Results.jacobian,2));
    RobJacob(1,:) = sum(Results.jacobian(1:EstimOpt.NCTMiss(1),:),1);
    for i = 2:EstimOpt.NP
        RobJacob(i,:) = sum(Results.jacobian(sum(EstimOpt.NCTMiss(1:(i-1)))+1:sum(EstimOpt.NCTMiss(1:i)),:),1);
    end
    RobustHess = RobJacob'*RobJacob;
    Results.ihess = Results.ihess*RobustHess*Results.ihess;
end
Results.std = sqrt(diag(Results.ihess));
if sum(EstimOpt.BActive == 0) > 0
    Results.std(EstimOpt.BActive == 0) = NaN;
end
Results.std(imag(Results.std) ~= 0) = NaN;
Results.R = [Results.bhat , Results.std , pv(Results.bhat , Results.std)];

EstimOpt.Params = length(b0);
if isfield(EstimOpt,'BActive')
    EstimOpt.Params = EstimOpt.Params - sum(EstimOpt.BActive == 0);
end

if isfield(Results_old,'POISS0') && isfield(Results_old.POISS0,'LL')
    Results.stats = [Results.LL; Results_old.POISS0.LL;  1-Results.LL/Results_old.POISS0.LL; NaN; ((2*EstimOpt.Params-2*Results.LL))/EstimOpt.NObs; ((log(EstimOpt.NObs)*EstimOpt.Params-2*Results.LL))/EstimOpt.NObs ;EstimOpt.NObs; EstimOpt.NP; EstimOpt.Params];
else
    Results.stats = NaN(9,1);
end

Results.INPUT = INPUT;
Results.EstimOpt = EstimOpt;
Results.OptimOpt = OptimOpt;

Results.DetailsA = zeros(EstimOpt.NVarA,4);
Results.DetailsA(:,[1 3 4]) = [Results.bhat(1:EstimOpt.NVarA),Results.std(1:EstimOpt.NVarA),pv(Results.bhat(1:EstimOpt.NVarA),Results.std(1:EstimOpt.NVarA))] ;

if EstimOpt.NB == 1
    Results.DetailsNB = zeros(EstimOpt.NVarNB,4);
    Results.DetailsNB(:,[1 3 4]) = [Results.bhat(EstimOpt.NVarA+1:EstimOpt.NVarA+EstimOpt.NVarNB),Results.std(EstimOpt.NVarA+1:EstimOpt.NVarA+EstimOpt.NVarNB),pv(Results.bhat(EstimOpt.NVarA+1:EstimOpt.NVarA+EstimOpt.NVarNB),Results.std(EstimOpt.NVarA+1:EstimOpt.NVarA+EstimOpt.NVarNB))] ;
elseif EstimOpt.NB == 2
    Results.DetailsNB = zeros(EstimOpt.NVarNB+1,4);
    Results.DetailsNB(:,[1 3 4]) = [Results.bhat(EstimOpt.NVarA+1:EstimOpt.NVarA+EstimOpt.NVarNB+1),Results.std(EstimOpt.NVarA+1:EstimOpt.NVarA+EstimOpt.NVarNB+1),pv(Results.bhat(EstimOpt.NVarA+1:EstimOpt.NVarA+EstimOpt.NVarNB+1),Results.std(EstimOpt.NVarA+1:EstimOpt.NVarA+EstimOpt.NVarNB+1))] ;   
end
if EstimOpt.Display == 1
    disp(' ');
    disp('Lambda parameters');
    disp(['var.', blanks(size(char(EstimOpt.NamesA),2)-2) ,'coef.      st.err.  p-value'])
    disp([char(EstimOpt.NamesA) ,blanks(EstimOpt.NVarA)', num2str(Results.DetailsA(:,1),'%8.4f') star_sig(Results.DetailsA(:,4)) num2str(Results.DetailsA(:,3:4),'%8.4f %8.4f')])

    if EstimOpt.NB  > 0 
        disp(' ');
        disp('Overdispersion parameters');
        disp(['var.', blanks(size(char(EstimOpt.NamesNB),2)-2) ,'coef.      st.err.  p-value'])
        disp([char(EstimOpt.NamesNB) ,blanks(EstimOpt.NVarNB + (EstimOpt.NB == 2))', num2str(Results.DetailsNB(:,1),'%8.4f') star_sig(Results.DetailsNB(:,4)) num2str(Results.DetailsNB(:,3:4),'%8.4f %8.4f')])
    end
end
%% Output

% Template1 = {'DetailsA'};
% Template2 = {'DetailsA'};
% Names.DetailsA = EstimOpt.NamesA;
% Heads.DetailsA = {'Covariates of lambda'};
% ST = {};
% 
% if EstimOpt.NB  > 0
%     Template1 = [Template1; 'DetailsNB'];
%     Template2 = [Template2; 'DetailsNB'];
%     Heads.DetailsNB = {'Covariates of alpha (overdispersion parameter)'};
%     Names.DetailsNB = EstimOpt.NamesNB;
%     ST = {'DetailsNB'};
% end
% Head = cell(1,2);
% if EstimOpt.NB == 0
%     Head(1,1) = {'Poisson regression'};
% else
%     Head(1,1) = {'Negative Binomial regression'};
% end
% if EstimOpt.Censored > 0
%     Head(1,2) = {num2str(EstimOpt.Censored, 'with censoring at %2.0f')};
% else
%     Head(1,2) = {'without censoring'};
% end
% 
% Tail = cell(16,2);
% Tail(2,1) = {'Model diagnostics'};
% Tail(3:16,1) = { 'LL at convergence' ; 'LL at constant(s) only'; strcat('McFadden''s pseudo-R',char(178));strcat('Ben-Akiva-Lerman''s pseudo-R',char(178))  ;'AIC/n' ;'BIC/n'; 'n (observations)'; 'r (respondents)';'k (parameters)';' ';'Estimation method';'Optimization method';'Gradient';'Hessian'};
% 
% Tail(3:11,2) = num2cell(Results.stats);
% 
% if any(INPUT.W ~= 1)
%     Tail(13,2) = {'weighted'};
% else
%     Tail(13,2) = {'maximum likelihood'};
% end
% 
% Tail(14,2) = {OptimOpt.Algorithm;};
% 
% if strcmp(OptimOpt.GradObj,'on')
%     if EstimOpt.NumGrad == 0
%         Tail(15,2) = {'user-supplied, analytical'};
%     else
%         Tail(15,2) = {['user-supplied, numerical ',num2str(OptimOpt.FinDiffType)]};
%     end
% else
%     Tail(15,2) = {['built-in, ',num2str(OptimOpt.FinDiffType)]};
%     
% end
% 
% outHessian = [];
% if isequal(OptimOpt.Algorithm,'quasi-newton')
%     outHessian='off, ';
%     switch EstimOpt.HessEstFix
%         case 0
%             outHessian = [outHessian, 'retained from optimization'];
%         case 1
%             outHessian = [outHessian, 'ex-post calculated using BHHH'];
%         case 2
%             outHessian = [outHessian, 'ex-post calculated using high-precision BHHH'];
%         case 3
%             outHessian = [outHessian, 'ex-post calculated numerically'];
%         case 4
%             outHessian = [outHessian, 'ex-post calculated analytically'];
%     end
% else
%     if strcmp(OptimOpt.Hessian,'user-supplied')
%         outHessian = 'user-supplied, BHHH, ';
%         
%     else
%         outHessian = ['built-in, ', num2str(OptimOpt.HessUpdate), ', '];
%     end
%     switch EstimOpt.HessEstFix
%         case 0
%             outHessian = [outHessian, 'retained from optimization'];
%         case 1
%             outHessian = [outHessian, 'ex-post calculated using BHHH'];
%         case 2
%             outHessian = [outHessian, 'ex-post calculated using high-precision BHHH'];
%         case 3
%             outHessian = [outHessian, 'ex-post calculated numerically'];
%         case 4
%             outHessian = [outHessian, 'ex-post calculated analytically'];
%     end
% end
% Tail(16,2) = {outHessian};

%% Tworzenie ResultsOut, drukowanie na ekran i do pliku .xls
Results.Dist = -ones(EstimOpt.NVarA,1);
% if EstimOpt.Display == 1
%     Results.R_out = genOutput(EstimOpt, Results, Head, Tail, Names, Template1, Template2, Heads, ST);
%     fullOrgTemplate = which('template.xls');
%     currFld = pwd;
%     if isfield(EstimOpt,'ProjectName')
%         fullSaveName = strcat(currFld,'\POISS_results_',EstimOpt.ProjectName,'.xls');
%     else
%         fullSaveName = strcat(currFld,'\POISS_results.xls');
%     end
% 
%     copyfile(fullOrgTemplate, 'templateTMP.xls')
%     fullTMPTemplate = which('templateTMP.xls');
%     excel = actxserver('Excel.Application');
%     %     try WbookCheck = excel.Workbooks('fullTMPTemplate');
%     %         Close(WbookCheck)
%     %         catch
%     %              excelWorkbook = excel.Workbooks.Open(fullTMPTemplate);
%     %     end
%     excelWorkbook = excel.Workbooks.Open(fullTMPTemplate);
%     excel.Visible = 1;
%     excel.DisplayAlerts = 0;
%     excelSheets = excel.ActiveWorkbook.Sheets;
%     excelSheet1 = excelSheets.get('Item',1);
%     excelSheet1.Activate;
%     column = size(Results.R_out,2);
%     columnName = [];
%     while column > 0
%         modulo = mod(column - 1,26);
%         columnName = [char(65 + modulo) , columnName];
%         column = floor(((column - modulo) / 26));
%     end
%     rangeE = strcat('A1:',columnName,num2str(size(Results.R_out,1)));
%     excelActivesheetRange = get(excel.Activesheet,'Range',rangeE);
%     excelActivesheetRange.Value = Results.R_out;
%     i = 1;
%     if isfield(EstimOpt,'xlsOverwrite') && EstimOpt.xlsOverwrite == 0
%         while exist(fullSaveName, 'file') == 2
%             if isempty(strfind(fullSaveName, '('))
%                 pos = strfind(fullSaveName, '.xls');
%                 fullSaveName = strcat(fullSaveName(1:pos-1),'(',num2str(i),').xls');
%             else
%                 pos = strfind(fullSaveName, '(');
%                 fullSaveName = strcat(fullSaveName(1:pos),num2str(i),').xls');
%             end
%             i = i+1;
%         end
%     end
%     excelWorkbook.ConflictResolution = 2;
%     SaveAs(excelWorkbook,fullSaveName);
%     excel.DisplayAlerts = 0;
%     excelWorkbook.Saved = 1;
%     Close(excelWorkbook)
%     Quit(excel)
%     delete(excel)
%     delete(fullTMPTemplate)
% end
end