function [r,ndot] = Rrate_NewDist(T,wNH3,R,in,kval,d)

% This function calculates the reation rate of ammonia in kg NH3/kg cat,s.
%
% The required input is given by
%   T       Temperature of the stage
%   wNH3    Mass fraction of ammonia of the stage
%   R       H2/N2 ratio of the stage
%   in      Array containing the inlet conditions for the pressure
%   kval    Structure containing parameters for the reaction rate constant
%   d       Catalyst activity

% Import of CASADI
import casadi.*

% Calculation of the mass fractions of the streams
if isa(T,'casadi.MX')
    m = MX.zeros(3,numel(wNH3));
else
    m = zeros(3,numel(wNH3));
end
m(1,:) = (1-wNH3').*2.*R'./(28+2*R'); 
m(2,:) = (1-wNH3').*28./(28+2*R'); 
m(3,:) = wNH3';

% Calculation of the mass and molar flow of the streams
mdot = m*in(1);
ndot = mdot./([2 28 17]'*ones(1,numel(wNH3)));

% Calculation of the molar fractions
n = ndot./(ones(3,1)*sum(ndot));

% Calculation of the reaction rate constants
k1  = kval.Afor*exp(-kval.Eafor./(kval.R*(T'+273))); 	% Forward reaction rate constant
k2  = kval.Abac*exp(-kval.Eabac./(kval.R*(T'+273)));	% Backward reaction rate constant

% Calculation of the partial pressure [bar]
pH2  = n(1,:)*in(2);
pN2  = n(2,:)*in(2); 
pNH3 = n(3,:)*in(2) ;

% Calculation of the reaction rate and transformation of it
r = k1.*pN2.*pH2.^1.5./pNH3 - k2.*pNH3./pH2.^1.5;   % [kmol N2/ m3 cat, h]
r = [3*2 28 -2*17]'*r/2200/3600;                    % [kg [H2 N2 NH3]/ kg cat, s]
r = (1+d)*4.75*r;                                   % Multiply by 4.75 to match industrial instability