%
% Unit fuel cell code, using the formulation developed in the UPC 
%    workshop lectures
% 
% This code is exercise #4 in those lectures
%
%------------------------

% While in code development, just in case 
clear param sol

% parameters

% physical parameters
param.R = 8.314; 
param.F = 96485;

% Electrochemistry fitted constants
param.alpha = 1;
param.E = 1.28;
param.Cref = 40.9;
param.i0 = 1e-7;
param.delta = 8;
param.Rm = 0.1;

% Operating conditions
param.T = 350;
param.iave = 1;
param.s = 1.1;
param.Ctot = 100;

% Numerical parameters
N= 100;
param.N = N;
param.h = 1/N;
Ntheta = 5;

%-------------------------

% Find theta=0 solution
% First oxygen channel concentration
Ctot = param.Ctot;
C0 = 0.21 * Ctot;
param.C0 = C0;
sol.C = zeros(N+1,1) + C0;

% Now cell voltage
F = param.F;
delta = param.delta;
E = param.E; 
R = param.R;
T = param.T;
alpha = param.alpha;
i0 = param.i0;
Cref = param.Cref;
Rm = param.Rm;
iave = param.iave;
I = iave;
U = E - R*T/alpha/F*(log(I/i0)- log((C0-delta*I)/Cref))-Rm*I;
sol.U = U;

% Now constant nitrogen flux 
s = param.s;
Qn = 0.79/0.21*s*iave/4/F;
param.Qn = Qn;

%--------------------------

% Now do continuation to find the solution 
% A more sophisticated code would use variable steps for the 
%     continuation process
vtheta = linspace(0,1,Ntheta);

for icont = 2:Ntheta
    theta = vtheta(icont);
    final =0;
    if icont == Ntheta
        final = 1;
    end
    [fail sol] = solve_unit(sol, param, theta, final);
    if fail == 1
        disp(sprintf('Solution Failure: trying increasing Ntheta'))
        disp(sprintf('Paused - hit control c'))
        pause
    end
end

h = param.h;
xc = ((1:(N+1))-1)*h;
xi = ((1:N)-0.5)*h;

figure(1) 
plot(xc, sol.C, 'k-')
axis([0 1 0 C0+1])
xlabel('x')
ylabel('Channel O2 concentration')

figure(2) 
plot(xi, sol.i, 'k-')
imax = C0/delta;
axis([0 1 0 imax])
xlabel('x')
ylabel('Current density')

disp(sprintf('Unit cell voltage: %d',sol.U));