function [x, xdot, u, out_data, equil_profiles, equil_psiRZmap, res_vec_j, res_vec_p, res_vec_ip] = RAPTOR_CHEASE_iterator(config, model, params, init, g, v, u, ...
         paramz, inputs, options)
% Find stationary state radial profiles with consitent MHD equilibrium.
% --
% CHEASE, a fixed boundary Grad-Shafranov solver, updates the MHD equilibrium,
% with current and pressure profiles imposed by RAPTOR.
% RAPTOR, a core transport solver, updates the radial profiles, 
% with geometry imposed by CHEASE.
% --
% Iteration until the radial profiles of the RAPTOR solution sufficiently match the
% radial profiles on which the underlying CHEASE solution is based.
% --
% The options structure contains
% -- options.free_fns -- with options: 'jpar_p_rhotor', 'jpar_p_rhopol', 'ttprime_pprime_rhopol'
% Selects current and pressure related profiles, defining the free functions
% p'(psi) and TT'(psi) in the Grad-Shafranov equation;
% furthermore, specifies whether these profiles are provided on rhopol or
% rhotor.
% -- options.res_prof1 with options: 'ttprime', 'q', 'jpar'
% -- options.res_prof2 with options: 'pprime', 'p'
% -- options.rm_eqfiles -- (boolean)
% -- options.verbosity -- with options 0 (no report at all), 1 (RAPTORss global
% information reported), 2 (RAPTORss iteration step information reported + 
% plots RAPTORss-CHEASE convergence), 3 (plots RAPTORss)

% Numerics
nr_iter = 3; res_tol = 1e-3;
grelax = 1;

out_data = cell(1, nr_iter+1);
equil_profiles = cell(1, nr_iter+1);
equil_psiRZmap = cell(1, nr_iter+1);

% Initial MHD equilibrium
addpath(fullfile(model.envopts.RAPTOR_path, 'equils', 'chease'));
load(config.equi.filenames{1});

raptor_grid = model.rgrid.rho;

equil_profiles{1} = get_equil_profiles(raptor_grid, chease_data);
equil_psiRZmap{1} = get_equil_psiRZmap(chease_data);

% Solve core transport: RAPTOR
if isfield(init,'x0')
  [x, xdot, u, out_data{1}] = stationary_state_solver(model, params, init, g, v, u, ...
                            paramz, inputs, options.verbosity, init.x0);
else
  [x, xdot, u, out_data{1}] = stationary_state_solver(model, params, init, g, v, u, ...
                            paramz, inputs, options.verbosity);
end
index = 1;
if options.verbosity >= 2
  plot_trap_and_equil_profiles(raptor_grid, equil_profiles, out_data, 1);
%   plot_geom_profiles(raptor_grid, model, g, 2);
end

nout_file = '';
res_vec_j = zeros(1, nr_iter);
res_vec_p = zeros(1, nr_iter);
res_vec_ip =  zeros(1, nr_iter);

while index <= nr_iter
  
  [res_j, res_p] = eval_res(equil_profiles, out_data, options.res_prof1, options.res_prof2, index);
  res_vec_j(index) = res_j;
  res_vec_p(index) = res_p;
  res_vec_ip(index) = chease_data.eqdsk.ip;
  
  if res_j < res_tol && res_p < res_tol
    if options.verbosity > 0
      fprintf('RAPTORss and CHEASE converged in %u iterations\n', index-1);
    end
    out_data = out_data(~cellfun('isempty', out_data));
    equil_profiles = equil_profiles(~cellfun('isempty', equil_profiles));
    res_vec_j = res_vec_j(res_vec_j~=0);
    res_vec_p = res_vec_p(res_vec_p~=0);
    res_vec_ip = res_vec_ip(res_vec_ip~=0);
    break
  else
    % Solve MHD equilibrium
    if index == 1
      B0 = chease_data.eqdsk.b0;
      R0 = chease_data.eqdsk.r0;
      r_contour = chease_data.eqdsk.rplas_equi;
      z_contour = chease_data.eqdsk.zplas_equi;
      RZ_contours =  [r_contour, z_contour]./R0;
    end
    ip = u(1);
    [g_upd, chease_data, nout_file] = find_consistent_CHEASE_equilibrium(out_data{index}, ip, config, model, ...
                                      B0, R0, RZ_contours, nout_file, options.free_fns, index);
    equil_profiles{index+1} = get_equil_profiles(raptor_grid, chease_data);
    equil_psiRZmap{index+1} = get_equil_psiRZmap(chease_data);
    g = grelax*g_upd + (1-grelax)*g;
    
    % Solve core transport
    if isfield(init,'x0')
      [x, xdot, u, out_data{index+1}] = stationary_state_solver(model, params, init, g, v, u, ...
        paramz, inputs, options.verbosity, init.x0);
    else
      [x, xdot, u, out_data{index+1}] = stationary_state_solver(model, params, init, g, v, u, ...
        paramz, inputs, options.verbosity);
    end
    
    index = index+1;
    
    if options.verbosity > 1
      plot_trap_and_equil_profiles(raptor_grid, equil_profiles, out_data, index);
%       plot_geom_profiles(raptor_grid, model, g, index);
    end
    
  end
  
  
end

% automatically remove EXPEQ, EQDSK, namelist files
if options.rm_eqfiles
  rmdir(fullfile(model.envopts.RAPTOR_path, 'data', 'CHEASE_updt_EQ'), 's');
  [~,uname] = unix('whoami < /dev/null');
  rmdir(strcat('/tmp/', uname, '/chease*'), 's'); 
  fprintf(strcat('chease* files in /tmp/', uname, ' removed \n'));
