%GT01.m

function f = GT01(t,X_BK,flag,d,T_BKs,F_Fs,u_0)

global Y2 Y YY tt

%Konstanter
    R = 8.3145;                     %J/(K*mol) fra SICD
    Vm = 23.690E-3;                 %m3/mol, molart volum ved 1 atm og 288.7 K
    eta_KOMPR = 0.9;                %Polytropisk kompressorvirkningsgrad (fra Bolland (høstprosjekt))
    eta_GT = 0.9;                   %Adiabatisk turbinvirkningsgrad (fra Bolland (høstprosjekt))
    delta_rH_ref = -803E3;          %J/mol, forbrenningsvarme fra delta_fH i SICD, for T_0=298 K, gjelder rx i BK
    delta_rcp = 11;                 %J/(K*mol), fra cp-verdier i SICD, gjelder rx i BK
    
    Mm = [32.0; 28.0; 44.0; 18.0; 16.0];    %g/mol  Mm =[O2; N2; CO2; H20; CH4] verdier fra SICD
    cp = [29; 29; 37; 34; 36];              %J/(K*mol) cp =[O2; N2; CO2; H20(g); CH4] verdier fra SICD
    x_LUFT = [0.21; 0.79; 0; 0; 0];         % mol% (Luft-sammensetning)
    x_CH4 = [0; 0; 0; 0; 1];                % mol% (Metanføde-sammensetning)
    
    cp_LUFT = cp'*x_LUFT;
    gamma_LUFT = cp_LUFT/(cp_LUFT-R);   %adiabateksponent
    gamma_CH4 = cp(5)/(cp(5)-R);        %adiabateksponent
    
    V_BK = 0.4;                         %m3 BK-volum
    
    %Priser
    P_el = 0.15;                        %NOK/kWh, (verdi fra Bolland)
    P_NG = 0.70;                        %NOK/Sm3, (verdi fra Bolland)
%--------------------------------------------------------------------------------    
 %Regulatorer 
 
    %Beregning av reguleringsavvik
    e_T_BK = T_BKs - X_BK(5);     
    e_F_F = F_Fs - X_BK(6);
    
    e_T_BK_int = X_BK(7);    %Integrert reguleringsavvik
    e_F_F_int = X_BK(8);    %Integrert reguleringsavvik
    
    Ws_CH4 = u_0(1);    %Settes her i tilfelle regulatorene er "slått av"
    Ws_LUFT = u_0(2);
if 1==1
    %Regulering av T_BK
    Kc_T_BK = 2845;
    tau_T_BK = 4.7E-3;
    Ws_LUFT_0 = u_0(2);

    Ws_LUFT = Ws_LUFT_0 - Kc_T_BK*(e_T_BK + e_T_BK_int/tau_T_BK);  %W, dersom T_BK > T_BKs ==> øker Ws_LUFT
    if Ws_LUFT >= 2*Ws_LUFT_0
        Ws_LUFT = 2*Ws_LUFT_0;
    elseif Ws_LUFT <= 0
        Ws_LUFT = 0;
    end
end

if 1==1
    %Regulering av F_F
    Kc_F_F = 55556;
    tau_F_F = 60;
    Ws_CH4_0 = u_0(1);
        
    Ws_CH4 = Ws_CH4_0 + Kc_F_F*(e_F_F + e_F_F_int/tau_F_F); %W, dersom F_F > F_Fs ==> reduseres Ws_CH4
    if Ws_CH4 >= 2*Ws_CH4_0
        Ws_CH4 = 2*Ws_CH4_0;
    elseif Ws_CH4 <= 0
        Ws_CH4 = 0;
    end
