%GTfmincon.m

function J = GTfmincon(uX,d,uXsc)
global Z

    %Tilstander og pådrag
    uX = uX.*uXsc;      %Skalerer uX tilbake til "normale" fysikalske størrelser
    u = uX(1:2);        %u(1) og u(2) er arbeid tilført henholdsvis metankompressor og luftkompr.
    X_BK = uX(3:7);
    
    %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
    delta_rcp = 11;                 %J/(K*mol), fra cp-verdier i SICD
    
    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)
    
    %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] 
    
    %Her beregnes flow gjennom kompressorene og temp. ut av kompressorene
        %Metankompressor
            %Flow og Temp_ut
            F_CH4 = (u(1)*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) + u(1)/(F_CH4*cp(5));     %K, temp. i strøm 2, d(3)=T_CH4
    
        %Luftkompressor
            %Flow og Temp_ut
            F_LUFT = (u(2)*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 kopressor)
            
            T_1 = d(4) + u(2)/(F_LUFT*cp_LUFT); %K, temp. i strøm 1, d(4)=T_LUFT    
            
    %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);
        
   f = [DnDt(1:4); DT_BKDt];   
   fsc = [1; 1; 1; 1; 1];       %Skaleringsvektor for f
   f = f./fsc;                  %Skalerer f
   
   %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
        
  %Netto el-produksjon
        sigma_el = Ws_GT + sum(u);                  %J/s
        
  %Kostnadsfunksjon
  J = -3600*((P_el/3600)*(-sigma_el/1000) - P_NG*Vm*F_CH4);  % NOK/h
  
  Z = [J; sigma_el; F_CH4; F_LUFT; p_BK; rho_3; uX];       %Denne tas med bare for å sjekke disse verdiene
  return