function [model,params,init,g,v,U] = build_RAPTOR_model(config)
% Builds RAPTOR model based on information supplied in config
% model: model structure, can not change after this call
% params: params structure, may change after this call
% g: geometry information
% v: other prescribed information
% U: example input matrix of the correct size

%% check inputs
check_config(config);

%%
model.tokamak = config.equi.tokamak;
model.realtime = config.realtime;

%% input validity check
if model.realtime
  assert(config.psi.dist.ndist+config.te.dist.ndist>0,'can not have ndist=0 for RT case')
end

%% Grid
ngauss = config.numerics.ngauss;
rho = config.grid.rhogrid';
[rhogauss,wgauss] = GaussGrid(rho,ngauss);

rgrid.rho = config.grid.rhogrid' ; % copy
rgrid.nrho = numel(config.grid.rhogrid);
rgrid.rhogauss = rhogauss;
rgrid.nrhogauss = numel(rhogauss);
rgrid.wgauss = wgauss;
rgrid.ngauss = ngauss;

model.rgrid = rgrid;

%% Spline matrices
% knots points for matrices
xkts = rho;
% spline order
sporder = config.numerics.sporder;

% Spline matrices for quantities with zero gradient at axis
% on rhogauss grid
[S,Sp,Spp,Sppp] = splineval(rhogauss,xkts,sporder,[1 0 0]);
Lamgauss = S'; Lampgauss = Sp'; Lamppgauss = Spp'; Lampppgauss = Sppp';

% on rhotor grid
[S,Sp,Spp,Sppp] = splineval(rho,xkts,sporder,[1 0 0]);
Lam = S'; Lamp = Sp'; Lampp = Spp'; Lamppp = Sppp';
nsp = size(Lam,2);

%% Alternative H mode spline basis
hmode_modeltype_imposed = strcmp(config.hmode.modeltype,'imposed');

if hmode_modeltype_imposed
  % Special splines for prescribed H-mode pedestal
  
  % knot points
  nknots = numel(xkts);
  rhoped = config.hmode.params.rhoped;
  xkts = linspace(0, rhoped, nknots)';
  
  % on rhogauss grid
  [S,Sp,Spp] = splineval_hmode(rhogauss,xkts,sporder);
  Lamgauss_hmode = S'; Lampgauss_hmode = Sp'; Lamppgauss_hmode = Spp';
  
  % on rhotor grid
  [S,Sp,Spp] = splineval_hmode(rho,xkts,sporder);
  Lam_hmode = S'; Lamp_hmode = Sp'; Lampp_hmode = Spp';
end
%% Profile calculations setup

model.psi = setup_psi(config.psi);
model.te = setup_tx(config.te);
model.ti = setup_tx(config.ti);
model.ne = setup_nx(config.ne);
model.ni = setup_nx(config.ni);
model.n1 = setup_nx(config.n1);
model.n2 = setup_nx(config.n2);
model.n3 = setup_nx(config.n3);
model.ze = setup_ze(config.ze);
model.vt = setup_vt(config.vt);
% alternative spline basis Te during H-mode
if hmode_modeltype_imposed
  core_indices = model.rgrid.rhogauss < config.hmode.params.rhoped;
  model.te = scale_Lam(model.te,Lamgauss_hmode,Lampgauss_hmode,Lamppgauss_hmode,Lam_hmode,Lamp_hmode,Lampp_hmode,core_indices,model.te.scal);
  if strcmp(model.ne.method, 'state')
    model.ne = scale_Lam(model.ne,Lamgauss_hmode,Lampgauss_hmode,Lamppgauss_hmode,Lam_hmode,Lamp_hmode,Lampp_hmode,core_indices,model.ne.scal);
  end
  if strcmp(model.ti.method, 'state')
    model.ti = scale_Lam(model.ti,Lamgauss_hmode,Lampgauss_hmode,Lamppgauss_hmode,Lam_hmode,Lamp_hmode,Lampp_hmode,core_indices,model.ti.scal);
  end
