clear;
close all;
clc;

% Task 2

% Number of nodes in the grid
N = 100;


% Collecting the grid used in th previous exercise (the nodes will be along
% the z-axis)
[z_hat, w, A, B] = FiniteDifferenceGrid (N);

% The dimensionless parameters
alpha = 10;
beta = 1000;

% Need initial guesses for fsolve, choosing values in the middle of range:
xA0 = 0.5*ones(N,1);

% Calling fsolve to find xA:
xA = fsolve(@(xA) sysEquations(xA, A, B, alpha, beta, N), xA0);

% Calculating xB from the fact that we have a binary mixture
xB = 1-xA;


%% Plotting the result
% Properties of the plot
figProps = struct('Color', [1 1 1], 'OuterPosition', [170, 170, 1000, 700]);
fontProps = struct('FontName', 'Calibri', 'FontSize', 18, 'FontWeight', 'bold');
% Initializing the plot
fig = figure(1);
ax = axes;
% Plotting the results
plot(z_hat, xA, 'Color', 'r', 'LineWidth', 2, 'Marker', 'o', 'MarkerEdgeColor', 'r', 'MarkerFaceColor', 'none', 'DisplayName', 'x_A');
hold on;
plot(z_hat, xB, 'Color', 'b', 'LineWidth', 2, 'Marker', 'o', 'MarkerEdgeColor', 'b', 'MarkerFaceColor', 'none', 'DisplayName', 'x_B');
% Tidying up
set(fig, figProps);
set(ax, 'FontSize', fontProps.FontSize);
set(xlabel("z' = z/L", 'interpreter', 'tex'), fontProps);
set(ylabel('Molar fractions, x_i'), fontProps);
set(title(['Mole fractions for \alpha = ', num2str(alpha), ' and \beta = ', num2str(beta)],'interpreter', 'tex'), fontProps);
xlim([0 1]);
box("on");
grid("on");
set(legend(), 'FontName', fontProps.FontName)
%Saving the figure
saveas(fig, 'Ex6_2_plot', 'jpg')

%% Functions used in the problem

% Creating a function for the f(xA) defined on my submission
function toReturn = f(x)
   toReturn = log(1/(1-0.5*x)); 
end

% Creating the equation system to be solved
function wantToBeZero = sysEquations(xA, dxA_matrix, ...
    d2xA2_matrix, param1, param2, numberOfNodes)
    
    %Creating empty array to optimize performance
    wantToBeZero = zeros(numberOfNodes, 1);

    % Computing the differentials values beforehand using the grid
    d2xA_dz2 = d2xA2_matrix*xA;
    dxA_dz = dxA_matrix*xA;

    %% Boundary conditions
    % At z=0
    wantToBeZero(1) = xA(1)-1;

    %% General solution for the middle nodes:
    for i = 2:(numberOfNodes-1)
        wantToBeZero(i) = d2xA_dz2(i) - param1*dxA_dz(i) - param2*f(xA(i));
    end
    % At z=1
    wantToBeZero(numberOfNodes) = dxA_dz(numberOfNodes);
end 