function Z = p1model( x, theta )
%
%   Z = p1model( x, theta )
%
% This is the model of the calculation of the impedance
% of a PBI fuel cell, obtained by means of a impedance
% test with variable voltage frequency.
% 
% This is the phase-1 model, assuming a constant value
% for k's and C's.
%
% The parameter vector theta contains ka, kb, Ca, Cb, Ri, 
% iC, Lext. SI units, assuming atmospheric pressure as the
% reference for the k's.
%
% The variable vector includes w, T, V, pk_O2, pk_N2, pk_H2O,
% pa_H2, pa_CO, pa_CO2, pa_H2O. SI units.
%
% The output vector, Z, contains the real and imaginary part
% of the impedance. Complex numbers are formally avoided
% for simplicity.
%
%
% Some constants are embedded in this function as they are 
% assumed not to vary in the experiments.
% 
% Test values:
 theta = [ 1E-12, 1E-3, 500, 10, 2E-4, 5, 1E-8 ]' ;
 x     = [ 100, 430, 0.6, 20000, 80000, 1325, 101325, 0, 0, 0 ]' ;

% $Id: p1model.m,v 1.3 2005/08/16 09:37:05 zenith Exp $

% Electrode parameters: porosity, tortuosity, thickness.
% Can hopefully be estimated by lab responsibles.
eps_a = 0.4 ;     % Anode porosity
eps_c = 0.4 ;     % Cathode porosity
tau_a = 7   ;     % Anode tortuosity
tau_c = 7   ;     % Cathode tortuosity
L_a   = 3E-4;     % Anode diffusion layer thickness, m
L_c   = 3E-4;     % Cathode diffusion layer thickness, m
A     = 4E-4;     % Cell area, m² 

% Electrochemical constants:
n_a    = 2     ;    % Electrons exchanged in anodic reaction
n_c    = 4     ;    % Electrons exchanged in cathodic reaction
alfa_a = 0.5   ;    % Anodic simmetry coefficient
alfa_c = 0.5   ;    % Cathodic simmetry coefficient
F      = 96485 ;    % Faraday constant, C/mol
R      = 8.314 ;    % Perfect-gas universal constant, J/mol/K

% Model parameters from input theta:
k_a   = theta( 1 );   % Anodic reaction constant, [?]
k_c   = theta( 2 );   % Cathodic reaction constant, [?]
C_a   = theta( 5 );   % Anodic capacitance, F/m²
C_c   = theta( 4 );   % Cathodic capacitance, F/m²
ri    = theta( 5 );   % Internal resistance, ohm m²
iC    = theta( 6 );   % Crossover current, A/m²
Lext  = theta( 7 );   % Inductance of the external circuit, H

% Model variables from input x:
w  = x( 1 )    ;    % Pulsation, rad/s
T  = x( 2 )    ;    % Temperature, K
V  = x( 3 )    ;    % Voltage, V
pc = x( 4:6 )  ;    % Cathode pressures (oxygen, nitrogen, water), Pa
pa = x( 7:10 ) ;    % Anode pressures (hydrogen, carbon monoxide, carbon dioxide, water), Pa

% Evaluating the pressure-diffusivity products 
% at the given temperature
pD_O2_H2O = 1.2649E-4 * T^1.75 ;
pD_H2O_N2 = 1.2231E-4 * T^1.75 ;
pD_O2_N2  = 9.7684E-5 * T^1.75 ;
pD_H2_CO  = 3.7444E-4 * T^1.75 ;
pD_H2_CO2 = 3.1551E-4 * T^1.75 ;
pD_H2_H2O = 4.3080E-4 * T^1.75 ;

if ( V > 0 ) & ( V < i2V( 0 ) )
    i = fsolve( @VZERO, 1000 ) ;
else 
    error('Unacceptable value for V, beyond open-circuit value.')
end

ir = i + iC ;

pR_H2  = pa( 1 ) - tau_a^2 / eps_a * L_a * R * T * (  pa( 2 )/pD_H2_CO + pa( 3 )/pD_H2_CO2 + pa( 4 )/pD_H2_H2O  ) * ir / F / n_a ;
pR_O2  = pc( 1 ) - tau_c^2 / eps_c * L_c * R * T * ( ( 2*pc( 1 ) + pc( 3 ) ) / pD_O2_H2O + pc( 2 ) / pD_O2_N2    ) * ir / F / n_c ;
pR_H2O = pc( 3 ) + tau_c^2 / eps_c * L_c * R * T * ( ( 2*pc( 1 ) + pc( 3 ) ) / pD_O2_H2O + 2*pc( 2 ) / pD_H2O_N2 ) * ir / F / n_c ;


