function [xk, xdotk, gk, uk, stapk, geopk,  trapk, exitflag, res, inewt,   ...
  dftilde_dx, dftilde_dxdot, dftilde_dxk, dftilde_dxkm1, dftilde_du, dftilde_dd, dftilde_dw,...
  eqdata,subiter_info, tran_info] = run_eq_tran_iterations(gkm1, xkm1, xdotkm1, vk,vdotk, dk, U, it, dtk, stapkm1, geopkm1, trapkm1, model, params)

% Iteration between equilibrium and transport. params.numerics.eq_transp_sub_iteration = 1 for no iterations.
% k == time index, ii == subiteration index.
% 'variablek' means at time index 'it'
% 'variablekii' means at time iteration 'it' and eq_transp iteration subii
% Iterated variables alphak will get the name alphakii during equilibrium and transport iteration loop at time it

% the equilbrium results, varagout = { eqresult, debugging flags }

%% Initialization

% Initialize iterated quantities
xkiim1 = xkm1;
gkiim1 = gkm1;
xdotkiim1 = xdotkm1;
stapkiim1 = stapkm1;
geopkiim1 = geopkm1;
trapkiim1 = trapkm1;


nsubii = params.numerics.eq_tran_iter.N_max_iteration;

%% Start loop
for subii = 1:nsubii
  
    %% Update equilibrium
    % Recompute equilibrium
    [gkii, eqdata{subii}] = update_equilibrium(model, params, xkiim1, gkiim1,vk,  subii, it);
    % Relax geometry
    gkii = gkii*params.numerics.eq_tran_iter.relax_alpha + (1-params.numerics.eq_tran_iter.relax_alpha)*gkiim1;
    gdotkii = 0*(gkii-gkiim1)/dtk; % Compute time derivative
    

  %% Run RAPTOR
  % RAPTOR predictive step k-1-> k
  % B.C. INPUT & FEEDBACK CONTROL ISSUES
  % sometimes input is not defined directly but only via a control law
  % in this case the dU_dp from the input parametrization is wrong
  % and needs to be fixed here
  ukii = RAPTOR_controllers(xkiim1,gkii,vk,U,it,model,params);
  
  if params.numerics.usemex
    error('Todo: fix mex call')
  else
    [xkii,xdotkii,dftilde_dx,dftilde_dxdot,dftilde_dxk,dftilde_dxkm1,dftilde_du,dftilde_dd,dftilde_dw,~,...
      stapkii,geopkii,trapkii,res,inewt,exitflag] = ...
      RAPTOR_predictive_step(it,dtk,xkm1, xdotkm1,gkii,gdotkii,vk,vdotk,ukii,dk,stapkiim1,geopkiim1,trapkiim1,model,params);
  end
  
     % Relax xkm xdot, stap, geop,trap
    alpha = params.numerics.eq_tran_iter.relax_alpha;
    xkii  = alpha*xkii + (1-alpha)*xkiim1;
    xdotkii = alpha*xdotkii + (1-alpha)*xdotkiim1;
    
    %% Convergence test and storing
    
    % Store number of iterations and exit flags
    tran_info{subii}.N_Newton_iterations_to_convergence = inewt;
    tran_info{subii}.RAPTOR_exitflag_store = exitflag;
    if subii==1
      tran_info{subii}.RAP_state_inc_norm = NaN;
      tran_info{subii}.RAP_out_inc_norm = NaN;
      isconverged = false;
    else
      tran_info{subii}.RAP_state_inc_norm= eval_RAP_state_inc_norm(xkii, xkiim1, model);
      tran_info{subii}.RAP_out_inc_norm = eval_RAP_out_inc_norm(xkii,gkii,xkiim1,gkiim1, vk,model);
      isconverged = check_all_converged(tran_info{subii},eqdata{subii},subii,it,params);
    end
    
    % Print iteration information
    w2cmdl(subii, it,model, xkii, gkii, tran_info{subii}, eqdata{subii}, params.debug.eq_tran_wtcmdl)
  
    % Exit condition
    doexit = isconverged || (subii == nsubii);
    
  %% Exit or update quantities
  if doexit
    break
  else
    % next subiteration
    xkiim1 = xkii;
    gkiim1 = gkii;
    xdotkiim1 = xdotkii;
    stapkiim1 = stapkii;
    geopkiim1 = geopkii;
    trapkiim1 = trapkii;
  end