end
%--------------------------------------------------------------------------------
%Brennkammer, BK (her beregnes trykket ut av kompressorene)
    %Antar all hold-up i GT skjer i brennkammeret.
    %Antar forbrenningen skjer momentant og fullstendig, og 100% omsetning av CH4
    
    n_BK = [X_BK(1:4); 0];               %mol i BK av komponentene
    T_BK = X_BK(5);                 %K 
    
    n_BKtot = sum(n_BK);            %mol i BK totalt
    x_BK = n_BK/n_BKtot;            % mol% (sammensetning i BK)
    
    p_BK = n_BKtot*R*T_BK/V_BK;     %Pa 
    
    ny = [-2; 0; 1; 2; -1];          %reaksjonskoeffisienter, [O2; N2; CO2; H2O; CH4] 
    
    %Kompressorberegninger
        %Metankompressor
            %Flow og Temp_ut
            F_CH4 = (Ws_CH4*eta_KOMPR)/(cp(5)*d(3)*((p_BK/d(1))^((gamma_CH4 - 1)/gamma_CH4) - 1));   %mol/s
            F_2 = F_CH4;                                            %mol/s (ingen hold-up i kompressor)
        
            T_2 = d(3) + Ws_CH4/(F_CH4*cp(5));     %K, temp. i strøm 2, d(3)=T_CH4
                          
       %Luftkompressor
            %Flow og Temp_ut
            F_LUFT = (Ws_LUFT*eta_KOMPR)/(cp_LUFT*d(4)*((p_BK/d(2))^((gamma_LUFT - 1)/gamma_LUFT) - 1)); %mol/s
            F_1 = F_LUFT;                                           %mol/s (ingen hold-up i kompressor)
            
            T_1 = d(4) + Ws_LUFT/(F_LUFT*cp_LUFT); %K, temp. i strøm 1, d(4)=T_LUFT
            
    F_Fny = F_LUFT / F_CH4;
    
    %Turbin (her beregnes massestrøm ut fra BK)
        p_3_0 = 3.0E6;            %Pa, verdi fra høstprosjektet
        rho_3_0 = 5.75;         %kg/m3, verdi beregnet fra høstprosjekt samt oppholdstid (Bolland)
        w_3_0 = (1.73E6)/3600;  %kg/s, verdi hentet fra høstprosjektet
        p_3 = p_BK;             %Pa
        rho_3 = (n_BK'*Mm)/(V_BK*1000);    %kg/m3    
        w_3 = w_3_0*sqrt((p_3*rho_3)/(p_3_0*rho_3_0));   %kg/s, (fra Hyllseth (perm))
        
    %Brennkammer
        T_ref = 298;                  %K, referansetemperatur for entalpievalueringer
        delta_rH_T_BK = delta_rH_ref + delta_rcp*(T_BK - T_ref);      %J/mol (CH4)
    
        F_3 = (w_3*1000)/(Mm'*x_BK);                                %mol/s
        %Materialbalanse   
        DnDt = F_LUFT*x_LUFT + F_CH4*x_CH4 - F_3*x_BK + F_CH4*ny;   %mol/s
        
        %Energibalanse
        DT_BKDt = (F_1*cp_LUFT*(T_1 - T_BK) + F_2*cp(5)*(T_2 - T_BK) + (-delta_rH_T_BK)*F_CH4)/(n_BK'*cp);
    
        %Molforhold
        DF_FDt = F_Fny - X_BK(6);
        
   %Turbin (her beregnes arbeidet som tas ut i turbinen)
        p_4 = 1E5;                                  %Pa
        gamma_3 = (cp'*x_BK)/((cp'*x_BK) - R);      %adiabateksponent
        Wsrev_GT = F_3*(cp'*x_BK)*T_BK*((p_4/p_3)^((gamma_3 - 1)/gamma_3) - 1); %J/s, reversibelt arbeid
        Ws_GT = Wsrev_GT*eta_GT;                     %J/s, arbeid som tas ut i GT
        
        cp_4 = cp'*x_BK;                            %J/(K*mol)
        T_4 = T_BK + Ws_GT/(F_3*cp_4);               %K, fra likn. i Skogestad (2000)
        
  %Netto el-produksjon
        sigma_el = Ws_GT + Ws_CH4 + Ws_LUFT;     %J/s, '+' er valgt fordi produsert el har negativt fortegn 
      
  %Kostnadsfunksjon
    J = -3600*((P_el/3600)*(-sigma_el/1000) - P_NG*Vm*F_CH4);  %NOK/s
 
  %Reguleringsavvik integreres med tid ved å legges inn i f
        De_T_BKDt = T_BKs - T_BK;
        De_F_FDt = F_Fs - X_BK(6);
                            
   f = [DnDt(1:4); DT_BKDt; DF_FDt; De_T_BKDt; De_F_FDt];
   
  %Y = [sigma_el; Ws_GT; F_CH4; F_LUFT; p_BK; rho_3; F_F];       %Denne tas med bare for å sjekke disse verdiene
  Y2 = [J; T_BK; p_BK; F_CH4; X_BK(6); T_4; w_3];
  Y = [Ws_CH4, Ws_LUFT, J];
  YY = [YY; Y];
  tt = [tt; t];
  return