end
% MHD configuration
model.ntm   = rmfield(config.ntm,'params');
model.saw   = config.saw.modeltype;

% disturbances for state observer
% parametrization in terms of spline matrices
[model.dist,nd,nw] = build_dist(model);

%% Indices
[model,nx,nv,uind_bc] = setup_indices(model,config);

%% Save dimensions for state x, kinetic profiles v, actuator input u, disturbance d, source matrix
dims.nx = nx;
dims.nv = nv;
dims.nd = nd;
dims.nw = nw;
dims.nsrc = 21; % must match the assignment of src in state_equation.m and unpack in src_unpack.m

%% Copy some parts directly from config
model.envopts = config.envopts; % environment
model.buildopts = config.buildopts;
model.atom = config.atom;
model.controllers = config.controllers;
model.diagnostics = config.diagnostics;

%% HCD module setup
assert(isa(config.echcd,'RAPTORmodule'));
assert(isa(config.nbhcd,'RAPTORmodule'));
assert(isa(config.lhhcd,'RAPTORmodule'));
assert(isa(config.ichcd,'RAPTORmodule'));

[model.echcd,params.echcd] = config.echcd.setup(model);
[model.nbhcd,params.nbhcd] = config.nbhcd.setup(model);
[model.lhhcd,params.lhhcd] = config.lhhcd.setup(model);
[model.ichcd,params.ichcd] = config.ichcd.setup(model);

[model,dims.nu] = setup_u_indices(model,uind_bc);

%% special others ones
model.runaways = config.runaways.modeltype;
model.hmode.modeltype = config.hmode.modeltype;
if hmode_modeltype_imposed
  % Impose fixed h-mode pedestal location
  model.hmode.rhoped = config.hmode.params.rhoped;
end
model.dims = dims;

%% other

model.cvp = config.cvp;

%% transport modules
assert(isa(config.chi_e,'RAPTORmodule'),'config.chi_e must be a RAPTORmodule');
assert(strcmp(config.chi_e.Type,'chi_e'),'config.chi_e is a RAPTORmodule of the wrong type');
assert(isa(config.vpdn_e,'RAPTORmodule'),'config.vpdn_e must be a RAPTORmodule');
assert(strcmp(config.vpdn_e.Type,'vpdn_e'),'config.vpdn_e is a RAPTORmodule of the wrong type');


[model.chi_e,params.chi_e] = config.chi_e.setup(model);
[model.vpdn_e,params.vpdn_e] = config.vpdn_e.setup(model);

model.sn_e = config.sn_e.modeltype;

%% additions
model.atom = config.atom; % to pass atom data to modules

%model.out = out;
model.np = 0; % disabled pending removal

% ped model (tbd)
model.ped = config.ped;

%% build parameter structure
params.tgrid = config.grid.tgrid;
params.numerics = config.numerics;
if (isfield(params.numerics, 'realtime')), params.numerics = rmfield(params.numerics , 'realtime'); end % This is no longer a parameter

assert(isa(config.palpha,'RAPTORmodule'));
assert(isa(config.pbrem,'RAPTORmodule'));
assert(isa(config.prad,'RAPTORmodule'));
assert(isa(config.pei,'RAPTORmodule'));

[model.palpha,params.palpha] =  config.palpha.setup(model);
[model.pbrem,params.pbrem]   =  config.pbrem.setup(model);
[model.prad,params.prad]     =  config.prad.setup(model);
[model.pei,params.pei]       =  config.pei.setup(model);

params.sn_e = config.sn_e.params;
if hmode_modeltype_imposed
  params.hmode.active = config.hmode.params.active;
else
  params.hmode = config.hmode.params;
end
params.saw = config.saw.params;
params.ntm = config.ntm.params;

params.neos = config.neos;

