%% Bremsstrahlung Power Loss
%   This module calculates, in a simplistic way, the radiation losses from
%   bremsstrahlung which is caused by electron acceleration by collissions
%   with the plasma ions. The expression for the power loss due to this
%   radiation was taken from Tokamaks (3rd edition) by Wesson, J. The
%   original expression was:
%   $$ 5.35*10^{-37} Z^{2}n_{e}n_{z}T_{e}^{1/2} [W/m^{-3}].$$
%   Here this expression is used in a different from by using Te
%   in eV instead of keV.
%

%% Default Parameters
function [varargout] = pbrem(varargin) %#codegen
coder.extrinsic('check_params');

%   params = P_brem                                    Gets default parameters
%   [Pbrem] = P_brem(stap,trap,model,module_params)
%   [Pbrem, dPbrem_dx] = P_brem(stap,trap,model,module_params)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Input processing
if nargin==0
  varargout{1} = 'src'; % type
  varargout{2} = [];
  return
elseif nargin == 2
  module_params = struct( ...
    'active', false ,   ... % Module is not-active as default
    'check', false      ... % flag for checking derivatives
    );
  
  mm.name = mfilename;
  
  varargout{1} = mm;
  varargout{2} = module_params;
  return %empty call, probably to get default structures
elseif nargin==5 % change this depending on number of inputs to module
    stap = varargin{1};
    model = varargin{4};
    pb = varargin{5}; % distribute inputs
else
    error('must call with 0 or 5 inputs');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize
Pbrem       = zeros(model.rgrid.nrhogauss,1);
dPbrem_dx   = zeros(model.rgrid.nrhogauss,model.dims.nx);

if pb.active == 1
    
    te = stap.te;
    ne = stap.ne;
    ze = stap.ze;
   
    te = abs(te); ne = abs(ne); % avoid < 0 at all costs 
    
    Pbrem   = 1.6918e-38.*ze.*ne.^2.*te.^(1/2); % prefactor from Wesson but now with Te in eV
    
    %% Derivatives
    if nargout >= 2
        dPbrem_dx = zeros(model.rgrid.nrhogauss,model.dims.nx); % init
        
        if derivative_is_required('te',model)
            dPbrem_dTe  = Pbrem./2./te;
            dPbrem_dx = dPbrem_dx + bsxfun(@times,dPbrem_dTe,stap.dte_dx);
        end
        
        if derivative_is_required('ne',model)
            dPbrem_dne  = 2.*Pbrem./ne;
            dPbrem_dx = dPbrem_dx + bsxfun(@times,dPbrem_dne,stap.dne_dx);
        end
    
        if derivative_is_required('ze',model)
            dPbrem_dze  = Pbrem./ze;
            dPbrem_dx = dPbrem_dx + bsxfun(@times,dPbrem_dze,stap.dze_dx);
        end
        
        
    end
end

varargout{1} = Pbrem;
varargout{2} = zeros(size(Pbrem));

if nargout>=3
    varargout{3} = dPbrem_dx;
    varargout{4} = zeros(size(dPbrem_dx));
end

return
