function varargout = palpha(varargin)%#codegen
coder.extrinsic('check_gradients');
% Alpha particle module for RAPTOR
% Syntax:
%   params = P_alpha; % get default parameters
%   [Palphae,Palphai,dPalphae_dx,dPalphae_du,dPalphai_dx,dPalphai_du] = P_alpha(stap,geop,u,model,params) % run module
%
% returns Alpha particle heating and derivatives. All outputs are on grid [rhogauss]
%

%Input processing
if nargin==0
  varargout{1} = 'src'; % type
  varargout{2} = [];
  return
elseif nargin == 2,
  %% DEFAULT PARAMETERS
  default_params = struct(...
    'active',false,     ... % do not use by default
    'calpha',1,         ... % power scaling parameter
    'check',false,      ... % optional check of derivatives
    'DTfraction',.5    ... % nT/(nD+nT) fraction of Tritium vs Deuterium
    );
  
  mm.name = mfilename;
  varargout{1} = mm;
  varargout{2} = default_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};
  pap = varargin{5}; % distribute inputs
else
  error('must call with 0 or 5 inputs');
end

%% Constants
pm      = 1.6725e-27;                                               % Proton mass [kg]
nm      = 1.6748e-27;                                               % Neutron mass [kg]

me      = 9.1091e-31;                                               % Electron mass [kg]
mD      = nm+pm;                                                    % Deuterium mass [kg]
mT      = nm+2*pm;                                                  % Tritium mass [kg]
mi      = (pap.DTfraction*mT + (1-pap.DTfraction)*mD);        % Mean plasma ion mass [kg]
ma      = 2*pm + 2*nm ;                                             % Alpha particle mass [kg]

Afrac   = 3.51e6/((3*pi^(1/2)/4)^(2/3)*(mi/me)^(1/3)*(ma/mi));
Ealpha  = 5.68e-13;                                                 % Alpha particle birth energy [J]

%%%% ALPHA PARTICLE HEATING MODULE

% allocate
Palphai = zeros(model.rgrid.nrhogauss,1);
Palphae = zeros(model.rgrid.nrhogauss,1);
if nargout>2
  dPalphai_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
  dPalphae_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
end

if pap.active
  
  % simple equation ~p^2
  
  te = stap.te;
  dte_dx = stap.dte_dx;
  %ne = eval_ne(x,g,v,model,true);
  
  ti = stap.ti;
  dti_dx = stap.dti_dx;
  ni = stap.ni;
  dni_dx = stap.dni_dx;
  
  nT = pap.DTfraction*ni;
  nD = (1-pap.DTfraction)*ni;
  
  y = Afrac./te;                                          % te should be given in eV
  
  if nargout>2
    [R, dR_dTi]     = reactionrate(ti);                 % ti should be given in eV
    [phi, dphi_dy]  = fastionenergyfraction(y);         % te should be given in eV;
  elseif nargout == 2
    [R] = reactionrate(ti);                             % ti should be given in eV
    [phi] = fastionenergyfraction(y);                   % te should be given in eV;
  end
  
  Bprofile    = pap.calpha*nT.*nD.*R;              % nT,nD should be given in m^3 s^-1
  Palphai     = Ealpha*phi.*Bprofile;
  Palphae     = Ealpha*(1-phi).*Bprofile;
  
  % derivatives
  if nargout>2
    
    % calculate T-derivatives
    dPalphai_dTe        = Ealpha*Bprofile.*dphi_dy.*(-y.^2/Afrac); % chain rule dphi_dy*dy_dTe
    dPalphai_dTi        = Ealpha*phi.*dR_dTi.*pap.calpha.*nT.*nD;
    dPalphae_dTe        = Ealpha*Bprofile.*dphi_dy.*(y.^2/Afrac);
    dPalphae_dTi        = Ealpha*(1-phi).*dR_dTi.*pap.calpha.*nT.*nD;
    % calculate n-derivatives
    dPalphai_dni        = Ealpha*phi.*pap.calpha.*(nT*(1-pap.DTfraction) + pap.DTfraction*nD).*R;
    dPalphae_dni        = Ealpha*(1-phi).*pap.calpha.*(nT*(1-pap.DTfraction) + pap.DTfraction*nD).*R;
    
    if strcmp(model.ti.method, 'state') && strcmp(model.ne.method, 'state')
      dPalphai_dTehat     = bsxfun(@times,dPalphai_dTe,dte_dx(:,model.te.xind));
      dPalphai_dTihat     = bsxfun(@times,dPalphai_dTi,dti_dx(:,model.ti.xind));
      dPalphai_dnehat     = bsxfun(@times,dPalphai_dni,dni_dx(:,model.ne.xind));
      dPalphae_dTehat     = bsxfun(@times,dPalphae_dTe,dte_dx(:,model.te.xind));
      dPalphae_dTihat     = bsxfun(@times,dPalphae_dTi,dti_dx(:,model.ti.xind));
      dPalphae_dnehat     = bsxfun(@times,dPalphae_dni,dni_dx(:,model.ne.xind));
      
      % assign to dPalpha_dx
      dPalphai_dx(:,model.te.xind) = dPalphai_dTehat;
      dPalphai_dx(:,model.ti.xind) = dPalphai_dTihat;
      dPalphai_dx(:,model.ne.xind) = dPalphai_dnehat;
      dPalphae_dx(:,model.te.xind) = dPalphae_dTehat;
      dPalphae_dx(:,model.ti.xind) = dPalphae_dTihat;
      dPalphae_dx(:,model.ne.xind) = dPalphae_dnehat;
      
    elseif strcmp(model.ti.method, 'state')
      dPalphai_dTehat     = bsxfun(@times,dPalphai_dTe,dte_dx(:,model.te.xind));
      dPalphai_dTihat     = bsxfun(@times,dPalphai_dTi,dti_dx(:,model.ti.xind));
      dPalphae_dTehat     = bsxfun(@times,dPalphae_dTe,dte_dx(:,model.te.xind));
      dPalphae_dTihat     = bsxfun(@times,dPalphae_dTi,dti_dx(:,model.ti.xind));
      
      % assign to dPalpha_dx
      dPalphai_dx(:,model.te.xind) = dPalphai_dTehat;
      dPalphai_dx(:,model.ti.xind) = dPalphai_dTihat;
      dPalphae_dx(:,model.te.xind) = dPalphae_dTehat;
      dPalphae_dx(:,model.ti.xind) = dPalphae_dTihat;
    
    else
      dPalphai_dTehat     = bsxfun(@times,dPalphai_dTe,dte_dx(:,model.te.xind)) + ...
        bsxfun(@times,dPalphai_dTi,dti_dx(:,model.te.xind)); % w.r.t. spline coefficients
      dPalphae_dTehat     = bsxfun(@times,dPalphae_dTe,dte_dx(:,model.te.xind)) + ...;
        bsxfun(@times,dPalphae_dTi,dti_dx(:,model.te.xind));
    
      % assign to dPalpha_dx
      dPalphai_dx(:,model.te.xind) = dPalphai_dTehat;
      dPalphae_dx(:,model.te.xind) = dPalphae_dTehat;
    end
