function varargout = build_plasma_state_characteristics(varargin)
% function v = build_plasma_state_characteristics(varargin)
% v contains those kinetic profiles and other quantities which are not 
% in the state and are pre-assigned. At present this is ne,ni,ti,ze (=zeff).
% 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 charnge
        '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};
elseif nargin==3
    model = varargin{1};
    config = varargin{2};
    profiles = varargin{3};
else
    error('must call with 0, 2 or 3 inputs');
end

switch nargin
    case 2
        % initial kinetic profiles
        v = zeros(model.dims.nv,1);

        % mhd placeholder
        v(model.mhd.vind_saw) = 0;
        v(model.ntm.vind) = zeros(size(model.ntm.vind));

        v(model.hmode.vind) = 0;
        
        % Take default values for boundary kappa, main magnetic field,
        % boundary condition for Te.
        v(model.kappa.ind,:) = model.equi.kappa;
        v(model.B0.ind,:) = model.equi.B0;
        v(model.te_BC.ind,:) = model.te.BCval;
        
        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)));
            otherwise
                error('not supported')
        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;
            otherwise
                error('not supported')
        end

        switch model.ze.method
            case 'direct'
                zehat = model.ze.Lamgauss\(init.ze0 + model.rgrid.rhogauss*(init.zee-init.ze0));
                v(model.ze.vind,:) = zehat;
            otherwise
                error('not supported')
        end    
        varargout{1} = v;
    case 3  
        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 = config.numerics.sporder;
        % knots points for spline matrix
        xkts = config.grid.rhogrid;       

        % Spline matrix for quantities with zero gradient at axis
        S = splineval(profiles.rho,xkts,sporder,[1 0 0]);
        Lam = S';
        
        % Get spline coefficients
        v(model.ne.vind,:) = Lam*model.ne.scal\profiles.ne;
        v(model.ni.vind,:) = Lam*model.ni.scal\profiles.ni;
        v(model.ti.vind,:) = Lam*model.ti.scal\profiles.ti;
        v(model.ze.vind,:) = Lam*model.ni.scal\profiles.ze;
        % out
        varargout{1} = v(model.ne.vind,:);
        varargout{2} = v(model.ni.vind,:);
        varargout{3} = v(model.ti.vind,:);
        varargout{4} = v(model.ze.vind,:);
end
return