%%
params.runaways = config.runaways.params;
params.debug = config.debug;

%% special cases
% set to use mex by default if built
if model.buildopts.mexify
  params.numerics.usemex = true;
end

%% Real-time simulation parameters
if model.realtime == true
  assert(all(abs(diff(diff(params.tgrid)))<1000*eps),'tgrid must be uniform for realtime')
  params.dt = diff(params.tgrid(1:2));
  %     params = fix_params_realtime(params);
end

%% Geometry
% Get spline matrices and spline coeffcients for equilibrium profiles and
% store in geom and g correspondingly.
[model.geom,g,model.dims.ng, config, params] = build_geom(rho,rhogauss,config, params);

% ------------------------------------------------------
% Build RAPTOR matrices
% Prepare FE  matrices and boundary conditions for state equations
% equations.
if ~hmode_modeltype_imposed
  [model.mats] = build_matrices(model,config,model.geom,g,rho,rhogauss);
else
  [model.mats, model.mats_hmode] = build_matrices(model,config,model.geom,g,rho,rhogauss);
end

%% Add to model the equilibrium in case it was updated in the previous call
model.equi = config.equi; % to pass R0,B0 data
model.equi = rmfield(model.equi,'filenames');
%% Init structure
% Init structure from config from which simple profiles are defined afterwards like ne0, newidth, etc
init = prepare_init(config);

%% Build kinetic profiles
% Construct kinetic profiles according to default settings.
v = build_kinetic_profiles(model,init);
if hmode_modeltype_imposed
  v(model.hmode.vind.te_rhoedge) = config.hmode.params.te_rhoedge;
  v(model.hmode.vind.activation) = config.hmode.params.active;
  v(model.te.BC.vind_value) = config.hmode.params.te_rhoped;
  if strcmp(model.ne.method, 'state')
    v(model.hmode.vind.ne_rhoedge) = config.hmode.params.ne_rhoedge;
    v(model.ne.BC.vind_value) = config.hmode.params.ne_rhoped;
  end
  if strcmp(model.ti.method, 'state')
    v(model.hmode.vind.ti_rhoedge) = config.hmode.params.ti_rhoedge;
    v(model.ti.BC.vind_value) = config.hmode.params.ti_rhoped;
  end
end

%% build example input matrix
U = zeros(model.dims.nu,numel(params.tgrid));
U(1,:) = init.Ip0*ones(1,numel(params.tgrid));

%% Special cases
% set to use mex by default if built
if model.buildopts.mexify
  params.numerics.usemex = true;
end