%     if pap.check
%        check_gradients(x,g,v,model,pap);
%     end
  end
  
end
% assign to varargout

varargout{1} = Palphae;
varargout{2} = Palphai;
if nargout>2
  varargout{3} = dPalphae_dx;
  varargout{4} = dPalphai_dx;
end

end


%% Coefficients
% The form of the reaction rate was extrapolated in the following
% reference:
% "Improved formulas for fusion cross-sections and thermal
% reactivities" - H.S. Bosch, G.M. Hale, Nucl. Fusion 32 611.
% The form of the reactivity was validated for temperature 0.2 keV < T <
% 100 keV. The reactivity is given in m^3/s when T is given in keV.
function varargout = reactionrate(ti)
%Timin    = 0.2;
%Timax    = 100;

TkeV        = 1e-3*ti;                       % conversion from ev to keV
%ind = (any(TkeV<=Timin | TkeV>=Timax));
%TkeV(ind) = 0;

Bg=     34.3827;
Mrc2=   1124656;
C1=     1.17302e-9;
C2=     1.51361e-2;
C3=     7.51886e-2;
C4=     4.60643e-3;
C5=     1.35e-2;
C6=     -1.06750e-4;
C7=     1.366e-5;

T2 = TkeV.*TkeV;
T3 = T2.*TkeV;
T4 = T3.*TkeV;
T5 = T4.*TkeV;
T6 = T3.*T3;

theta       =   TkeV./(1-(TkeV .* ( C2+TkeV.*(C4+TkeV.*C6)))./(1+TkeV.*(C3 + TkeV.* (C5+TkeV.*C7))));
e           =   (Bg.^2./4./theta).^(1/3);
R           =   C1.*theta.*sqrt(e./(Mrc2.*T3)).*exp(-3.*e);
%R(ind)      =   0;
varargout{1}=   R.*1e-6;                    % Conversion from cm^3/s to m^3/s.

if nargout == 2
  dR_dT = (2*C1.*exp(-3.*e).*sqrt(e./(Mrc2.*T3)).*(-4 + TkeV.*(9*C2-8*C3)+TkeV.^2.*(4*C2*C3-4*C3^2+14*C4-8*C5)+...
    T3.*(9*C3*C4-C2*C5-8*C3*C5+19*C6-8*C7)+T4.*(4*C4*C5-4*C5^2+14*C3*C6-6*C2*C7-8*C3*C7)+...
    T5.*(9*C5*C6-C4*C7-8*C5*C7)+T6.*(4*C6*C7-4*C7^2)+6.*e.*(1+2.*C3.*TkeV+TkeV.^2.*(-C2*C3+C3^2+C4+2*C5)+...
    2.*T3.*(-C2*C5+C3*C5+C6+C7)+T4.*(-C4*C5+C5^2+C3*C6-3*C2*C7+2*C3*C7)+2.*T5.*C7.*(C5-C4)+C7.*T6.*(C7-C6))))...
    ./(12.*(1 + TkeV.*(-C2+C3+TkeV.*(-C4+C5-C6.*TkeV+C7.*TkeV))).^2);
  
  dR_dT = dR_dT.*1e-9;                % Conversion to m^3 s^-1 eV^-1
  %  dR_dT(ind,:) = 0;
  varargout{2} = dR_dT;
  
end
end