%This script solves the parameter space in which ignition is achievable for given values of
%rho, which is the ratio of helium energy confinement time to general plasma confinement time
%fz, which is the fraction of impurities
%Z, which is the impurity charge state.
%Additionally, the script plots the ignition parameter space for pure DT fusion for the sake of comparison.
%Running the script requires solve_fthe.m, which is used for solving the helium fraction.


%INPUT PARAMETERS
%You can play with these parameters
%%%%%%%
rhovec = [5, 5, 5, 5];  % Give the rho values (>0) to be plotted, e.g. rho = 5
fzvec = [0, 4e-5, 8e-5, 1e-4]; %Give the impurity concentration vector, e.g. fz = 1e-4.
Zvec = [50, 50, 50, 50]; % Give the impurity charge vector, e.g. Z=6 for carbon.
radiationflag = 1;% Give the radiation flag (1 for radition included, 0 for radiation neglected).
colvec = [0,0,0; 1,0,0; 0,0,1; 1,0.4,0];%Give the colors of the curves, just cosmetic.


%SCRIPT STARTS DOING ITS THING
%You don't have to change anything in here.
%%%%%%%%%%%
for j=1:length(rhovec) % looping over the rho values
    rho=rhovec(j);
    fz = fzvec(j);
    Z = Zvec(j);
    for i=4.5:0.1:500; % looping over the temperatures from 4.5 to 100 keV
       % Solving the f_He values
       [x_out,dtc] = solve_fhe(i,rho,fz,Z,radiationflag);
       % Excluding values outside 0<f_he<0.5
       inds1 = x_out<0.51;
       x_out=x_out(inds1);
       inds = length(x_out);
       % Plotting the obtained solutions (2 points)
       for k=1:inds
           y = 4*x_out(k)/(rho*(1-2*x_out(k))^2*dtc);
           loglog(i,y,'s','LineWidth',2,'color',colvec(j,:));
           hold on
       end
    end
end

%PLOTTING STUFF
%%%%%%%%%
disp('Plotting intensifies...')
% Setting up the figure
set(gca,'FontSize',14);
axis square;
axis([4 1000 1e20 1e23]);
%set(gca,'FontWeight','bold');
%set(gca,'FontSize',18);
set(gca,'TickLength',[0.025 0.025]);
xlabel('T [keV]');
ylabel('n\tau [m^{-3}s]');
set(gca,'XTick',[1e0 1e1 1e2 1e3])

% Plot the pure DT values
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];
x = DT_cross_sect_table(:,1);
dtc = DT_cross_sect_table(:,2);

ya = 12*x*1e3*1.6022e-19./(dtc*3.5*1e6*1.6022e-19 - 4*1.71*1e-38*sqrt(x*1e3));
ind = ya > 0;
loglog(x(ind),ya(ind),'m^','LineWidth',2)
grid on
hold off