%GT01sim.m
clear all
global Y2 Y YY tt

%load results/GTopt/optGT
load optGT

format compact
format long g

%Pådrag 
    %u = [Ws_CH4kompressor; Ws_LUFTkompressor]
    %u = [8.2E5; 2.6E8];     % J/s (verdier hentet fra høstprosjektet)
    u = uXGT(1:2);
    u = u(:);
    %u(1) = 1.1.*u(1);   %Sprang i pådrag
    u_0 = u;
    
    %Forstyrrelser
    %d = [p_CH4; p_LUFT; T_CH4; T_LUFT]
    d = [20E5; 1E5; 283; 288];  %Pa (d(1:2)) og K (d(3:4)) (verdier hentet fra Skogestad(2000))
    
%Initialbetingelser i brennkammeret (verdier funnet fra optimalisering)
    %X_BK =[n_O2; n_N2; n_CO2; n_H20(g); n_CH4; T_BK]
    X_BK = uXGT(3:7);   %Utelater n_CH4
    
    %Feil fra reguleringsønsket legges her inn som "tilstander"
    T_BKs = 1773;           %[K], settpunkt for temp.reg. i BK
    F_Fs = 34.0;       %Settpunkt for molforhold mellom fødestrøm LUFT og CH4 
    F_F = 34.0; %Molforhold for initielle fødestrømmer hentet fra opt02
    e_T_BK_int = T_BKs - X_BK(5);  %Integrert reguleringsavvik
    e_F_F_int = F_Fs - F_F;       %Integrert reguleringsavvik          
    
    X_BK = [X_BK; F_F; e_T_BK_int; e_F_F_int];
    X_BK = X_BK(:);              

% ODE15s
    t = 4000;      %s
    [tid,X_BK] = ode15s('GT01',[0 t],X_BK,[],d,T_BKs,F_Fs,u_0)

% Figurer
if 1==0
figure(1)
    subplot(321)
    plot(tid, X_BK(:,1)), ylabel('mol O2'), title('Hold-up BK'), %xlabel('Tid [s]');
    subplot(322)
    plot(tid, X_BK(:,2)), ylabel('mol N2'), title('Hold-up BK'), %xlabel('Tid [s]');
    subplot(323)
    plot(tid, X_BK(:,3)), ylabel('mol CO2'), title('Hold-up BK'), %xlabel('Tid [s]');
    subplot(324)
    plot(tid, X_BK(:,4)), ylabel('mol H2O'), title('Hold-up BK'), %xlabel('Tid [s]');
    subplot(325)
    plot(tid, X_BK(:,5)), xlabel('Tid [s]'), ylabel('Kelvin'), title('Temperatur BK');
end
%Beregner p_BK for forskjellige tidspunkt og plotter det mot tid
    V_BK = 0.4;     %m3
    R = 8.314;      %J/(K*mol)
    n_BKtot = zeros(length(tid),1);        %Gir kolonnevektor av riktig størrelse
    p_BK = zeros(length(tid),1);           %Gir kolonnevektor av riktig størrelse
for i = 1:length(tid)
    n_BKtot(i,1) = sum(X_BK(i,1:4));            %mol, har tidligere satt n_CH4=0
    p_BK(i,1) = n_BKtot(i,1)*R*X_BK(i,5)/V_BK;    %Pa
end

figure(1)
subplot(326)
    plot(tid,p_BK), xlabel('Tid [s]'), ylabel('P(BK), [Pa]'), title('Trykk i BK');

if 1==0
figure(5)
subplot(211)
plot(tid, X_BK(:,7)), ylabel('Avvik, Temp. [K]'), title('Integrert Reguleringsavvik T(BK)'), %xlabel('Tid [s]');
subplot(212)
plot(tid, X_BK(:,8)), xlabel('Tid [s]'), ylabel('Avvik, Molforhold'), title('Integrert Reguleringsavvik F/F');
end

if 1==1
tid1 = tt(1);
yy = YY(1,:);
    for i = [2:length(tt)];
        if tt(i-1) < tt(i);
        tid1 = [tid1; tt(i)];
        yy = [yy; YY(i,:)];   
        end
    end
end
figure(2)
    subplot(321)
    plot(tid1, yy(:,1)), ylabel('Ws,CH4, [W]')%, title('Forstyrrelse T(LUFT) = 288 K')
    subplot(322)
    plot(tid, X_BK(:,6)), ylabel('Fluft/Fch4')
    subplot(323)
    plot(tid1, yy(:,2)), ylabel('Ws,LUFT, [W]')
    subplot(324)
    plot(tid, X_BK(:,5)), ylabel('T BK, [K]')
    subplot(325)
    plot(tid,p_BK), xlabel('Tid [s]'), ylabel('P(BK), [Pa]')
    subplot(326)
    plot(tid1, yy(:,3)), ylabel('J, [NOK/h]'), xlabel('Tid [s]')
    
if 1==0
figure(3)
    subplot(211)
    plot(tid, X_BK(:,6)), ylabel('Fluft/Fch4')%, title('Settpunktsendring for Fluft/Fch4 fra 34.0 til 34.5');
    subplot(212)
    plot(tid1, yy(:,1)), ylabel('Ws,CH4, [W]'), xlabel('Tid [s]')

figure(4)
subplot(211)
plot(tid, X_BK(:,5)), ylabel('T BK, [K]')%, title('Settpunktsendring for T(BK) fra 1773 til 1850 K'), %xlabel('Tid [s]');
subplot(212)
plot(tid1, yy(:,2)), xlabel('Tid, [s]'), ylabel('Ws,LUFT, [W]'), %title('pCH4 = 24 bara');
end
Y2 
    
