function [varargout] = neHmodeEvo(varargin)
% Parametrization for density profile evolution during L-H transition
module_params = struct( ...
    'LHtime', 60, ... % time of LH transition
    'LHduration',10, ...    % duration of LH transition
    'HLtime', 60, ... % time of LH transition
    'HLduration',10, ...    % duration of LH transition
    'peaking',1.3,...  % peaking from pedestal to center (flat=1)
    'neL_t',1e19,...    % constant (or time evolution of) density during L mode
    'neH_t',10e19...    % constant (or time evolution of) density during H mode
    );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Input processing

if nargin == 0,
    varargout{1} = module_params;
    return %empty call, probably to get default structures
elseif nargin==3;       % change this depending on number of inputs to module
    [model,params,neparams] = deal(varargin{:}); % distribute inputs
    
else
    error('must call with 0 or 3 inputs');
end

%% density  profile shapes
% ne as in H-mode, normalized to ne(0) = 1
irhoped = find(model.rgrid.rhogauss>=.92,1,'first');

neH_prof = zeros(numel(model.rgrid.rhogauss),1); % init
neH_prof(1:irhoped) = (neparams.peaking-1)*(1-linspace(0,1,irhoped).^2) + 1.0;
neH_prof(irhoped:end) = linspace(1,.4,numel(neH_prof(irhoped:end)));

%plot(ne_hmode_0)

% ne in Lmode as fitted from CRONOS, normalized to ne(0) = 1
neL_prof = 1 - 0.4*model.rgrid.rhogauss.^2;
% neL_prof = neH_prof;

%% L and H mode density evolutions in time
neL_t = neparams.neL_t.*ones(size(params.tgrid)); % time trace of ne evolution for L mode
neH_t = neparams.neH_t.*ones(size(params.tgrid)); % time trace of ne evolution for H mode

%ne_lmode = rampfun(params.tgrid,20,1.1e19,100,3e19);
%ne0(indLHtrans:end)  = ne0(indLHtrans)*rampfun(params.tgrid(indLHtrans:end),neparams.LHtrans,1,neparams.LHstop,.9e20/ne0(indLHtrans));

%% L-H transition
% linear interpolation with ne = lambda*ne_Hmode + (1-lambda)*ne_Lmode
rampfun = @(t,tmin,ymin,tmax,ymax) max(ymin,min((ymax-ymin)/(tmax-tmin)*(t-tmin),ymax-ymin)+ymin); % anonymous function for ramps
lambda = rampfun(params.tgrid,neparams.LHtime,0,neparams.LHtime+neparams.LHduration,1) - ...
         rampfun(params.tgrid,neparams.HLtime,0,neparams.HLtime+neparams.HLduration,1);

ne_LH = bsxfun(@times,lambda,neH_prof*neH_t) + bsxfun(@times,(1-lambda),neL_prof*neL_t);

%% output
varargout{1} = ne_LH;