% optional mexifying
if model.buildopts.mexify
  RAPTOR_mexify(model,params);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NESTED Functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function psi = setup_psi(psiconfig)
    % Psi equation
    psi = psiconfig; % copy some stuff directly
    % scaling
    psi.scal = 1;         psi.eqscal = 0.1;
    
    psi.nsp = nsp;
    % Other auxiliary quantities for psi
    % Psi equation
    % Matrix for computing 1/rho dpsi/drho  (use l'hopitale for central 0/0)
    LamPsipor = [Lampp(1,:);...
      bsxfun(@times,1./rho(2:end),Lamp(2:end,:))];
    LamPsiporgauss = bsxfun(@times,1./rhogauss,Lampgauss);
    
    % matrices to compute d/drho (1/rho) dpsi/drho (needed for diota/drho)
    LamPsipprgauss = bsxfun(@times,1./rhogauss,Lamppgauss) - bsxfun(@times,1./(rhogauss.^2),Lampgauss);
    LamPsippr = (bsxfun(@times,1./rho,Lampp) - bsxfun(@times,1./(rho.^2),Lamp));
    % first few points are wrong due to axis (singularity)
    % don't bother using L'hopital's rule: involves third order derivatives!
    % so just assign to same as next point
    LamPsippr(1,:) = LamPsippr(2,:);
    irr  = find(rhogauss < rho(2)); % find rhogauss rows for first interval: to clip to same value
    LamPsipprgauss(irr,:) = repmat(LamPsippr(2,:),length(irr),1);
    
    
    % evaluate at rhogrid points
    Psi_scal = diag(psi.scal*ones(psi.nsp,1));
    psi.Lam    = Lam    * Psi_scal;
    psi.Lamp   = Lamp   * Psi_scal ;
    psi.Lampp  = Lampp  * Psi_scal ;
    psi.Lamppp = Lamppp * Psi_scal ;
    psi.Lampor = LamPsipor * Psi_scal;
    psi.Lamppr = LamPsippr * Psi_scal;
    
    % to evaluate at rhogauss points
    psi.Lamgauss    = Lamgauss    * Psi_scal;
    psi.Lampgauss   = Lampgauss   * Psi_scal;
    psi.Lamppgauss  = Lamppgauss  * Psi_scal;
    psi.Lampppgauss = Lampppgauss * Psi_scal;
    psi.Lamporgauss = LamPsiporgauss * Psi_scal;
    psi.Lampprgauss = LamPsipprgauss * Psi_scal;
  end

  function tx = setup_tx(txconfig)
    %% ti equation
    tx = txconfig;
    
    tx = init_BC_settings(tx);
    
    [tx.scal,tx.eqscal] = temperature_equation_scalings(tx);
    
    tx.Lam = Lam*tx.scal;
    tx.Lamp = Lamp*tx.scal;
    tx.Lampp = Lampp*tx.scal;
    
    tx.Lamgauss = Lamgauss*tx.scal;
    tx.Lampgauss = Lampgauss*tx.scal;
    
    tx.nsp = size(tx.Lam,2);
  end

  function nx = setup_nx(nxconfig)
    %% n1,n2... equation
    nx = nxconfig;
    nx = init_BC_settings(nx);
    
    [nx.scal,nx.eqscal] = density_equation_scalings(nxconfig);
    nx.Lam = Lam*nx.scal;
    nx.Lamp = Lamp*nx.scal;
    nx.Lampp = Lampp*nx.scal;
    
    nx.Lamgauss = Lamgauss*nx.scal;
    nx.Lampgauss = Lampgauss*nx.scal;
    nx.Lamppgauss = Lamppgauss*nx.scal;
    
    nx.nsp = size(nx.Lam,2);
  end

  function ze = setup_ze(zeconfig)
    %% ze equation
    ze = zeconfig;
    ze.scal = 1; ze.eqscal = 1;
    
    ze.Lam        =   Lam;
    ze.Lamp       =   Lamp;
    ze.Lampp      =   Lampp;
    ze.Lamgauss   =   Lamgauss;
    ze.Lampgauss  =   Lampgauss;
    ze.Lamppgauss =   Lamppgauss;
    ze.nsp = size(ze.Lam,2);
  end

  function vt = setup_vt(vtconfig)
    %% vtor equation
    vt = vtconfig;
    vt.scal = 1e3; vt.eqscal = 1;
    
    vt.Lam        =   Lam;
    vt.Lamp       =   Lamp;
    vt.Lamgauss   =   Lamgauss;
    vt.Lampgauss  =   Lampgauss;
    vt.nsp = size(vt.Lam,2);
  end
end

function [model,nu] = setup_u_indices(model,uind_bc)
nu = uind_bc;
model.echcd.uind = model.echcd.uind+1;
nu = nu + numel(model.echcd.uind);

model.nbhcd.uind = model.nbhcd.uind+nu;
nu = nu + numel(model.nbhcd.uind);

model.lhhcd.uind = model.lhhcd.uind+nu;
nu = nu + numel(model.lhhcd.uind);

model.ichcd.uind = model.ichcd.uind+nu;
nu = nu + numel(model.ichcd.uind);
end


function check_config(config)
% various checks of proper configuration

%% check correct number of qnze species (determined from quasineutrality & zeff)
n_qnze=(sum([...
  strcmpi(config.ne.method,'qnze');
  strcmpi(config.ni.method,'qnze');
  strcmpi(config.n1.method,'qnze');
  strcmpi(config.n2.method,'qnze');
  strcmpi(config.n3.method,'qnze');
  strcmpi(config.ze.method,'qnze');]));

assert(n_qnze == 2 || n_qnze == 0,'must set 2 species or ze to ''qnze'' to satisfy quasineutrality and zeff')

assert(numel(unique([config.atom.Zi,config.atom.Z1,config.atom.Z2,config.atom.Z3]))==4,...
  'atomic charge should be different for all species')

% check that, if we solve for ne via quasineutrality & zeff, we don't also
% scale other profiles from this profile
if strcmpi(config.ne.method,'qnze')
  assert(~any([strcmpi(config.ni.method,'nescal')
    strcmpi(config.n1.method,'nescal')
    strcmpi(config.n2.method,'nescal')
    strcmpi(config.n3.method,'nescal')]),'if config.method.ne=''qnze'', no other profile can be ''nescal''');
end

%idem for ni
if strcmpi(config.ni.method,'qnze')
  assert(~any([strcmpi(config.ni.method,'niscal')
    strcmpi(config.n1.method,'niscal')
    strcmpi(config.n2.method,'niscal')
    strcmpi(config.n3.method,'niscal')]),'if config.method.ni=''qnze'', no other profile can be ''niscal''');
end

% check that if ze = 'qnze', then neos can not be implicit
if strcmp(config.ze.method,'qnze') && (strcmp(config.ne.method,'state') || strcmp(config.ni.method,'state'))
  if config.neos.implicit
    str = 'neos may not be implicit if ze=''qnze'', because dneos/dzeff is not implemented yet\n Set config.neos.implicit=false to fix this';
    warning('checkConfig:neosImplicitWithZeStateDependence',str);
  end
end

%% check for config field structure correctness

% special obsoleteness info
for channels = {'psi','te','ti','ne','ni'}
  assert(~isfield(config.(channels{:}),'BCval'),'BCval usage no longer supported, set BC.defaultValue instead')
end

%% generic check
config_default = construct_defaultconfig(config);

assert(structcmp(config,config_default,-1),'RAPTOR:ConfigStructureInvalid',...
  'config structure does not match default structure')

end

function config_default = construct_defaultconfig(config)

module_fields = {'chi_e','vpdn_e'};
model_fields = {'echcd','nbhcd','ichcd','lhhcd','hmode','saw','prad','runaways','sn_e'};
special_fields = {'ntm','controllers'};

% construct default config with the same modules as the input config structure
config_default = RAPTOR_config;
fields = fieldnames(config);

for ifield=1:numel(fields)
  myfield = fields{ifield};
  
  if ismember(myfield,special_fields)
    config_default.(myfield) = config.(myfield); % don't test these ones at all yet
    continue
  end
  
  if isfield(config.(myfield),'params')
    % if this field corresponds to an 'old' module setup
    %fprintf('old module %s\n',myfield);
    if ismember(myfield,module_fields)
      module_name = config.(myfield).model.name;
    elseif ismember(myfield,model_fields)
      module_name = config.(myfield).modeltype;
    else
      module_name = '';
    end
    config_default.(myfield) = eval(sprintf('%s(''%s'')',myfield,module_name));
  elseif isa(config.(myfield),'RAPTORmodule')
    moduleName = config.(myfield).Name;
    config_default.(myfield) = RAPTORmodule(moduleName);
  else
    assert(isfield(config_default,myfield),'RAPTOR:ConfigFieldInvalid',...
      'Field ''%s'' is invalid: not part of default RAPTOR config',myfield)
  end
end
end

function [model,nx,nv,nu] = setup_indices(model,config)
%% Indices
vind = 0; xind = 0; uind = 0;
for myeqs = {'psi','te','ti','ne','ni','ze','n1','n2','n3','vt'}
  myeq = myeqs{:};
  if strcmp(model.(myeq).method,'state')
    model.(myeq).xind = xind + (1:model.(myeq).nsp);
    xind = xind + model.(myeq).nsp;
    model.(myeq).vind = [];
    
    switch model.(myeq).BC.valueSource
      case 'v'
        vind = vind+1;
        model.(myeq).BC.vind_value = vind;
      case 'u'
        uind = uind+1;
        model.(myeq).BC.uind = uind;
      otherwise
        % nothing
    end
    
    switch model.(myeq).BC.rhoSource
      case 'v'
        vind = vind+1;
        model.(myeq).BC.vind_rho = vind;
      otherwise
        % nothing
    end
    
  else
    model.(myeq).vind = vind + (1:model.(myeq).nsp);
    vind = vind + model.(myeq).nsp;
    model.(myeq).xind = [];
  end
end

% hmode index: next one in v0
model.hmode.vind.activation = vind + 1; % always present
vind=vind+1;
switch config.hmode.modeltype
  % Specific vind handling for various options of H-mode module
  case 'imposed'
    model.hmode.vind.te_rhoedge = vind+1;
    vind=vind+1;
    if strcmp(model.ne.method,'state')
      model.hmode.vind.ne_rhoedge = vind+1;
      vind=vind+1;
    end
    if strcmp(model.ti.method,'state')
      model.hmode.vind.ti_rhoedge = vind+1;
      vind=vind+1;
    end
  case 'gradient'
    model.hmode.vind.rhoped = vind + 1;
    vind=vind+1;
    
    model.hmode.vind.gteped = vind + 1;
    vind = vind+1;
    
    model.hmode.vind.gtiped = vind + 1;
    vind = vind+1;
    
    model.hmode.vind.gneped = vind + 1;
    vind = vind+1;
end

model.mhd.vind_saw   = vind + 1; % add to vk
vind=vind+1;

% special case for NTM: may be in state, but need seed island trigger via v
% in any case.
if strcmp(model.ntm.method,'state')
  model.ntm.xind = xind + (1:model.ntm.n_ntm);
  xind = xind + model.ntm.n_ntm;
  model.ntm.vind = vind + (1:model.ntm.n_ntm);
else
  model.ntm.xind = [];
  model.ntm.vind = vind + (1:model.ntm.n_ntm);
end
vind = vind + model.ntm.n_ntm; % update counter

nv = vind;
nx = xind;
nu = uind;

end

function tx = init_BC_settings(tx)
% Make BC settings compatible with solver settings
if strcmp(tx.method,'state')
  if strcmp(tx.BC.rhoSource,'v');
    % remove rho value in structure
    tx.BC = rmfield(tx.BC,'rho');
  end
else    % needs no BC
  tx.BC = 'n/a';
end

end

function modelxx = scale_Lam(modelxx,Lamgauss_hmode,Lampgauss_hmode,Lamppgauss_hmode,Lam_hmode,Lamp_hmode,Lampp_hmode,core_indices,scal)
modelxx.Lamgauss_hmode = Lamgauss_hmode*scal;
modelxx.Lampgauss_hmode = Lampgauss_hmode*scal;
modelxx.Lamppgauss_hmode = Lamppgauss_hmode*scal;
modelxx.Lam_hmode = Lam_hmode*scal;
modelxx.Lamp_hmode = Lamp_hmode*scal;
modelxx.Lampp_hmode = Lampp_hmode*scal;
modelxx.T_LH = modelxx.Lamgauss_hmode(core_indices, 1:end-1) \ modelxx.Lamgauss(core_indices, :);
modelxx.T_HL = modelxx.Lamgauss \ modelxx.Lamgauss_hmode;
end