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

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

% 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 = 2.0;
param.Ctot = 100;

% Stack parameters
M = 7; % There are 2M-1 cells in the stack
param.M = M; 
param.sanom = 1.2; % anomalous stoich of centre cell #1
param.lambda = 100; % bipolar plate resistivity

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

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

% Find theta=0 solution for the background unit cell 
% 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;

% Convenient place to here, although will only be used in the stack
%    computations later 

sanom = param.sanom;
Qna = 0.79/0.21*sanom*iave/4/F;
param.Qna = Qna;

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

% Now do continuation to find the base unit cell solution 

vtheta = linspace(0,1,Ntheta);

disp(sprintf('Begin unit cell computation at base conditions'))
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('Unit Cell 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')

fprintf('Unit cell voltage: %d \n \n',sol.U);

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

% Now initialize the stack computation

Cs = zeros((N+1),M);
for m=1:M
    Cs(:,m) = sol.C;
end
Us = zeros(N,M)+ sol.U;
sols.Cs = Cs;
sols.Us = Us;

% Use continuation to solve the stack problem

disp(sprintf('Begin stack computation'))

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

Us = sols.Us;
is = sols.is;

figure(3)
plot(xi,Us(:,1),'k-',xi,Us(:,2),'k-',xi,Us(:,3),'k-',xi,Us(:,M),'k-')

figure(4)
plot(xi,is(:,1),'k-',xi,is(:,2),'k-',xi,is(:,3),'k-',xi,is(:,M),'k-')