end
end

function [error_j, error_p] = eval_res(equil_profiles, out_data, res_prof1, res_prof2, index)
switch res_prof1
  case 'ttprime'
    j_scal = norm(out_data{index}.ttprime * (2*pi));
    error_j = norm(equil_profiles{index}.ttprime - out_data{index}.ttprime * (2*pi))/j_scal;
  case 'q'
    j_scal = norm(out_data{index}.q);
    error_j = norm(equil_profiles{index}.q - out_data{index}.q)/j_scal;
  case 'jpar'
    j_scal = norm(out_data{index}.jpar);
    error_j = norm(equil_profiles{index}.jpar - out_data{index}.jpar)/j_scal;
end
switch res_prof2
  case 'p'
    p_scal = norm(out_data{index}.p);
    error_p = norm(equil_profiles{index}.p - out_data{index}.p)/p_scal;
  case 'pprime'
    p_scal = norm(out_data{index}.pprime * (2*pi));
    error_p = norm(equil_profiles{index}.pprime - out_data{index}.pprime * (2*pi))/p_scal;
end
end

function plot_trap_and_equil_profiles(raptor_grid, equil_profiles, out_data, index)
% figure(2*index-1)
figure(index)
subplot(3,2,1); plot(raptor_grid, equil_profiles{index}.pprime, 'b', raptor_grid, out_data{index}.pprime * (2*pi), 'r--');
legend({'chease', 'raptor'}); title('pprime')
subplot(3,2,2); plot(raptor_grid, equil_profiles{index}.ttprime, 'b', raptor_grid, out_data{index}.ttprime * (2*pi), 'r--');
legend({'chease', 'raptor'}); title('ttprime')
subplot(3,2,3); plot(raptor_grid, equil_profiles{index}.p, 'b', raptor_grid, out_data{index}.p, 'r--');
legend({'chease', 'raptor'}); title('p')
subplot(3,2,4); plot(raptor_grid, equil_profiles{index}.q, 'b', raptor_grid, out_data{index}.q, 'r--');
legend({'chease', 'raptor'}); title('q')
subplot(3,2,5); plot(raptor_grid, equil_profiles{index}.jpar, 'b', raptor_grid, out_data{index}.jpar, 'r--');
legend({'chease', 'raptor'}); title('jpar')
% 2pi factor because of difference in psi definition RAPTOR and CHEASE
end

% function plot_geom_profiles(raptor_grid, model, g, index)
% figure(2*index)
% subplot(3,2,1); plot(raptor_grid, model.geom.LamG*g(model.geom.ind_gg0)); title('g0');
% subplot(3,2,2); plot(raptor_grid, model.geom.LamG*g(model.geom.ind_gg1)); title('g1');
% subplot(3,2,3); plot(raptor_grid, model.geom.LamG*g(model.geom.ind_gg2)); title('g2');
% subplot(3,2,4); plot(raptor_grid, model.geom.LamG*g(model.geom.ind_gg3)); title('g3');
% subplot(3,2,5); plot(raptor_grid, model.geom.LamG*g(model.geom.ind_F)); title('F');
% subplot(3,2,6); plot(raptor_grid, model.geom.LamG*g(model.geom.ind_Vp)); title('Vp');
% end

function equil_profiles = get_equil_profiles(raptor_grid, chease_data)

chease_rhonrm_eqdsk = chease_data.eqdsk.rhotor_norm;
chease_rhonrm_cols = chease_data.profiledata.rho_tor_norm.data;

chease_p = chease_data.eqdsk.p;
chease_q = chease_data.eqdsk.q;
chease_pprime = chease_data.eqdsk.pprime;
chease_ttprime = chease_data.eqdsk.FFprime;
B0 = chease_data.eqdsk.b0; R0 = chease_data.eqdsk.r0; mu0 = 4e-7*pi;
chease_jpar = chease_data.profiledata.j_par_eq_av_j_dot_b_av_over_b0.data * (B0/(mu0*R0));

equil_profiles.p = interpos(chease_rhonrm_eqdsk, chease_p, raptor_grid);
equil_profiles.q = interpos(chease_rhonrm_eqdsk, chease_q, raptor_grid);
equil_profiles.pprime = interpos(chease_rhonrm_eqdsk, chease_pprime, raptor_grid);
equil_profiles.ttprime = interpos(chease_rhonrm_eqdsk, chease_ttprime, raptor_grid);
equil_profiles.jpar = interpos(chease_rhonrm_cols, chease_jpar, raptor_grid);
end

function equil_psiRZmap = get_equil_psiRZmap(chease_data)
equil_psiRZmap.rmesh = chease_data.eqdsk.rmesh;
equil_psiRZmap.zmesh = chease_data.eqdsk.zmesh;
equil_psiRZmap.psi = chease_data.eqdsk.psi;
equil_psiRZmap.rplas = chease_data.eqdsk.rplas;
equil_psiRZmap.zplas = chease_data.eqdsk.zplas;
equil_psiRZmap.raxis = chease_data.eqdsk.raxis;
equil_psiRZmap.zaxis = chease_data.eqdsk.zaxis;
equil_psiRZmap.psiaxis = chease_data.eqdsk.psiaxis;
equil_psiRZmap.psiedge = chease_data.eqdsk.psiedge;
end