clc; clear; close all; load gibbseng.mat;

%% Task 3 (B)
% This script attempts to solve n5 using the built in Matlab symbolic
% toolbox function "solve". 

% IMPORTANT:
% This file requires the matlab "gibbseng.mat" file and "GibbsEngCalcs"
% function to work. 

%% Variables

carbon=8;
hydrogen=18;
oxygen=0;
molar_air_stoich=12.5;
temperature=(100)';
pressure=(20)';
fuelequiv=(0.8)';
Error = 1e-12;


%% Symbolic toolbox: solve
% Notes: This system seems to work well at solving one set of temperature, 
% pressure and FER variables at a time. When a for loop was integrated, 
% the solver took far to long to converge for all temperature, pressure
% and FER values, so was removed. 

[kco2,kh2o] = GibbsEngCalcs(temperature,gibbseng);

syms n1 n2 n3 n4 n5 n6 a ntsum n5ntvalue

eqn1 = n3==(carbon.*fuelequiv)./(kco2.*a+1);
eqn2 = n4==((hydrogen./2).*fuelequiv)./(kh2o.*a+1);
eqn3 = n1==carbon.*fuelequiv-n3;
eqn4 = n2==(hydrogen./2).*fuelequiv-n4;
eqn5 = n5==(molar_air_stoich.*2+oxygen.*fuelequiv-2.*n1-n2-n3)./2;
eqn6 = n6==molar_air_stoich.*3.76;
eqn7 = ntsum==n1+n2+n3+n4+n5+n6;
eqn8 = a==sqrt(n5ntvalue.*pressure);
eqn9 = Error==n5ntvalue-n5./ntsum;

awns = solve([eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8,eqn9],[n1,n2,n3,n4,n5,n6,a,ntsum,n5ntvalue]);

n1sol = double(awns.n1(4));
n2sol = double(awns.n2(4));
n3sol = double(awns.n3(4));
n4sol = double(awns.n4(4));
n5sol = double(awns.n5(4));
n6sol = double(awns.n6(4));
n5ntvaluesol = double(awns.n5ntvalue(4));

fprintf('Solver Awnsers: n1 = %.3f, n2 = %.3f, n3 = %.3f, n4 = %.3f, n5 = %.3f, n6 = %.3f, n5nt Value %.4f\n',n1sol,n2sol,n3sol,n4sol,n5sol,n6sol,n5ntvaluesol)