function S = sources(T,P_aux,p)

%%
a = [-21.38, -25.20, -7.101e-2 1.938e-4 4.925e-6 -3.984e-8];
alp = 0.2935;
E_alpha = 4.5E6; % [eV]
q_e = 1.60217657E-19; % electron charge [coulomb] or Boltzmann constant [J/eV]

V = (2*pi*p.R_0)*(pi*p.kappa*p.a^2); % plasma volume [m^3]
epsilon = p.a/p.R_0;

sigmav = 1e-6*exp(a(1)./((T/1000).^alp) + a(2) + a(3)*(T/1000) + ...
         a(4)*(T/1000).^2 + a(5)*(T/1000).^3 + a(6).*(T/1000).^4);

P_alpha = V * p.f_DT/(1+p.f_DT)^2*q_e*E_alpha*p.n^2*sigmav;
P_ohm = (5.6E4/(1-1.31*epsilon^0.5+0.46*epsilon))*...
    ((p.R_0*p.I_p.^2)./(p.a^2*p.kappa*(T/1000).^(3/2)));

tau_e = 0.145*p.I_p^(0.93)*...
              p.R_0^(1.39)*...
              p.a^(0.58)*...
              p.kappa^(0.78)*...
              (p.n/1E20)^(0.41)*...
              p.B_0^(0.15)*...
              p.A^(0.19)*...
              ((P_aux+P_alpha+P_ohm)/1E6).^(-0.69);

pressure = 2*p.n*q_e*T;
          
%%
S.S_aux = P_aux/V;
S.S_ohm = P_ohm/V;
S.S_alpha = P_alpha/V;
S.S_rad = 5.35E3*p.Z_eff*(p.n/1E20)^2*sqrt(T/1000);
S.S_cond = 3/2*pressure./tau_e;
S.tau_e = tau_e;
S.pressure = pressure;

end