function varargout  = reaseprec(un,param)

%Calculate all mass balance derivatives and other balance equations (over
%purge split and recycle mixing) from un (combined vector which contains,
%in order, U, nr and ncol)

un = param.scale_un.*un ; % Get unscaled values of un back

%Extract inputs
u = un(1:11) ; %Inputs to process

Tr = u(1) ; %#ok<*NASGU> % Reactor temperature
RF = u(2)  ;% Flow to reactor
x_rf = u(3:5);  % Composition of reactor feed
RP = u(6) ; % Flow out of reactor
L = u(7)  ;% Column reflux rate
V = u(8)  ;% Column boilup rate
FA = u(9) ; % Fresh A feed rate
P = u(10)  ;% Purge ratio (purge rate/distillate rate)
RCY = u(11) ; % Recycle rate



nr = un(12:12+param.NC-1) ; %Component holdups in reactor
Mr = sum(nr) ;
ncol = un((12+param.NC):end) ; %Holdups in column. Remember this is total holdups and holdups of NC - 1 lightest components.

reactor  %Runs the reactor calculations

ucol = [L V] ;
dcol = [RP 1 xr'] ;

[dndt_col B D xB xD] = column(ncol,ucol,dcol,param) ;


%Collect the mole balance derivatives (for reactor and column)
dndt = [dndt_r ;
        dndt_col];

%Then the other equality constraints (mass balance over purge split, mass balance over feed/recycle mix and total )

xfa = param.xfa ;
res_alg = [(D-RCY-D*P) ; (RF*x_rf - RCY*xD - FA*xfa) ; (RF - RCY - FA)];

eqn = [dndt ; res_alg] ;

switch param.output
    case('Solve_equation') % For use with fsolve
        vars = [Tr Mr V xB(2) FA P]' ;
        iu = 1:6;
        specs = vars(iu) - param.specs(iu) ;
        eqn = [eqn ; specs]  ;
        varargout={eqn} ;
    case('Solve_optimizer') % For optimization with 0 degrees of freedom
        vars = [Tr Mr V xB(2) FA P]' ;
        specs = vars - param.specs ;
        eqn = [eqn ; specs]  ;
        con = [];
        varargout={con,eqn} ;        
    case('Opt_1') %  Optimization with Fmax specified and F < Fmax
        con = [0.9 - xB(2) ;
            Mr - param.Mrmax ;
            Tr - param.Trmax ;
            V - param.Vmax ;
            FA - param.Fmax ;
            -D
            -B
            P*D - param.pmax];
        con = con.*[1 1e-4 1e-2 1e-1 1 1 1 1]';
        varargout={con,eqn} ;
	case('Opt_2') % Optimization with F specified directly
        eqn = [eqn ; (FA - param.Fmax)]  ;
        con = [0.9 - xB(2) ;
            Mr - param.Mrmax ;
            Tr - param.Trmax ;
            V - param.Vmax ;
            FA - param.Fmax ;
            -D
            -B
            P*D - param.pmax];
        con = con.*[1 1e-4 1e-2 1e-1 1 1 1 1]';
        varargout={con,eqn} ;
end