clc;
close all;
clear all;

%Fuel cell reaction H2 + 0.5*O2 => H2O Gibbs Free energy [Source FC Explained fig2.4]

%Liquid H2O
T_tab_liq = [25 80];
T_ql = 25:5:80;
DeltaG_tabl = [-247.2 -228.2]; %kJ/mol
DeltaGl = interp1(T_tab_liq,DeltaG_tabl,T_ql); %kJ/mol
%Gas H2O
T_tabg = [100 200 400 600 800 1000];
T_qg = 100:10:1000;
DeltaG_tabg = [-225.2 -220.4 -210.3 -199.6 -188.6 -177.4]; %kJ/mol
DeltaGg = interp1(T_tabg,DeltaG_tabg,T_qg); %kJ/mol

DeltaG0 = interp1(T_tab_liq,DeltaG_tabl,70)*1000; %[J/mol of H2] standard Gibbs free energy of reaction with liquid water as bioproduct at 70C
Tcell = 70+273.15; %[k] Cell temperature
n = 2; %[-] Number of electrons involved in the electrochemical reaction
F = 96485; %[C] Charge of 1 mole of electrones
R = 8.314; %[J/mol] Universal gas constant

P_an = 101325; %[Pa] Anode pressure
P_cat = 101325; %[Pa] Cathode press

P_sat = exp(23.1963-3816.44./(Tcell-46.13)); % [Pa] Saturation pressure of water vapor
p_H2 = (P_an-P_sat)/P_an; %[-] H2 partial pressure with RH100%
p_O2 = (P_cat-P_sat)*0.21/P_cat; %[-] O2 partial pressure
p_H2O = 1; %Water partial pressure assuming liquid water


j_cell = 0.01:0.01:1;  %[A/cm^2] Cell current density
%Reversible voltage [V]
E = -DeltaG0/(n*F);
%Open circuit voltage (OCV)
p_RATIO = (p_H2*p_O2^0.5)/p_H2O;
Vocv = E + (R*Tcell)/(n*F)*log(p_RATIO);  %[V] open circuit voltage
%Activation overvoltage
alpha = 0.2;  %[-] Simmetry factor
j_0 = 0.01;    %[-] 
S_Tafel = (R*Tcell)/(alpha*n*F);%Tafel slope
Vact = S_Tafel*log((j_cell)/j_0);%[V]
%Ohmic overvoltage
ASR = 0.1; %[Ohm*cm^2]
Vohm = j_cell*ASR;%[V]
%Concentration overvoltage
c = 0.032;
j_cell_L = 1;%[A/cm^2]
Vcon = c*log(j_cell_L./(j_cell_L-j_cell));%[V]
%Cell voltage
V = Vocv-Vact-Vohm-Vcon;%[V]
A_cell = 10*10;%cm2
ncell = 50;
P_cellcm2 = j_cell.*V;
I_cell = j_cell*A_cell;
V_stack = A_cell*ncell*V;
P_cell = A_cell*j_cell.*V;
P_stack = P_cell*ncell;

%Plotting results
figure(1)
plot(j_cell,V,j_cell,Vact,j_cell,Vohm,j_cell,Vcon,j_cell,P_cellcm2),legend('Vcell','Vact','Vohm','Vcon','P_cell')
xlabel('Current Density [A/cm2]')
ylabel('Voltage [volt]  // Power [watt]')
figure(2)
plot(T_ql,DeltaGl,T_qg,DeltaGg),legend('DeltaG(g)','DeltaG(l)')
xlabel('Temperature [degreeC]')
ylabel('DeltaG')
figure(3)
plot(I_cell,V_stack,I_cell,P_stack),legend('V_{stack}','P_{stack}')
xlabel('Current [A]')
ylabel('Voltage [volt] (Blue) // Power [watt] (red)')