function f = phimod(input)
%This file computes the cost function of the model for the optimization problem
%(Electrical efficiency).
%INPUT - Inputs - Flow rates of H2, O2 and i (current density)
%OUTPUT - cost to be minimized, states (y) and Ucell at the optimum.
Ncells = 5;
F = 96485;
LHV = 242000; %Lower heating value in J/mol
x0 = [1050]; %Starting guess for fsolve
options = optimset('Tolfun',1e-8,'display','off');
%%%%%%%%%Blower details%%%%%%%%%%%%
%The blower gives a pressure drop of 26 mbar for 60 ml/min/cm^2 of air.
%This is converted into SI units in terms for mol/sec of air. The blower is
%assumed to have a 40% efficiency.
CONV_BLOWER = 4e4*26/33;
BLOWER_EFF = 0.4;
[y,fval,flag,output] = fsolve(@(x) model(x,input),[x0],options);
[fv,U,Nout,power] = model(y,input);
blower = CONV_BLOWER*(input(2)/5*4.76)^2/BLOWER_EFF*5;
cost = -(power-blower)/(input(1)*LHV);
f=cost;
return