function varargout = prad_line(varargin) %#codegen
coder.extrinsic('check_gradients');
%% Line Radiation Power Loss
% This module is written to include the Power loss due to line radiation in
% the tokamak plasma. At this stage this Power loss is very simplistically
% realised in RAPTOR by taking of a certain percentage of the total power
% supply:
% Prad = Constant * Te0;
% Plus a contribution from radiation peaks
%
% Ptot = Ptot - Prad;
%
%   params = prad_line        % Gets default parameters
%   [Prad] = prad_line(Poh,Paux,...,...,...,...,module_params)
%
% Radiation peaks are quite simplystically calculated using:
%     factAt1 * ne(1) * exp(-(pr.TeAt1 - te).^2./(pr.WiAt1.^2));

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

if nargin == 0,

    module_config = struct(...
    'modeltype' , 'linerad'  ... % default module
    ); % load default params
  
    varargout{1} = 'prad'; % type
    varargout{2} = module_config; 						 
    return %empty call, probably to get default structures
elseif nargin==2
						 
    module_config = varargin{1};
    module_params = struct( ...
    'check', false, ...     % Check for derivatives
    'active', true,   ... % Module is active as default
    'fact', 1250,      ... % Prefactor Pbulk
    'factAt1', 2e-16, ...
    'TeAt1', 250, ... % temperature for first radiation peak
    'WiAt1', 200, ... % width (in Te space) of first radiation peak
    'factAt2', 2e-17, ... % amount of radiation (Prad = fact*ne*exp(Te-T
    'TeAt2', 2750, ...
    'WiAt2', 2000 ...
    );
    mm.name = mfilename;
    mm.modeltype = module_config.modeltype;
  
    varargout{1} = mm;
    varargout{2} = module_params; % return appropriate parameter structure
    return
elseif nargin==4;       % change this depending on number of inputs to module
    stap = varargin{1};
    %it = varargin{2};
    model = varargin{3};
    pr = varargin{4}; % distribute inputs
else
    error('must call with 0 or 5 inputs');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize
nrg = numel(model.rgrid.rhogauss);

Prad        = zeros(nrg,1);
dPrad_dx   =  zeros(nrg,model.dims.nx);

if pr.active ~= 0
    te  = stap.te;
    ne  = stap.ne;
    te0 = te(1);

    % avoid negative sqrt(imag numbers)
    te = abs(te); ne = abs(ne); te0 = abs(te0);

    Pbulk   = sqrt(1e-3*te0/11)*pr.fact*ones(size(Prad));
    PradAt1 = pr.factAt1*ne(1)*exp(-(pr.TeAt1 - te).^2./(pr.WiAt1.^2));
    PradAt2 = pr.factAt2*ne(1)*exp(-(pr.TeAt2 - te).^2./(pr.WiAt2.^2));
    Prad    = Pbulk + PradAt1 + PradAt2;

    if nargout>1

        dPradAt1_dTe = 2/(pr.WiAt1^2)*(pr.TeAt1-te).*PradAt1;
        dPradAt2_dTe = 2/(pr.WiAt2^2)*(pr.TeAt2-te).*PradAt2;

        dPradAt_dTehat = bsxfun(@times,dPradAt1_dTe + dPradAt2_dTe,model.te.Lamgauss);
        dPbulk_dTehat  = (Pbulk./(2*te0))*model.te.Lamgauss(1,:);

        dPrad_dTehat   = dPradAt_dTehat + dPbulk_dTehat;

        dPrad_dx(:,model.te.xind) = dPrad_dTehat;

        %if pr.check==1
         %   check_gradients(Prad,dPrad_dx,x,g,v,model,pr)
        %end
    end

end
varargout{1} = Prad;
if nargout>1
    varargout{2} = dPrad_dx;
end

end
