% nonisothermal_CSTR_Da_scan.m
% This MATLAB program performs arc-length
% continuation to plot phi_A and theta
% vs. Da at fixed values of the other parameters.
% Actually, the system steps in log space of the Damkoehler
% number.
% K.J. Beers. MIT ChE. 9/24/03
function iflag_main = nonisothermal_CSTR_Da_scan();
iflag_main = 0;

% First, input the system parameters
disp('nonisothermal_CSTR_Da_scan.m:');
disp('Enter parameters ...');
disp(' ');
disp('Da is Damkoehler #, inportance of rxn to convection.');
Da_0 = input('Enter lowest Da number : '); Da_0_log = log10(Da_0);
Da_1 = input('Enter highest Da number : '); Da_1_log = log10(Da_1);
disp(' '); disp('beta is dimensionless heat of reaction');
beta = input('Enter fixed beta value : ');
disp(' '); disp('chi is dimensionless cooling heat transfer strength');
chi = input('Enter fixed chi value : ');
disp(' '); disp('gamma is dimensionless activation energy');
gamma = input('Enter fixed gamma value : ');
disp(' '); disp('theta_c is dimensionless temp of coolant');
theta_c = input('Enter fixed theta_c value : ');

ArcLenParam.atol = 1e-6;
    
% generate linear path in parameter space by changing
% Damkoehler number logarithmically
Param_0 = [Da_0_log; beta; chi; gamma; theta_c];
Param_1 = [Da_1_log; beta; chi; gamma; theta_c];
    
% set initial guess for solution at Param_0
x0 = [1; 1];
    
% call arc-length continuation routine to generate curve
fun_name = 'nonisothermal_CSTR_calc_f';
[x,Param,lambda,f_norm,stab] = ...
    arclength_continuation(fun_name,Param_0,Param_1,x0,ArcLenParam);

% extract Da numbers along graph
Da = 10.^Param(1,:);

% Next, identify where stability changes along curve
% record values for first point
istab_point(1) = 1;
count = 1;
for k=2:length(stab)
    if(stab(k) ~= stab(k-1))  % if change in stability
        count = count + 1;
        istab_point(count) = k;
    end
end
% record last point
istab_point(count+1) = length(stab);
% Now, go along curves and plot each section, using
% a solid line for stable sections and a dashed line
% for unstable section
figure;
for isection=1:count
    istart = istab_point(isection);
    iend = istab_point(isection+1);
    if(stab(istart))  % if stable
        subplot(2,1,1);
        semilogx(Da(istart:iend),x(1,istart:iend),'-');
        hold on;
        subplot(2,1,2);
        semilogx(Da(istart:iend),x(2,istart:iend),'-');
        hold on;
    else  % if unstable
        subplot(2,1,1);
        semilogx(Da(istart:iend),x(1,istart:iend),'-.');
        hold on;
        subplot(2,1,2);
        semilogx(Da(istart:iend),x(2,istart:iend),'-.');
        hold on;
    end
    % for intermediate points, demark points at which
    % stability changes by a circle
    if(istart ~= 1)
        subplot(2,1,1);
        semilogx(Da(istart),x(1,istart),'o');
        subplot(2,1,2);
        semilogx(Da(istart),x(2,istart),'o');
    end
end
% add labels and title
subplot(2,1,1);
xlabel('Da');  ylabel('\phi_A');
title_phrase = ['\phi_A, \theta vs. Da: ', ...
        '\beta = ', num2str(beta), ...
        ', \chi = ', num2str(chi), ...
        ', \gamma = ', num2str(gamma), ...
        ', \theta_c = ', num2str(theta_c)];
title(title_phrase);
subplot(2,1,2)
xlabel('Da');  ylabel('\theta');

% now make separate plot of infinity norm of f vs. lambda as
% check that the arclength continuation worked correctly.
figure;  semilogy(lambda,f_norm);
xlabel('lambda');  ylabel('||f||_{infinity}');
title('Infinity norm of f along curve');
axis tight;

iflag_main = 1;
return;