%RSR_PLOT_STUFF

%Run the rsr model with specifications on the variables assumed to be
%active constraints, and plot the remaining variables of interest

%Clear stuff to prepare for the fun
close all
clear all
clc

%Load parameters and initial values
load param.mat
load un0.mat

param.scale_un = abs(un0) ;

param.specs = [390 11000 30 0.900 1 0.2]' ; %Tr, Mr, V,xB(2), FA, P/D 
un0_opt = un0./param.scale_un ;
param.prices = [1 0.07 -2 -0.5]' ;

spec = 20:0.5:30 ; % Here, specify a vector of values to be used for one of the specifications.
d = zeros(size(spec)); % Distillate flowrate, mol/s
b = zeros(size(spec)); % Bottoms (main product) flowrate, mol/s
l = zeros(size(spec)); % Reflux flowrate, mol/s
v = zeros(size(spec)); % Boilup flowrate, mol/s
R = zeros(size(spec)); % Recycle flowrate, mol/s
FR = zeros(size(spec)); % Reactor feed flowrate, mol/s
PR = zeros(size(spec)); % Reactor exit flowrate, mol/s
X = zeros(size(spec)); % Reactant conversion %
S = zeros(size(spec)); % Relative reaction selectivity s_A/s_B
J = zeros(size(spec)); % Objective function


for i=1:length(spec)
    param.specs(3) = spec(i) ; %first tell which specification is to be replaced
    param.output =  'Solve_equation' ; %Use fmincon with 0 degrees of freedom to solve equations
    options_I = optimset('Algorithm','levenberg-marquardt','MaxIter',1e3,'MaxFunEvals',1e7,'Display','off','TolX',1e-5,'TolFun',1e-5);
    [unopt,fval,exitflag] = fsolve(@reaseprec,un0_opt,options_I,param) ;
    un = unopt.*param.scale_un ;
    %Now we need to calculate stuff.
    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) ;
    d(i) = D ;
    b(i) = B ;
    l(i) = L ;
    v(i) = V ;
    R(i) = RCY ;
    FR(i) = RF ;
    PR(i) = RP ;
    X(i) = conv_A ;
    S(i) = s_AB ;
    J(i) = rsr_obj(unopt,param) ;
    disp(['Calculation number ' num2str(i) ' completed'])
end


figure(1),plot(spec,d),ylabel('Distillate flowrate');
figure(2),plot(spec,b),ylabel('Bottoms flowrate');
figure(3),plot(spec,l),ylabel('Reflux flowrate');
figure(4),plot(spec,v,'LineWidth',3.0),xlabel('Feed rate (mol/s)','FontSize',20),ylabel('Column boilup (mol/s)','FontSize',20),set(gca,'LineWidth',2.0,'FontSize',16.0)
figure(5),plot(spec,R),ylabel('Recycke flowrate');
figure(6),plot(spec,FR),ylabel('Reactor feed flowrate');
figure(7),plot(spec,PR),ylabel('Reactor exit flowrate');
figure(8),plot(spec,X),ylabel('% conversion in reactor','FontSize',20);
figure(9),plot(spec,S),ylabel('Relative selectivity in reactor');
figure(10),plot(spec,J),ylabel('Objective function value');