function varargout = build_kinetic_profiles(varargin)
% function v = build_kinetic_profiles(varargin)
%  v = build_kinetic_profiles(model,init); % specify via init
%  v = build_kinetic_profiles(model,params,profiles); % specify directly
%
% v contains those kinetic profiles and other quantities which are not
% solved for (i.e. not part of the state) are pre-assigned. 
% At present this is ne,ni,ti,ze (=zeff) and sometimes impurity profiles.
% and also H-mode, NTM, sawtooth crash information.
% If nargin==2 then plasma params defined via init structure and varargout==v.
% If nargin==3 then plasma params defined via profiles structure which has
% to have the following fields: electron density .ne, ion density .ni, ion
% temperature .ti, effective charge .ze, normalized radial grid for the
% profiles .rho. There are the following output parameters:
% v(model.ne.vind,:), v(model.ni.vind,:), v(model.ti.vind,:),
% v(model.ze.vind,:).

if nargin == 0,
    %% DEFAULT PARAMETERS
    module_params = struct(...
        'ne', [], ... % electron density
        'ni',[],  ... % ion density
        'ti',[],... % ion temperature
        'ze',[],... % effective charge
        'n1',[],... % minor species No 1
        'n2',[],... % minor species No 2
        'n3',[],... % minor species No 3
        'rho',[]... % normalized radial grid for the profiles
        );
    
    varargout{1} = module_params;
    return % empty call, probably to get default structures
elseif nargin==2; % change this depending on number of inputs to module
    model = varargin{1};
    init = varargin{2};
    srcmode = 'init';
elseif nargin==3; % change this depending on number of inputs to module
    model = varargin{1};
    params = varargin{2};
    profiles = varargin{3};
    srcmode = 'profiles';
else
    error('must call with 0 or 2 or 3 inputs');
end


switch srcmode
    case 'init'
        % initial kinetic profiles
        v = zeros(model.dims.nv,1);

        switch model.ti.method
            case 'direct'
                v(model.ti.vind,:) = model.ti.Lamgauss\(init.ti0*exp(-model.rgrid.rhogauss.^2/(init.tiwidth^2)));
            case 'tescal'
                % scaling profile
                tiscal = init.ti0scal + model.rgrid.rhogauss*(init.tiescal-init.ti0scal);
                v(model.ti.vind,:) = model.ti.Lamgauss\tiscal;
        end

        switch model.ne.method
            case 'direct'
                v(model.ne.vind,:) = model.ne.Lamgauss\(init.ne0*exp(-model.rgrid.rhogauss.^2/(init.newidth^2)));
        end

        switch model.ni.method
            case 'direct'
                v(model.ni.vind,:) = model.ni.Lamgauss\(init.ni0*exp(-model.rgrid.rhogauss.^2/(init.niwidth^2)));
            case 'nescal'
                % scaling profile
                niscal = init.ni0scal + model.rgrid.rhogauss*(init.niescal-init.ni0scal);
                v(model.ni.vind,:) = model.ni.Lamgauss\niscal;
            case 'calc'
                % do nothing
            case 'qnze'
                % calculated from quasineutrality and Zeff
            otherwise
                error('not supported')
        end

        switch model.ze.method
            case 'direct'
                vthat = model.ze.Lamgauss\(init.ze0 + model.rgrid.rhogauss*(init.zee-init.ze0));
                v(model.ze.vind,:) = vthat;
            case 'qnze'
                % calculated from quasineutrality and Zeff
            otherwise
                error('not supported')
        end   
        
        % for ion species n1,n2,n3...
        for ii=1:3
            nxmodel = get_nxmodel(model,ii);
            switch nxmodel.method
                case {'nescal','niscal'}
                    sc = get_nxscal(init,ii);
                    scal = sc*ones(model.rgrid.nrhogauss,1);
                    v(nxmodel.vind,:) = nxmodel.Lamgauss\scal;
                case 'qnze'
                    % calculated from quasineutrality and Zeff
                case 'direct'
                    v(nxmodel.vind) = 0; % zero by default when directly prescribed
                otherwise
                    error('not supported')
            end
        end
        
        switch model.vt.method
            case 'direct'
                vthat = model.vt.Lamgauss\(init.vt0 + model.rgrid.rhogauss.^2*(init.vte-init.vt0));
                v(model.vt.vind,:) = vthat;
            otherwise
                error('not supported')
        end 
        
        v = set_v_BC(v,model);
        
        varargout{1} = v;
    case 'profiles'  
        v = zeros(model.dims.nv,size(profiles.ne,2));
        
        % Construct spline matrix to get spline coefficients for kinetic
        % profiles and store them in v.
        % spline order
        sporder = params.numerics.sporder;
        % knots points for spline matrix
        xkts = model.rgrid.rho;       

        % Spline matrix for quantities with zero gradient at axis
        S = splineval(profiles.rho,xkts,sporder,[1 0 0]);
        Lam = S';
        
	assert(size(Lam,1)>=size(Lam,2),['insufficient radial grid', ...
		    ' points on prescribed profile grid w.r.t.', ...
		     ' RAPTOR grid. Increase resolution of input', ...
		      ' grid or decrease resolution of RAPTOR grid']);
        % Get spline coefficients
        v(model.ne.vind,:) = model.ne.scal*Lam\profiles.ne;
        v(model.ni.vind,:) = model.ni.scal*Lam\profiles.ni;
        v(model.ti.vind,:) = model.ti.scal*Lam\profiles.ti;
        v(model.ze.vind,:) = model.ni.scal*Lam\profiles.ze;  
        
        % Set boundary conditions
        v = set_v_BC(v,model);
        
        varargout{1} = v;

end
return

function  nxmodel = get_nxmodel(model,ii)
% # codegen
coder.extrinsic('warning');

if ii==1
    nxmodel = model.n1;
elseif ii==2
    nxmodel = model.n2;
elseif ii==3
    nxmodel = model.n3;
else
    warning('wrong model index')
    nxmodel = model.n1; % dummy for codegen
end
return
    
function scal = get_nxscal(init,ii)
% # codegen
coder.extrinsic('warning');

if ii==1
    scal = init.n1scal;
elseif ii==2
    scal = init.n2scal;
elseif ii==3
    scal = init.n3scal;
else
    warning('wrong model index')
    scal = 1; % dummy for codegen
end
return

function v = set_v_BC(v,model)

defaultRho = 1;

% set boundary conditions from model values
for channels = {'te','ti','ne','ni'}
    mychannel = channels{:};
    if ~strcmp(model.(mychannel).BC,'n/a')
        if strcmp(model.(mychannel).BC.valueSource,'v')
            v(model.(mychannel).BC.vind_value) = model.(mychannel).BC.defaultValue;
        end
        if strcmp(model.(mychannel).BC.rhoSource,'v')
            v(model.(mychannel).BC.vind_rho) = defaultRho;
        end

    end
end
return
