% DYNAMIC FUNCTION AND STATE EQUATIONS FOR THE 4:1 HEN

function xprime = Dynamic(t,X,U,N,HXindex)

% Defining the outlet varibles
Th_out = X(1:N);
Twall = X(N+1:2*N);
Tc_out = X(2*N+1:3*N);

% Defining inlet parameters from Simulink
Th_in(1) = U(1);
mh_in = U(2);
Tc_in(1) = U(3);
m0_in = U(4);

% Calling parameters from Data.m file
[T0,Th1,Th2,Th3,Th4,Th5,...
            m0,m1,m2,mh1,mh2,mh3,mh4,mh5,...
            rho_0,hc,Cp0,...
            Vwall,rho_wall,Cp_wall,...
            P1,P2,P3,P4,P5] = Data;
        
if HXindex == 1
    Cph = 2;
    wh = Cph*mh_in;
    rho_h = rho_0;
    hh = 1.2*hc;
    U = (hh*hc)/(hh+hc);
    Vhot = mh_in/rho_h;
    Vcold = m0_in/rho_0;
    w0 = m0_in*Cp0;
    Ai = 19; 
        
 
elseif HXindex == 2
    Cph = 2;
    wh = Cph*mh_in;
    rho_h = rho_0;
    hh = 1.42*hc;
    U = (hh*hc)/(hh+hc);
    Vhot = mh_in/rho_h;
    Vcold = m0_in/rho_0;
    w0 = m0_in*Cp0;
    Ai = 29.5;
    

elseif HXindex == 3
    Cph = 2;
    wh = Cph*mh_in;
    rho_h = rho_0;
    hh = 1.389*hc;
    U = (hh*hc)/(hh+hc);
    Vhot = mh_in/rho_h;
    Vcold = m0_in/rho_0;
    w0 = m0_in*Cp0;
    Ai = 43.7;
    
 
elseif HXindex == 4
    Cph = 2;
    wh = Cph*mh_in;
    rho_h = rho_0;
    hh = 0.70*hc;
    U = (hh*hc)/(hh+hc);
    Vhot = mh_in/rho_h;
    Vcold = m0_in/rho_0;
    w0 = m0_in*Cp0;
    Ai = 103;


else HXindex == 5
    Cph = 2;
    wh = Cph*mh_in;
    rho_h = rho_0;
    hh = 1.43*hc;
    U = (hh*hc)/(hh+hc);
    Vhot = mh_in/rho_h;
    Vcold = m0_in/rho_0;
    w0 = m0_in*Cp0;
    Ai = 38.3;
    
 
end


% STATE EQUATIONS

% Hot stream
dThotdt(1) = (Th_in(1)-Th_out(1)-((U*Ai)/(wh*N))*(Th_out(1)-Twall(N))*...
    (mh_in*N)/(rho_h*Vhot));

% Wall
dTwalldt(1) = (hh*(Th_out(N)-Twall(1))-hc*(Twall(1)-Tc_out(1)))*...
    (Ai/(rho_wall*Cp_wall*Vwall));

% Cold stream
dTcolddt(1) =(Tc_in(1)-Tc_out(1)-((U*Ai)/(w0*N))*(Tc_out(1)-Twall(1)))*...
    ((m0_in*N)/(rho_0*Vcold));


for i = 2:N
    j = N-i+1;
    dThotdt(i) = (Th_out(i-1)-Th_out(i)-((U*Ai)/(wh*N))*...
        (Th_out(i)-Twall(j))*(mh_in*N)/(rho_h*Vhot));
end

for j = 2:N
    i = N-j+1;
    dTwalldt(j) = (hh*(Th_out(i)-Twall(j))-hc*(Twall(j)-Tc_out(j)))*...
        (Ai/(rho_wall*Cp_wall*Vwall));
    dTcolddt(j)=(Tc_out(j-1)-Tc_out(j)-((U*Ai)/(w0*N))*...
        (Tc_out(j)-Twall(j))*((m0_in*N)/(rho_0*Vcold)));
end

xprime = [dThotdt, dTwalldt, dTcolddt];
