% Det dynamiske problemet
%--------------------------------------------------------------------------
% Funskjonene tar inn tidspunktet t, tilstandene T og dataene D og
% returnerer de tidsderiverte. I tillegg lagres en del interne variabler
% globalt.
%--------------------------------------------------------------------------
function dT = model(t,T,D)
global u P w Tventil y

% Gir lokale navn til variablene
%--------------------------------------------------------------------------
Tw_NG1A = T(2:D.n+1);
T_NG1A = [T(1)*ones(D.n,1),T(D.n+2:2*D.n+1),T(2*D.n+2:3*D.n+1),T(3*D.n+2:4*D.n+1),T(4*D.n+2:5*D.n+1)];
T_SW1 = [T(5*D.n+2),T(5*D.n+3)];
T_SW2 = [T(5*D.n+4),T(5*D.n+5)];
ml_NG1A = T(5*D.n+6);
mg_NG1A = T(5*D.n+7);
ml_SW2 = T(5*D.n+8);
mg_SW2 = T(5*D.n+9);
T_tank = T(5*D.n+10:5*D.n+12)';
T_NG1B = [T(5*D.n+13)*ones(D.n,1),T(6*D.n+14:7*D.n+13),T(7*D.n+14:8*D.n+13),T(8*D.n+14:9*D.n+13),T(9*D.n+14:10*D.n+13)];
Tw_NG1B = T(5*D.n+14:6*D.n+13);
ml_NG1B(1) = T(10*D.n+14);
mg_NG1B(1) = T(10*D.n+15);
mg_SW1 = T(10*D.n+16);
mg_tank = T(10*D.n+17:10*D.n+19)';
ml_tank3 = T(10*D.n+20);
I = T(10*D.n+21:10*D.n+23);
%--------------------------------------------------------------------------

% Trykkberegninger
%--------------------------------------------------------------------------
% Gassfraksjon ut hvert kontrollvolum
f_NG1A = [1,1,0,1,1];
f_NG1B = [1,1,0,1,1];
f_SW1 = [1,0];
f_SW2 = [0,0];
f_tank = [1,1,0];

% Gass volumer
Vg_NG1A(1) = D.V_NG1A(1)-ml_NG1A/D.rhol_NG1A(1);
Vg_NG1B(1) = D.V_NG1B(1)-ml_NG1B/D.rhol_NG1B(1);
Vg_SW2(1) = D.V_SW2(1)-ml_SW2/D.rhol_SW2(1);
Vg_tank = [D.V_tank(1:2),D.V_tank(3)-ml_tank3/D.rhol_tank(3)];

% Trykk fra ideell gasslov
P_NG1A(1) = mg_NG1A*D.R*T_NG1A(1,1)/(Vg_NG1A(1)*D.MW_NG1A(1));
P_NG1B(1) = mg_NG1B(1)*D.R*T_NG1B(1,1)/(Vg_NG1B(1)*D.MW_NG1B(1));
P_SW1(1) = mg_SW1*D.R*T_SW1(1)/(D.V_SW1(1)*D.MW_SW1(1));
P_SW2(1) = mg_SW2*D.R*T_SW1(1)/(Vg_SW2*D.MW_SW2(1));
P_tank = mg_tank.*D.R.*T_tank./(Vg_tank.*D.MW_tank);
P_NG1A(3) = P_tank(3);
P_NG1B(3) = P_NG1A(3);
P = [P;P_NG1A(1),P_NG1A(3),P_NG1B(1),P_NG1B(3),P_SW1(1),P_SW2(1),P_tank];
%--------------------------------------------------------------------------

% Forstyrrelse i mengde naturgass
%--------------------------------------------------------------------------
if t<6000
    f1=0;
    f2=0;
elseif t>=6000
    f1=-0.25;
    f2=-0.25*(1-exp(-(t-6000)/300));
end
%--------------------------------------------------------------------------

% Setpunktsendring
%--------------------------------------------------------------------------
if t>=100
    y.s(1) = D.P_tank(1);
    y.s(2) = D.T_ut_NG1A(2);
    y.s(3) = D.T_ut_NG1B(2);
