function x0 = RAPTOR_initial_conditions(model,init,g,v)% #codegen
% x0 = RAPTOR_initial_conditions(varargin)
% Generates initial condition for state (x0).

coder.extrinsic('fprintf');
coder.extrinsic('error');
coder.extrinsic('warning');

assert(iscolumn(g),'input g must be a vector, e.g. g(:,1)');
assert(iscolumn(v),'input v must be a vector, e.g. g(:,1)');

if isfield(init,'new')
    warning('Field init.new is not allowed. Update your init structure w.r.t. the new default in config');
end

%%% INITIAL CONDITIONS FOR STATE

%% Psi
%%% assign via jtor0
if init.jpow > 0
    jtor0 = 1-(model.rgrid.rhogauss).^init.jpow; % polynomial current density profile
else
    jtor0 = (1-model.rgrid.rhogauss).^(-init.jpow); % polynomial current density profile
end
% % if init.j_peaked == false
% %   jtor0(model.rgrid.rhogauss<.3)=linspace(jtor0(find(model.rgrid.rhogauss>=.3,1,'first'))/2, jtor0(find(model.rgrid.rhogauss>=.3,1,'first')), sum(model.rgrid.rhogauss<.3));
% % end

xdummy = zeros(model.psi.nsp+model.te.nsp,1);
[~,djtor_dx] = eval_jtor(xdummy,g,v,model,true);
psihat0 = [djtor_dx(:,model.psi.xind);model.psi.Lam(end,:)]\[jtor0;0]; % Psihat for this jtor

Ip0 = init.Ip0;
% correct initial condition for consistency with initial input
if strcmp(model.psi.BC.type,'PsiOH');
    % must match boundary condition, let PsiOH(0) = 0; (by choice)
    %Lext*Ip(t0,rhoe) + psi(t0,rhoe) = PsiOH(t0) = 0;
    psihat0 = psihat0 * Ip0/(model.out.CIp(end,:)*psihat0); % set right Ip
    Psiedge0 = -profiles.eq_init.Lext * model.out.CIp(end,:)*psihat0;
    psihat0 = psihat0-psihat0(end)+Psiedge0; % put to zero at edge
elseif strcmp(model.psi.BC.type,'Ip')
    % scale current of initial Psihat to match U(1)
    Ipl = eval_Ipl([psihat0;0*psihat0],g,v,model,false);
    psihat0 = psihat0 * Ip0/Ipl(end);
    psihat0 = psihat0 + init.psib0; % add offset for desired boundary value
else
    error('unknown Psi BC type')
end


%% Te
% Initial profile for te
use_core_spline_base = hmode_spline_base_active(v,model);
te0 = init.te0*exp(-model.rgrid.rhogauss.^2/(init.tewidth^2));
tehat0 = build_initial_state(te0, model, 'te', use_core_spline_base);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Optional states
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initial profile for Ti
if strcmp(model.ti.method,'state')
    ti0 = init.ti0*exp(-model.rgrid.rhogauss.^2/(init.tiwidth^2));
    tihat0 = build_initial_state(ti0, model, 'ti', use_core_spline_base);
else
    tihat0 = [];
end

% Initial profile for ne
if strcmp(model.ne.method,'state')
    ne0 = init.ne0*exp(-model.rgrid.rhogauss.^2/(init.newidth^2));
    nehat0 = build_initial_state(ne0, model, 'ne', use_core_spline_base);
else
    nehat0 = [];
end

% Initial profile for ni
if strcmp(model.ni.method,'state')
    nihat0 = model.ne.Lamgauss\(init.ni0*exp(-model.rgrid.rhogauss.^2/(init.niwidth^2)));
    ni_edge = (model.ni.Lam(end,:)*nihat0);
    if  ni_edge ~= model.ni.BCval
        nihat0 = model.ni.Lamgauss\(model.ni.Lamgauss*nihat0+((model.rgrid.rhogauss.^2)*(model.ni.BCval-ni_edge)));
    end
else
    nihat0 = [];
end

%% NTM
% Initial NTM island width states are zero
if strcmp(model.ntm.method,'state')
    ntm0 = zeros(model.ntm.n_ntm,1);
else
    ntm0 = [];
end

%% checks
if ~use_core_spline_base
    te_bc = model.te.Lamgauss(end,:)*tehat0;
else
    core_indices = model.rgrid.rhogauss < model.hmode.rhoped;
    te_bc = model.te.Lamgauss_hmode(core_indices(end),1:end-1)*tehat0;
end
if norm(te_bc)<0; error('Zero Te_edge not compatible with nonzero j_edge'); end

%% assign
xind = [model.psi.xind,model.te.xind,model.ti.xind,model.ne.xind,model.ni.xind,model.ntm.xind];
x0 = zeros(numel(xind),1);
x0([model.psi.xind,model.te.xind,model.ti.xind,model.ne.xind,model.ni.xind,model.ntm.xind],1) = ...
    [psihat0;tehat0;tihat0;nehat0;nihat0;ntm0];

return

function xxhat0 = build_initial_state(xx0, model, xxID, use_core_spline_base)
switch xxID
  case 'te'
    modelxx = model.te;
  case 'ti'
    modelxx = model.ti;
  case 'ne'
    modelxx = model.ne;
end
if ~use_core_spline_base
    xxhat0 = modelxx.Lamgauss\xx0;
else
    % initial profile fitted by compressed splines on reduced domain [0, rhoped]
    core_indices = model.rgrid.rhogauss < model.hmode.rhoped;
    xxhat0 = modelxx.Lamgauss_hmode(core_indices,1:end-1)\xx0(core_indices);
end
return