function eta = i2eta( ir )

        % Calculating the exchange current densities
        i0_a = k_a * pR_H2^alfa_a * n_a * F ;
        i0_c = k_c * pR_O2^( 0.5 * alfa_c ) * pR_H2O^( 1 - alfa_c ) * n_c * F ;

        % Computing first guesses with the Tafel formula

        if ( ir < i0_a ) 
            eta_tafel( 1 ) = 0 ;
        else
            eta_tafel( 1 ) = R * T / alfa_a / n_a / F * log( ir / i0_a ) ;
        end

        if ( ir < i0_c ) 
            eta_tafel( 2 ) = 0 ;
        else
            eta_tafel( 2 ) = R * T / alfa_c / n_c / F * log( ir / i0_c ) ;
        end
end

eta = i2eta( ir ) ;

R.a = eta( 1 ) / ir / A ;
R.c = eta( 2 ) / ir / A ;
Ri  = ri / A ;

Z( 1 ) = Ri + R.a / ( 1 + w^2 * R.a^2 * C_a^2 ) + R.c / ( 1 + w^2 * R.c^2 * C_c^2 ) ;
Z( 2 ) = w * Lext - w * R.a^2 * C_a / ( 1 + w^2 * R.a^2 * C_a^2 ) - w * R.c^2 * C_c / ( 1 + w^2 * R.c^2 * C_c^2 ) ;

    function E = VZERO( i )
    E = V - i2V( i ) ;
    end


    function V = i2V( i )
    % Calculating the reaction-site concentrations
    % of relevant species:
    ir = i + iC ;

    pR_H2  = pa( 1 ) - tau_a^2 / eps_a * L_a * R * T * (  pa( 2 )/pD_H2_CO + pa( 3 )/pD_H2_CO2 + pa( 4 )/pD_H2_H2O  ) * ir / F / n_a ;
    pR_O2  = pc( 1 ) - tau_c^2 / eps_c * L_c * R * T * ( ( 2*pc( 1 ) + pc( 3 ) ) / pD_O2_H2O + pc( 2 ) / pD_O2_N2    ) * ir / F / n_c ;
    pR_H2O = pc( 3 ) + tau_c^2 / eps_c * L_c * R * T * ( ( 2*pc( 1 ) + pc( 3 ) ) / pD_O2_H2O + 2*pc( 2 ) / pD_H2O_N2 ) * ir / F / n_c ;
    
    eta = i2eta( ir ) ;
    V = E0( T ) - sum(eta) - ri * ir  ;

        function eta = i2eta( ir )

        % Calculating the exchange current densities
        i0_a = k_a * pR_H2^alfa_a * n_a * F ;
        i0_c = k_c * pR_O2^( 0.5 * alfa_c ) * pR_H2O^( 1 - alfa_c ) * n_c * F ;

        % Computing first guesses with the Tafel formula

        if ( ir < i0_a ) 
            eta_tafel( 1 ) = 0 ;
        else
            eta_tafel( 1 ) = R * T / alfa_a / n_a / F * log( ir / i0_a ) ;
        end

        if ( ir < i0_c ) 
            eta_tafel( 2 ) = 0 ;
        else
            eta_tafel( 2 ) = R * T / alfa_c / n_c / F * log( ir / i0_c ) ;
        end

        eta = fsolve( @etaloop, eta_tafel ) ;

            function error = etaloop( eta )
            
            error( 1 ) = ir - butlerVolmer( eta(1), i0_a, alfa_a, n_a ) ;
            error( 2 ) = ir - butlerVolmer( eta(2), i0_c, alfa_c, n_c ) ;
            
                function ir = butlerVolmer( eta, i0, alfa, n )
                % Calculates the reaction current based on the catalytic overvoltage
                ir = i0 * ( exp( alfa * n * F / R / T * eta ) - exp( - ( 1 - alfa ) * n * F / R / T * eta ) ) ;
                end
            end
        end
    end
end
