%% function varargout = prad_adas(varargin)
%  Computes radiated power from impurity Z:
%  Pz (W/m^3) = ne (m^-3) nz (m^-3) Lz (W m^3)
%  with Lz the cooling rate taken from ADAS

function varargout = prad_adas(varargin) %#codegen
coder.extrinsic('check_gradients');

if nargin == 0,
    module_config = struct(...
    'modeltype' , 'adas'  ... % 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
    );
    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};
    trap = varargin{2};
    model = varargin{3};
    pr = varargin{4};    
			    
else
    error('must call with 0 or 4 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;
    n1 = stap.n1;
    n2 = stap.n2;
    n3 = stap.n3;
  
    nztot = [n1,n2,n3];

    tabmat = Lz_zave_HDT0_default;

    fn = fieldnames(tabmat);

    for k=1:numel(fn)
       if(tabmat.(fn{k}).Z == model.atom.Z1)
          ind1 = k;
          elseif(tabmat.(fn{k}).Z == model.atom.Z2)
             ind2 = k;
          elseif(tabmat.(fn{k}).Z == model.atom.Z3)
             ind3 = k;
       end
    end

    assert(~isnan(ind1) && ~isnan(ind2) && ~isnan(ind3),...
    'No ADAS data for one or several impurities')
	    
    teimp=tabmat.(fn{ind1}).data(:,1); % Same Te grid for all species
    lz1=tabmat.(fn{ind1}).data(:,2);
    lz2=tabmat.(fn{ind2}).data(:,2);
    lz3=tabmat.(fn{ind3}).data(:,2);

% Interpolation on RAPTOR Te
    rz1=10 .^ pchip(log10(teimp),log10(lz1),log10(te/1e3));
    rz2=10 .^ pchip(log10(teimp),log10(lz2),log10(te/1e3));
    rz3=10 .^ pchip(log10(teimp),log10(lz3),log10(te/1e3));
    rztot = [rz1 rz2 rz3];

    Prad    = ne .* sum(nztot.*rztot,2);
    drz_dte = gradient(rztot')./gradient(repmat(te,1,size(rztot,2))');
    drz_dte = drz_dte';

    if nargout>1
        dPrad_dx = zeros(model.rgrid.nrhogauss,model.dims.nx); % init
        
        if derivative_is_required('te',model)
            dPrad_dTe  = ne.*sum(nztot.*drz_dte,2);
            dPrad_dx = dPrad_dx + bsxfun(@times,dPrad_dTe,stap.dte_dx);
        end
        
        if derivative_is_required('ne',model)
            dPrad_dne  = sum(nztot.*rztot,2);
            dPrad_dx = dPrad_dx + bsxfun(@times,dPrad_dne,stap.dne_dx);
        end

    end

end

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