% ----------------------------------------------------------------------- %
% This function solves the a_0 + a_1f_He + a_2f_He^2 + a_3f_He^3=0, with
% the given temperature [keV], rho = tau_alpha/tau_energy, impurity
% concentration fz, and impurity charge state Z. Use the radiationflag to
% investigate the effect of the radiation:
% radiationflag = 0 - radiation not included
% radiationflag = 1 - radiation included
% The function returns the f_He values in the vector out.
% The used DT cross section is returned in the variable dt_val.
% ----------------------------------------------------------------------- %

function [out,dt_val]= solve_fhe(T,rho,fz,Z,radiationflag)
% The DT fusion cross-section given in table (temperature, <sigma v>)
DT_cross_sect_table = [1, 5.5e-27; 2, 2.6e-25; 5, 1.3e-23; 10, 1.1e-22;...
    20, 4.2e-22; 50, 8.7e-22; 100, 8.5e-22; 200, 6.3e-22; 500, 3.7e-22;...
    1000, 2.7e-22];

% The average radiated power function values for carbon (Z<7) and tungsten
% (Z > 6) in temperatures 1 - 100 keV [page 209, J. Wesson, Tokamaks, 
% 2nd edition].
if(Z<7)
    R_rad_Z_val = 1e-34;
elseif(Z>6)
    R_rad_Z_val = 1e-31;
end
% Find the correct location in the table

dt_val = exp(spline(log(DT_cross_sect_table(:,1)),log(DT_cross_sect_table(:,2)),log(T)));

% The coefficients are given here
% In order to investigate the impact of the different terms, take a look at
% the paper [D. Reiter, NF 1990], and modify the a_x coefficients
% accordingly.
% All the units are in [keV]
c_bremss = 1.71*1e-38;
C = 1.6022e-19;
E_alpha = 3.5*1e3; % in keV
R_bremss_fuel = 1e-3*C^(-1)*c_bremss*sqrt(1000*T);
R_bremss_helium = 4*1e-3*C^(-1)*c_bremss*sqrt(1000*T);
R_rad_Z = 1e-3*C^(-1)*R_rad_Z_val;

a_0 = -(3/2)*T*(fz^3*(-Z^3+Z^2)+fz^2*(4*Z^2-2*Z)...
    +fz*(-5*Z+1)+2);
a_1 = -(3/2)*T*(fz^2*(4*Z-5*Z^2)+fz*(14*Z-4)-9) ...
    + (E_alpha/rho)*(fz^2*Z^2-2*fz*Z+1) ...
    + radiationflag*(4/(dt_val*rho))*((fz*Z-1)*R_bremss_fuel - fz*R_rad_Z);
a_2 =  -(3/2)*T*(fz*(-8*Z+4)+12) ...
    +(E_alpha/rho)*(fz*4*Z-4) ...
    + radiationflag*(4/(dt_val*rho))*(2*R_bremss_fuel-R_bremss_helium);
a_3 = 6*T+4*E_alpha/rho;


x_out = roots([a_3 a_2 a_1 a_0]);
x_out(imag(x_out) ~=0) = [];
out = x_out;