end
%--------------------------------------------------------------------------

% Step i pådrag (til regulatortuning)
%--------------------------------------------------------------------------
% u.n=[u.n;1];
% u.z=[u.z;0.5,0.51];
% f=0;


% Regulering fra ligning \ref{eq:PI}
%--------------------------------------------------------------------------
Kc = -1*[0.000157,0.067,0.078];
taui = [6.6,234,240];
e = y.s-[P_tank(1),T_NG1A(D.n,2),T_NG1B(D.n,2)];
u.n = [u.n;u.n(1,1)+Kc(1)*(e(1)+I(1)/taui(1))];
k=find(u.n<0);u.n(k)=0;
u.z = [u.z;u.z(1,:)+Kc(2:3).*(e(2:3)+I(2:3)'./taui(2:3))];
k=find(u.z>1);u.z(k)=1;k=find(u.z<0);u.z(k)=0;
siste = length(u.n);
%--------------------------------------------------------------------------


% Strømberegninger
%--------------------------------------------------------------------------
% Kompressorstrømmer
Q_komp = u.n(siste)*D.Q_komp; % \ref{eq:Vn}
w_komp = Q_komp.*P_tank(1:2).*D.MW_tank(1:2)./(D.R*T_tank(1:2)); % \ref{eq:kompFlow}
Ws = w_komp.*D.Cpg_tank(1:2).*T_tank(1:2).*(([P_tank(2),P_SW1(1)]./P_tank(1:2)).^(D.R/D.Cpgm_NG1A(1))-1)./D.eff; % \ref{eq:wKompressor}
% Ventilstrømmer \ref{eq:Ventil}
wvalve_NG1A = D.Cv(1)*u.z(siste,1)*sqrt(P_NG1A(3)-P_NG1A(1));
wvalve_NG1B = D.Cv(2)*u.z(siste,2)*sqrt(P_NG1B(3)-P_NG1B(1));

% Strøm ut av kontrollvolumene. Konstante eller fra ligning \ref{eq:Flow}.
w_NG1A(1) = D.k_flow_NG1A*sqrt((P_NG1A(1)-P_tank(2))*(P_NG1A(1)/T_NG1A(1)-P_tank(2)/T_tank(2)));
w_NG1A(2) = D.w_NG1A(2)+f1*D.w_NG1A(2);
w_NG1A(3) = wvalve_NG1A+wvalve_NG1B;
w_NG1A(4:5) = D.w_NG1A(4:5)+f2*D.w_NG1A(4:5);
w_NG1B(1) = D.k_flow_NG1B*sqrt((P_NG1B(1)-P_tank(1))*(P_NG1B(1)/T_NG1B(1)-P_tank(1)/T_tank(1)));
w_NG1B(2) = w_NG1A(2);
w_NG1B(3) = wvalve_NG1B;
w_NG1B(4:5) = w_NG1A(4:5);
w_SW1 = [D.k_flow_SW1*sqrt((P_SW1(1)-P_SW2(1))*(P_SW1(1)/T_SW1(1)+P_SW2(1)/T_SW2(1))),D.w_SW1(2)];
w_SW2 = [D.k_flow_SW2*sqrt((P_SW2(1)-P_tank(3))),D.w_SW2(2)];
w_tank = [w_komp,w_NG1A(3)];
w = [w;w_NG1A,w_NG1B,w_SW1(1),w_SW2(1),w_tank];

% Fordamping/kondensering innen hvert kontrollvolum
wvap_NG1A(1) = w_NG1A(1);
wvap_SW2(1) = -w_SW2(1);
wvap_NG1B(1) = w_NG1B(1);
wvap_tank3 = 0;
%--------------------------------------------------------------------------

% Varmeoverføring
%--------------------------------------------------------------------------
% Interne temperaturer
Tventil_NG1A = 1/(1/D.T0_sat_NG1A-log(P_NG1A(1)/D.P0_sat_NG1A)*D.R/(D.DHvapm_sat_NG1A+(D.Cpl_NG1A(1)-D.Cpg_NG1A(1))*(298-D.T0_sat_NG1A))); % \ref{eq:Tventil}
Tventil_NG1B = 1/(1/D.T0_sat_NG1B-log(P_NG1B(1)/D.P0_sat_NG1B)*D.R/(D.DHvapm_sat_NG1B+(D.Cpl_NG1B(1)-D.Cpg_NG1B(1))*(298-D.T0_sat_NG1B)));
Tventil = [Tventil;Tventil_NG1A,Tventil_NG1B];
T_NG1A(:,1)= [T_NG1A(1,1):(Tventil_NG1A-T_NG1A(1,1))/(D.n-1):Tventil_NG1A]'; % Konstruerer lineært profil
T_NG1B(:,1)= [T_NG1B(1,1):(Tventil_NG1B-T_NG1B(1,1))/(D.n-1):Tventil_NG1B]';
T_komp2 = Ws/(w_komp(2)*D.Cpg_tank(2))+T_tank(2); % \ref{eq:TKompressor}
DTlm_SW1 = ((T_komp2-T_SW1(2))-(T_SW1(1)-D.T_inn_SW1(2)))/log((T_komp2-T_SW1(2))/(T_SW1(1)-D.T_inn_SW1(2)));; % \ref{eq:DTlm}
% Varmeoverføring
Q_NG1A = ones(D.n,1)*D.hA_NG1A/D.n.*(Tw_NG1A*ones(1,5)-T_NG1A); % \ref{eq:QLinProf}
Q_NG1B = ones(D.n,1)*D.hA_NG1B/D.n.*(Tw_NG1B*ones(1,5)-T_NG1B);
Q_SW1(1) = -D.UA_SW1*DTlm_SW1; % \ref{eq:QSW1}
Q_SW1(2) = -Q_SW1(1);
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% Energibalanser
%--------------------------------------------------------------------------

% Entalpi ut av hvert kontrollvolum (om ikke annet er indikert).
% Ligninger gitt i bilag  \ref{bilag:Entalpi}
%--------------------------------------------------------------------------
h_NG1A = ones(20,1)*f_NG1A.*(ones(20,1)*D.Cpg_NG1A.*(T_NG1A-298)+ones(20,1)*D.DHvap_NG1A)+...
    ones(20,1)*(ones(1,5)-f_NG1A).*(ones(20,1)*D.Cpl_NG1A).*(T_NG1A-298);
h_NG_inn = D.Cpg_NG1A(2)*(D.T_inn_NG1A(2)-298)+D.DHvap_NG1A(2);
h_MK_inn = D.Cpg_NG1A(4)*(D.T_inn_NG1A(4)-298)+D.DHvap_NG1A(4);
h_UK_inn = D.Cpg_NG1A(5)*(D.T_inn_NG1A(5)-298)+D.DHvap_NG1A(5);
h_NG1B = ones(20,1)*f_NG1B.*(ones(20,1)*D.Cpg_NG1B.*(T_NG1B-298)+ones(20,1)*D.DHvap_NG1B)+...
    ones(20,1)*(ones(1,5)-f_NG1B).*(ones(20,1)*D.Cpl_NG1B).*(T_NG1B-298);
h_SW1 = f_SW1.*(D.Cpg_SW1.*(T_SW1-298)+D.DHvap_SW1)+(ones(1,2)-f_SW1).*D.Cpl_SW1.*(T_SW1-298);
h_SW2 = f_SW2.*(D.Cpg_SW2.*(T_SW2-298)+D.DHvap_SW2)+(ones(1,2)-f_SW2).*D.Cpl_SW2.*(T_SW2-298);
h_SW_inn = D.Cpl_SW1(2)*(D.T_inn_SW1(2)-298);
h_tank = f_tank.*(D.Cpg_tank.*(T_tank-298)+D.DHvap_tank)+(ones(1,3)-f_tank).*D.Cpl_tank.*(T_tank-298);
h_komp = h_tank(1:2)+Ws./w_komp;

% $dT/dt$
%--------------------------------------------------------------------------
% NG1A. Ligning \ref{eq:EBNG1A}
dT(1,1) = (wvalve_NG1A*h_NG1A(D.n,3)-w_NG1A(1)*h_NG1A(1,1)+sum(Q_NG1A(:,1)))...
    ./(mg_NG1A(1)*D.Cpg_NG1A(1)+ml_NG1A(1)*D.Cpl_NG1A(1));% 1
dT(2:D.n+1,1) = -sum(Q_NG1A,2)*D.n/(D.Cps_w(1)*D.m_w(1));% Vegg
dT(D.n+2:2*D.n+1,1) = (w_NG1A(2)*[h_NG_inn;h_NG1A(1:D.n-1,2)]-w_NG1A(2)*h_NG1A(:,2)+Q_NG1A(:,2))*D.n...
    ./(D.mg_NG1A(2)*D.Cpg_NG1A(2)+D.ml_NG1A(2)*D.Cpl_NG1A(2));%2
dT(2*D.n+2:3*D.n+1,1) = (w_tank(3)*[h_tank(3);h_NG1A(1:D.n-1,3)]-w_NG1A(3)*h_NG1A(:,3)+Q_NG1A(:,3))*D.n...
    ./(D.mg_NG1A(3)*D.Cpg_NG1A(3)+D.ml_NG1A(3)*D.Cpl_NG1A(3));%3
dT(3*D.n+2:4*D.n+1,1) = (w_NG1A(4)*[h_MK_inn;h_NG1A(1:D.n-1,4)]-w_NG1A(4)*h_NG1A(:,4)+Q_NG1A(:,4))*D.n...
    ./(D.mg_NG1A(4)*D.Cpg_NG1A(4)+D.ml_NG1A(4)*D.Cpl_NG1A(4));%4
dT(4*D.n+2:5*D.n+1,1) = (w_NG1A(5)*[h_UK_inn;h_NG1A(1:D.n-1,5)]-w_NG1A(5)*h_NG1A(:,5)+Q_NG1A(:,5))*D.n...
    ./(D.mg_NG1A(5)*D.Cpg_NG1A(5)+D.ml_NG1A(5)*D.Cpl_NG1A(5));%5
% SW1. Ligning \ref{eq:EBSW1}
dT(5*D.n+2,1) = (w_komp(2)*h_komp(2)-w_SW1(1)*h_SW1(1)+Q_SW1(1))./(D.mg_SW1(1)*D.Cpg_SW1(1)+D.ml_SW1(1)*D.Cpl_SW1(1));%6
dT(5*D.n+3,1) = (w_SW1(2)*(h_SW_inn-h_SW1(2))+Q_SW1(2))./(D.ml_SW1(2)*D.Cpl_SW1(2));%7
% SW2. Satt konstante i steden for å bruke rask regulering
dT(5*D.n+4,1) = 0;
dT(5*D.n+5,1) = 0;
% Tanker. Liging \ref{eq:EBtanker}
dT(5*D.n+10,1) = (w_NG1B(1)*h_NG1B(1)-w_tank(1)*h_tank(1))/(D.ml_tank(1)*D.Cpl_tank(1)+D.mg_tank(1)*D.Cpg_tank(1));%1
dT(5*D.n+11,1) = (w_komp(1)*h_komp(1)+w_NG1A(1)*h_NG1A(1,1)-w_tank(2)*h_tank(2))/(D.ml_tank(2)*D.Cpl_tank(2)+D.mg_tank(2)*D.Cpg_tank(2));%2
dT(5*D.n+12,1) = (w_SW2(1)*h_SW2(1)-w_tank(3)*h_tank(3))/(D.ml_tank(3)*D.Cpl_tank(3)+D.mg_tank(3)*D.Cpg_tank(3));%3
% NG1B. Ligingn \ref{eq:EBNG1B}
dT(5*D.n+13,1) = (wvalve_NG1B*h_NG1B(D.n,3)-w_NG1B(1)*h_NG1B(1,1)+sum(Q_NG1B(:,1)))...
    ./(mg_NG1B(1)*D.Cpg_NG1B(1)+ml_NG1B(1)*D.Cpl_NG1B(1));% 1
dT(5*D.n+14:6*D.n+13) = -sum(Q_NG1B,2)*D.n/(D.Cps_w(2)*D.m_w(2));% Vegg
dT(6*D.n+14:7*D.n+13) = (w_NG1A(2)*[h_NG1A(D.n,2);h_NG1B(1:D.n-1,2)]-w_NG1B(2)*h_NG1B(:,2)+Q_NG1B(:,2))*D.n...
    ./(D.mg_NG1B(2)*D.Cpg_NG1B(2)+D.ml_NG1B(2)*D.Cpl_NG1B(2));%2
dT(7*D.n+14:8*D.n+13) = ((w_NG1A(3)-wvalve_NG1A)*[h_NG1A(D.n,3);h_NG1B(1:D.n-1,3)]-w_NG1B(3)*h_NG1B(:,3)+Q_NG1B(:,3))*D.n...
    ./(D.mg_NG1B(3)*D.Cpg_NG1B(3)+D.ml_NG1B(3)*D.Cpl_NG1B(3));%3
dT(8*D.n+14:9*D.n+13) = (w_NG1A(4)*[h_NG1A(D.n,4);h_NG1B(1:D.n-1,4)]-w_NG1B(4)*h_NG1B(:,4)+Q_NG1B(:,4))*D.n...
    /(D.mg_NG1B(4)*D.Cpg_NG1B(4)+D.ml_NG1B(4)*D.Cpl_NG1B(4));%4
dT(9*D.n+14:10*D.n+13) = (w_NG1A(5)*[h_NG1A(D.n,5);h_NG1B(1:D.n-1,5)]-w_NG1B(5)*h_NG1B(:,5)+Q_NG1B(:,5))*D.n...
    ./(D.mg_NG1B(5)*D.Cpg_NG1B(5)+D.ml_NG1B(5)*D.Cpl_NG1B(5));%5

% $dm/dt$
%--------------------------------------------------------------------------
% NG1A. Ligning \ref{eq:MBNG1A}
dml_NG1A(1) = wvalve_NG1A-wvap_NG1A(1);
dmg_NG1A(1) = wvap_NG1A(1)-w_NG1A(1);
% SW2. Ligning \ref{eq:MBSW2}
dml_SW2(1) = -wvap_SW2(1)-w_SW2(1);
dmg_SW2(1) = w_komp(2)+wvap_SW2(1);
% NG1B. Ligning \ref{eq:MBNG1B}
dml_NG1B(1) = wvalve_NG1B-wvap_NG1B(1);
dmg_NG1B(1) = wvap_NG1B(1)-w_NG1B(1);
% SW1. Ligning \ref{eq:MBSW1}
dmg_SW1 = w_komp(2)-w_SW1(1);
% Tanker. Ligning \ref{eq:MBtanker}
dmg_tank = [w_NG1B(1)-w_tank(1),w_komp(1)+w_NG1A(1)-w_komp(2),wvap_tank3];
dml_tank3 = w_SW2(1)-wvap_tank3-w_tank(3);
%--------------------------------------------------------------------------

dT(5*D.n+6,1)=dml_NG1A;
dT(5*D.n+7,1)=dmg_NG1A;
dT(5*D.n+8,1)=dml_SW2;
dT(5*D.n+9,1)=dmg_SW2;
dT(10*D.n+14) = dml_NG1B(1);
dT(10*D.n+15) = dmg_NG1B(1);
dT(10*D.n+16) = dmg_SW1;
dT(10*D.n+17:10*D.n+19) = dmg_tank';
dT(10*D.n+20)=dml_tank3;
dT(10*D.n+21:10*D.n+23) = e;

% Returnerer de deriverte
%--------------------------------------------------------------------------
dT = dT;
%--------------------------------------------------------------------------