function dT_eVdt = thermal_energy_balance(t,T,time,P_in,p)

q_e = 1.60217657E-19; % electron charge [coulomb] or Boltzmann constant [J/eV]
P_aux = interp1(time,P_in,t); % Interpolate the input signal (time,u) at time t
S = sources(abs(T),P_aux,p);

dT_eVdt = 1/(3*p.n*q_e)*(S.S_alpha + S.S_ohm + S.S_aux - S.S_rad - S.S_cond);   % no factor 2 because T=Te=Ti and n=ne=ni are tacitly assumed such that ptot=ne Te+ni Ti=2nT
end