clc
clear all
close all

% ht = .076; %3 inches in meters
% rhob = 727; %Density of gasoline
% rhot = 670; %Density of feed
% MWb = 103; %Molecular weight of gasoline
% MWg = 48.14; %Molecular weight of propane/butane mix
% MWt = (MWb+MWg)/2; %Avg molecular weight
% Dk = 2.5; %Diameter of column
% cb = 2.22; %Specific heat of gasoline
B = 110.92; %Molar flow rate of bottoms
%F = 15480; %Mass flow rate of feed
D = 92.76; %Molar flow rate of distillate
xF = .341; %Liquid fraction of feed
%Increasing xF to 0.65 helps with the final value of the bottoms and the
%distillate, but does not help with the overall shape
LF = 104.25; %Molar liquid flow rate of feed
VF = 98.52; %Molar vapor flow rate of feed
% R = 5; %Reflux ratio
L = 75.64; %Molar flow rate of liquid condensate below feed
% Hnb = .386; % Forgot how we calculated this, I'll add later
Mb = 31.11; %Holdup on bottom tray
M = 5.8; %Holdup on every tray
Md = 13.07; %Holdup on distallate drum
% lamdab = -1.2552; %Latent heat of vaporization for distillate
% tb = 144; %Temperature of bottoms
% tf = 118; %Temperature of feed
% lamda = lamdab*tf; %Latent heat of vaporization at feed temperature
V = 66.34; %Molar flow rate of vapor above feed
alpha = 5.68; %Used for finding VLE
yF = alpha.*xF./(1+(alpha-1).*xF);

%V and L are to be added to VF and LF above and below the feed,
%respectively. Otherwise, the flow rate should just be VF and LF


%Model the equation as x' = A1*x + A2*y(x) + A3
%All values are taken from model equations in the paper
A1 = zeros(16,16);
A2 = zeros(16,16);

A1(1,1) = -B/Mb;
A1(1,2) = (L+LF)/Mb;
for i = 2:7
    A1(i,i) = -(L+LF)/M;
    A1(i,i+1) = (L+LF)/M;
end
A1(8,8) = -(L+LF)/M;
A1(8,9) = L/M;
for i = 9:15
    A1(i,i) = -L/M;
    A1(i,i+1) = L/M;
end
A1(16,16) = -(L+D)/Md;

A2(1,1) = -V/Mb;
for i = 2:8
    A2(i,i-1) = V/M;
    A2(i,i) = -V/M;
end
A2(9,8) = V/M;
A2(9,9) = -(V+VF)/M;
for i = 10:15
    A2(i,i-1) = (V+VF)/M;
    A2(i,i) = -(V+VF)/M;
end
A2(16,15) = (V+VF)/Md;

A3 = zeros(16,1);
A3(8) = LF*xF/M;
A3(9) = VF*yF/M;

x0 = zeros(16,1);

fun = @(x)[(dot(A1(1,:),x)+dot(A2(1,:),alpha.*x./(1+(alpha-1).*x))+A3(1));...
    (dot(A1(2,:),x)+dot(A2(2,:),alpha.*x./(1+(alpha-1).*x))+A3(2));...
    (dot(A1(3,:),x)+dot(A2(3,:),alpha.*x./(1+(alpha-1).*x))+A3(3));...
    (dot(A1(4,:),x)+dot(A2(4,:),alpha.*x./(1+(alpha-1).*x))+A3(4));...
    (dot(A1(5,:),x)+dot(A2(5,:),alpha.*x./(1+(alpha-1).*x))+A3(5));...
    (dot(A1(6,:),x)+dot(A2(6,:),alpha.*x./(1+(alpha-1).*x))+A3(6));...
    (dot(A1(7,:),x)+dot(A2(7,:),alpha.*x./(1+(alpha-1).*x))+A3(7));...
    (dot(A1(8,:),x)+dot(A2(8,:),alpha.*x./(1+(alpha-1).*x))+A3(8));...
    (dot(A1(9,:),x)+dot(A2(9,:),alpha.*x./(1+(alpha-1).*x))+A3(9));...
    (dot(A1(10,:),x)+dot(A2(10,:),alpha.*x./(1+(alpha-1).*x))+A3(10));...
    (dot(A1(11,:),x)+dot(A2(10,:),alpha.*x./(1+(alpha-1).*x))+A3(11));...
    (dot(A1(12,:),x)+dot(A2(12,:),alpha.*x./(1+(alpha-1).*x))+A3(12));...
    (dot(A1(13,:),x)+dot(A2(13,:),alpha.*x./(1+(alpha-1).*x))+A3(13));...
    (dot(A1(14,:),x)+dot(A2(14,:),alpha.*x./(1+(alpha-1).*x))+A3(14));...
    (dot(A1(15,:),x)+dot(A2(15,:),alpha.*x./(1+(alpha-1).*x))+A3(15));...
    (dot(A1(16,:),x)+dot(A2(16,:),alpha.*x./(1+(alpha-1).*x))+A3(16))];

tstep = 3600; %Number of steps
tfinal = 4; %Total run time in hours
%Runge-Kutta Parameters
a1 = 1/2;
a2 = 1/2;
p = 1;
q = 1;
h = tfinal/tstep; %Timestep in hours
xold = x0; %Use initial conditions
xfinal = zeros(16,tstep+1);
xfinal(:,1) = x0;
time = linspace(0,tfinal,tstep+1);
for i = 1:tstep
    xnew = xold + a1*h*fun(xold)+a2*h*fun(xold+h*fun(xold)); %One interation of Runge-Kutta to get new value
    xold = xnew; %Reset old value
    xfinal(:,i+1) = xnew; %Store new value
end
x1final = xfinal(1,:);
x2final = xfinal(2,:);
x3final = xfinal(3,:);
x4final = xfinal(4,:);
x5final = xfinal(5,:);
x6final = xfinal(6,:);
x7final = xfinal(7,:);
x8final = xfinal(8,:);
x9final = xfinal(9,:);
x10final = xfinal(10,:);
x11final = xfinal(11,:);
x12final = xfinal(12,:);
x13final = xfinal(13,:);
x14final = xfinal(14,:);
x15final = xfinal(15,:);
x16final = xfinal(16,:);
plot(time,x1final,time,x2final,time,x3final,time,x4final,time,x5final,...
    time,x6final,time,x7final,time,x8final,time,x9final,time,x10final,...
    time,x11final,time,x12final,time,x13final,time,x14final,time,...
    x15final,time,x16final)
legend('Bottoms', 'Tray 1', 'Tray 2', 'Tray 3', 'Tray 4', 'Tray 5',...
    'Tray 6', 'Tray 7', 'Tray 8', 'Tray 9', 'Tray 10', 'Tray 11',...
    'Tray 12', 'Tray 13', 'Tray 14', 'Distillate', 'location',...
    'northwest', 'NumColumns', 2)