end % END equilibrium and transport iteration

% Pass final outputs
xk = xkii;
gk = gkii;
xdotk = xdotkii;
stapk = stapkii;
geopk = geopkii;
trapk = trapkii;
uk = ukii;

% Store debugging flags
subiter_info.N_eq_tran_to_convergence = subii;
if subii == params.numerics.eq_tran_iter.N_max_iteration
    subiter_info.exceed_eq_transp_subiter_flag = 1;
else
    subiter_info.exceed_eq_transp_subiter_flag = 0;
end



end
%---------------------------------------------------------------------------------------------------

function w2cmdl(subii, it,model, xkii, gkii, tran_info, eqdata, eq_tran_wtcmdl)

%% [TODO] this works only for LIUQE now, should be generalized
if mod(subii, eq_tran_wtcmdl) == 0 && strcmp(model.equi.source,'liuqe')
  if subii==1
    fprintf('%-20s: %5s, %9s, %9s, %9s, %9s, %9s, %9s, %9s, %9s, %9s\n','Residual norms','ii','RAPpsiTe','RAPppttp', 'LIUFx' , 'LIUgeom','Ip LY','Ip LX','Ip RA', 'contour_iter', 'Picard iter');
  end
  Ipl_RAPTOR = eval_Ipl(xkii, gkii,[],model,false);
  Ip_RAPTOR=Ipl_RAPTOR(end);
  fprintf('%-20s  %5d: %5.3e, %5.3e, %5.3e, %5.3e, %5.4e, %5.4e, %5.4e, %d, %d \n','',subii, ...
    tran_info.RAP_state_inc_norm, tran_info.RAP_out_inc_norm,...
    eqdata.eq_info.eq_state_inc_norm, eqdata.eq_info.eq_out_inc_norm,...
    eqdata.data_out.Ip, eqdata.data_in.Ip, ...
    Ip_RAPTOR,...
    NaN,  eqdata.data_out.kit);
end

end

function converged = check_all_converged(tran_info,eqdata,subii,it,params)
converged = all([ tran_info.RAP_state_inc_norm,...
  tran_info.RAP_out_inc_norm, ...
  eqdata.eq_info.eq_state_inc_norm, ...
  eqdata.eq_info.eq_out_inc_norm ]...
  < params.numerics.eq_tran_iter.convergence_threshold);

end


function norm_out = eval_RAP_out_inc_norm(xknew,gknew,xkold, gkold, vk, model)

gauss = 1; 

ppnew = eval_pprime(xknew,[],vk,model,gauss);
ttpnew  = eval_ttprime(xknew,gknew,vk,model,gauss);
ppold = eval_pprime(xkold,[],vk,model,gauss);
ttpold = eval_ttprime(xkold, gkold,vk,model,gauss);

% the sum of the weight have to be =1 
wpp = .5;
wttp= .5;

gridsize = size(ppnew);
% masking point
mask = true(gridsize);
mask(end) = false;

% Compute the global norm
norm_out = wpp*norm(ppnew(mask)-ppold(mask))/norm(ppold(mask)) + wttp*norm(ttpnew(mask)-ttpold(mask))/norm(ttpold(mask));

end

function norm_out = eval_RAP_state_inc_norm(xknew, xkold, model)

% Evaluate on gaussian points
gauss = 1;

psinew = eval_psi(xknew, [], [], model, gauss);
tenew = eval_te(xknew, [], [], model, gauss);
psiold = eval_psi(xkold, [], [], model, gauss);
teold = eval_te(xkold, [], [], model, gauss);

% Sum of the weights have to be 1
wpsinorm = 0.5; 
wtenorm = 0.5;

% Global norm
norm_out = wpsinorm*norm(psinew-psiold)/norm(psiold) + wtenorm*norm(tenew-teold)/norm(teold);


end