function varargout = saw(varargin) %#codegen

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Input processing
if nargin == 0,
    %% DEFAULT PARAMETERS
    module_params = struct(...
    'modeltype' , 'none' , ... % default module
    'params',struct('active',false)); % load default params
    varargout{1} = module_params;
    return %empty call, probably to get default structures
elseif nargin == 1
    % get default parameters for each case
    def_params.modeltype = varargin{1}; 
    switch def_params.modeltype
        case 'none'
            def_params.params = struct('active',false);
        case 'dummy'
            def_params.params = dummy;
        case 'porcelli'
            def_params.params = porcelli;
        otherwise
            error('invalid sawtooth model type ''%s''',def_params.modeltype)
    end
    varargout{1} = def_params; % return appropriate parameter structure
    return
elseif nargin==7;
    x = varargin{1};
    g = varargin{2};
    v = varargin{3};
    u = varargin{4};
    it = varargin{5};
    model  = varargin{6};
    saw_params = varargin{7}; % distribute inputs
else
    error('must call with 0 or 7 inputs');
end

%% Compute output depending on model
switch model.saw 
    case 'none'
        xnew = x;
    case 'dummy'
        xnew = dummy(x,g,v,u,it,model,saw_params);
    case 'porcelli'
        xnew = porcelli(x,g,v,u,it,model,saw_params);
    otherwise
        error('unknown sawtooth model')
end


varargout{1} = xnew;

return

%% % dummy (example) module
function varargout = dummy(varargin)

if nargin == 0,
    module_params = struct(...
        'qcrit',1);
    varargout{1} = module_params;
    return %empty call, probably to get default structures
elseif nargin==7;
    x = varargin{1};
    g = varargin{2};
    v = varargin{3};
    u = varargin{4};
    it = varargin{5};
    model  = varargin{6};
    saw_params = varargin{7}; % distribute inputs
else
    error('must call with 0 or 7 inputs');
end

xnew = x; % init;
tegauss = model.te.LamTegauss*x(model.te.xind);
iotagauss = model.out.Ciotagauss*x(model.psi.xind);
qgauss = 1./iotagauss; qmin=min(qgauss);
ne_rho = model.te.LamTegauss * v(model.ne.vind);


if qmin<saw_params.qcrit
    %%
    % crash Te profile
    tegaussnew=tegauss;
    tegaussnew(model.rgrid.rhogauss<0.3) = min(tegauss(model.rgrid.rhogauss<0.3));
    
    % crash q profile
    qgaussnew=qgauss; % init
    qgaussnew(qgauss<1) = 1;
    % reassign state
    subplot(211)
    plot(model.rgrid.rhogauss,[tegauss,tegaussnew]);
    subplot(212)
    plot(model.rgrid.rhogauss,[qgauss,qgaussnew]);
    
    psib = model.psi.Lam(end,:)*x(model.psi.xind); % keep psib fixed
    xnew(model.psi.xind) = [model.out.Ciotagauss;model.psi.Lam(end,:)]\[(1./qgaussnew);psib];
    xnew(model.te.xind)  = model.te.LamTegauss\tegaussnew;

end
varargout{1} = xnew;
return