function [gkii, eqdata] = update_equilibrium(model, params, xkiim1, gkiim1,vk,  subii, it)
% Evaluate the equilibrium
persistent rst LYkii

switch model.equi.source
  case 'liuqe'
    %%% [TODO: put this in a separate function update_equilibrium_LIUQE?] %%%
    % Initialize rst flag
    if subii==1 && it ==1
      rst = 1;
      % rst = 1 first time slice
      % rst = 1 when LIQUE crashed
    end
    
    if rst % If rst==1, restart with standard LIUQE
      prescribe_bf = false; % use standard basis functions
    else % Use modified basis functions
      [pprime, ttprime, rhopsiN] = calc_pp_ttp( xkiim1,gkiim1,vk,model);
      prescribe_bf = true;
    end
    
    % Extract data
    LXk = liuxk(params.equi.LX,it);
    % Solve equilibrium
    if prescribe_bf % Prescribed basis functions
      LYtmp = liut_bf(params.equi.L,LXk,pprime,ttprime);
    else % Standard basis function
      fprintf(' \n \n \n call LIUQE with standard bf \n \n \n')
      LYtmp = liut(params.equi.L, LXk,'pq', params.equi.L.pq,'rst',1,'fixed_pprime_ttprime',0);
    end
    
    % Validate results
    if ~isempty(LYtmp) % Results available: update
      % Eval norms for convergence
      if subii == 1
        eq_state_inc_norm = NaN;
        eq_out_inc_norm = NaN;
      else
        [eq_state_inc_norm,eq_out_inc_norm] = eval_LIUQE_inc_norm(LYtmp, LYkii);
        
      end
      LYkii = LYtmp;
      rst = 0;
    else %if results not available, keep previous LYii and reset LIUQE at next subiteration
      err_str = sprintf('LIUQE failed at subiteration %d and time index %d \n', subii, it);
      fprintf(err_str);
      rst = 1;
      % Reset norms for convergence
      eq_state_inc_norm = NaN;
      eq_out_inc_norm = NaN;
    end
    
    % Store number of iterations and exit flags
    if ~rst
      %flags_eq_transp_debug.N_geometrical_post_processing_iterations_to_convergence(subii, it) = LYkii.zz_subiter;
      eq_info.N_Picard_iterations_to_convergence = LYkii.kit;
      if LYkii.kit >= params.equi.L.P.itert
        eq_info.exceed_Picard_iterations_flag = 1;
      else
        eq_info.exceed_Picard_iterations_flag = 0;          
      end
      
    end
    eq_info.Standard_LIUQE_called_flag = rst;
    eq_info.eq_state_inc_norm = eq_state_inc_norm;
    eq_info.eq_out_inc_norm = eq_out_inc_norm;
    
    
    %  --------------    OUTPUTS ------------------------
    [gkii,~] = LY2RAPTOR(params.equi.L.pQ.^2',LYkii,model.geom,gkiim1);
    eqdata.data_out = LYkii;
    eqdata.data_in = LXk;
    eqdata.eq_info = eq_info;
    
    
    %------------- The other functions have to be fixed still ----------
    
    
    
  case 'liuqe_std'
    error('broken at the moment, need to fix the norms for convergence')
    
    
    
  case 'chease'
    error('broken at the moment')
    gk_new = update_g_from_CHEASE(xkm1,gk,vk,xdot,u,stapk,geopk,trapk,it,dt,model,params);
    % update G for the rest of the foreseeable future
    G(:,it:end) = repmat(gk_new,1,size(G,2)-it+1); %use this
  case 'mock'
    % test case to test interfaces without using equilibrium code.
    [gkii, eqdata] = update_equilibrium_mock(gkiim1);
    eqdata.eq_info.eq_state_inc_norm = NaN;
    eqdata.eq_info.eq_out_inc_norm = NaN;
end

end


function [gkii, eqdata] = update_equilibrium_mock(gkiim1)
% mock version for testing
gkii = gkiim1; % just keep the same gk
eqdata = struct();
end


%%% TO BE MOVED TO SOME update_LIUQE_equilibrium.m
function [norm_state,norm_out] = eval_LIUQE_inc_norm(LYii, LYiim1)

% norm on psi(R,Z) increment
norm_state = norm(LYii.Fx-LYiim1.Fx)/norm(LYii.Fx);


Q0Qnew = LYii.Q0Q;
Q1Qnew = LYii.Q1Q;
Q2Qnew = LYii.Q2Q;
Q3Qnew = LYii.Q3Q;
Q4Qnew = LYii.Q4Q;
TQnew = LYii.TQ;
iqQnew = LYii.iqQ;
psiAnew = LYii.FA;
psiBnew = LYii.FB;

Q0Qold = LYiim1.Q0Q;
Q1Qold = LYiim1.Q1Q;
Q2Qold = LYiim1.Q2Q;
Q3Qold = LYiim1.Q3Q;
Q4Qold = LYiim1.Q4Q;
TQold = LYiim1.TQ;
iqQold = LYiim1.iqQ;
psiAold = LYiim1.FA;
psiBold = LYiim1.FB;

% weights for convergence
wQ0Q = 1/9;
wQ1Q = 1/9;
wQ2Q = 1/9;
wQ3Q = 1/9;
wQ4Q = 1/9;
wTQ = 1/9;
wiqQ = 1/9;
wpsiA = 1/9;
wpsiB = 1/9;

% Masking point if not used in the fitting 
mask = true(size(Q0Qnew));
mask(end) = false;


% output norm
norm_out = wQ0Q*norm(Q0Qnew(mask) - Q0Qold(mask))/norm(Q0Qold(mask))     + wQ1Q*norm(Q1Qnew(mask) - Q1Qold(mask))/norm(Q1Qold(mask))     +...
           wQ2Q*norm(Q2Qnew(mask) - Q2Qold(mask))/norm(Q2Qold(mask))     + wQ3Q*norm(Q3Qnew(mask) - Q3Qold(mask))/norm(Q3Qold(mask))     +...
           wQ4Q*norm(Q4Qnew(mask)-Q4Qold(mask))/norm(Q4Qold(mask))       + wTQ*norm(TQnew(mask)- TQold(mask))/norm(TQold(mask))          +...
           wiqQ*norm(iqQnew(mask) - iqQold(mask))/norm(iqQold(mask))     + wpsiA*norm(psiAnew - psiAold)/norm(psiAold) +...
           wpsiB*norm(psiBnew - psiBold)/norm(psiBold);

end