%% Imports for leakage calculation

datafile = "C:\Users\ehtor\Downloads\Datasett_1_K6.xlsx";

O2aveHot = trimmean(xlsread(datafile, 'O2', 'D3:D21600'), 10)/100; %Average value of O2-sensors before Eco
O2Cold = trimmean(xlsread(datafile, 'O2', 'E3:E21600'), 10)/100; %Value of O2-sensor at the chimney

[num, txt, timeDate_raw] = xlsread(datafile, 'O2', 'A3:A21600'); %Reads date and time into the variable timeDate_raw
for i = 781:1440:21598
    timeDate_raw(i) = strcat(timeDate_raw(i), ' 00:00:00'); %Fixing bug in date and time import
end
dateFormat = 'dd/mm/yyyy HH:MM:SS'; %Specify format of date and time import
dateConv = datenum(timeDate_raw, dateFormat); %Convert text string into date numbers

%% Calculate the percentage of air that has leaked into the flue gas, x_air
%Need O2aveHot, O2Cold

%x_air = n_air/n_flue, n_air + n_origFlue = n_flue, 
%O2aveHot*n_origFlue + 0.21*n_air = O2Cold*n_flue
%O2aveHot*(n_flue-n_air) + 0.21*n_air = O2Cold*n_flue
%n_air*(0.21-O2aveHot) = n_flue*(O2Cold-O2aveHot)
%x_air = n_air/n_flue = (O2Cold-O2aveHot)/(0.21-O2aveHot)

x_air = (O2Cold-O2aveHot)/(0.21-O2aveHot);

%% Imports and specifications for Energy loss calculations

T_kjelrom = trimmean(xlsread(datafile, 'Sheet1', 'AG3:AG21600'), 10); %Temperature outside the boiler
T_elfilt = trimmean(xlsread(datafile, 'Sheet1', 'M3:M21600'), 10); %Temperature after EL-filter

Cp_air = (Cp_air_T(T_elfilt+273.15) + Cp_air_T(T_kjelrom+273.15))/2; %[J/K*mol]
%Assume Cp_air is constant and equal to the average of the Cp at 36 and 148 degC

%n_air = x_air*n_flue, n_flue = m_flue/M_flue
M_flue = dot([4.5 11.28 62.47 21 0.74]*(1-x_air)/100+[0.21 0 0.78 0 0.01]*x_air, [32 44.01 28.01 18.02 39.95])/1000;
%[kg/mol}
m_flue = 29.3*1.2922; %[kg/s] Flue gas flow. Density calculated by using 
n_air = x_air*m_flue/M_flue; %[mol/s]

%% Energy loss because of heating of extra air, E_heatLoss

E_heatLoss = n_air*Cp_air*(T_elfilt-T_kjelrom)/(10^6); %[MJ/s]=[MW]

%% Economic loss because of heating of extra air, Monthly_Loss

El_price = mean(xlsread("C:\Users\ehtor\OneDrive\Dokumenter\MasterNorskeSkog\ELKjel.xlsx", 'Sheet1', 'E1529:E1544'));
El_price_year = mean(xlsread("C:\Users\ehtor\OneDrive\Dokumenter\MasterNorskeSkog\ELKjel.xlsx", 'Sheet1', 'E1266:E1630'));
%Average electricity price over the data period (21.feb-8.mar)
Monthly_Loss = E_heatLoss*24*31*El_price;
Yearly_Loss = E_heatLoss*24*365*El_price_year;

%% Imports and specifications for Shunt pump calcultaions
Tw_in = trimmean(xlsread(datafile, 'Sheet1', 'W3:W21600'), 10); %[degC] Water after mix with recycle, before eco 
Tw_feed = trimmean(xlsread(datafile, 'Sheet1', 'U3:U21600'), 10); %[degC] Water temp before mix
Tf_out = trimmean(xlsread(datafile, 'Sheet1', 'AF3:AF21600'), 10); %[degC] Flue temp after eco

m_feed = trimmean(xlsread(datafile, 'Sheet1', 'AE3:AE21600'), 10); %[kg/s] Feed water flow
m_recycle = 6.1; %[kg/s] Water recycle flow

Cp_w = 4.35; %[kJ/kg*K] Assume Cp_w = constant throughout temperature range (ca. 110-215 degC)
Cp_f = (Cp_f_T(380+273.15, x_air) + Cp_f_T(Tf_out+273.15, x_air))/2; %[kJ/kg*K]
%Assume Cp_f is constant and equal to the average of the Cp at 380 and 250 degC

%% Calculate outlet temperature of water, Tw_out
%Need Tw_in, Tw_feed, m_feed, m_recycle

%(Tw_out-Tw_feed)*Cp_w*m_recycle = (Tw_in-Tw_feed)*Cp_w*(m_recycle + m_feed)
%Tw_out*m_recycle = Tw_in*m_recycle + (Tw_in-Tw_feed)*m_feed

Tw_out = Tw_in + (Tw_in-Tw_feed)*m_feed/m_recycle;

%% Calculate inlet temperature of Flue gas, Tf_in
%Need: Cp_w, Cp_f, m_feed+m_recycle, m_flue, Tw_in, Tw_out, Tf_out

%(Tw_out-Tw_in)*Cp_w*(m_feed + m_recycle) = (Tf_in-Tf_out)*Cp_f*m_flue

Tf_in = Tf_out + (Tw_out-Tw_in)*(Cp_w/Cp_f)*(m_feed + m_recycle)/m_flue;

%% Calculate UA
%Need: Cp_w, m_feed+m_recycle, Tw_in, Tw_out, Tf_in, Tf_out

%(Tw_out-Tw_in)*Cp_w*(m_feed + m_recycle) =
%UA*((Tf_in-Tw_out)-(Tf_out-Tw_in))/ln((Tf_in-Tw_out)/(Tf_out-Tw_in))

UA = (Tw_out-Tw_in)*Cp_w*(m_feed + m_recycle)/deltaTm(Tf_in, Tf_out, Tw_in, Tw_out);

%% Calculate theoretical Tw_out and Tf_out when shunt pump is turned off
%Need: Tw_feed, Tf_in, Cp_w, Cp_f, m_feed, m_flue

syms Tw_out_new Tf_out_new

eqns = [Cp_w*m_feed*(Tw_out_new-Tw_feed) == UA*deltaTm(Tf_in, Tf_out_new, Tw_feed, Tw_out_new), ...
    Cp_w*m_feed*(Tw_out_new-Tw_feed) == Cp_f*m_flue*(Tf_in-Tf_out_new)];
%eqns = [Tw_out_new == 5, Tf_out_new == Tw_out_new + 10]
   
vars = [Tw_out_new, Tf_out_new];
init_guess = [Tw_out, Tf_out];

[T_water, T_flue] = vpasolve(eqns, vars, init_guess)

%% Energy loss because of heating of extra air, E_heatLoss_shunt

E_heatLoss_shunt = m_flue*Cp_f*(Tf_out-T_flue)*10^(-3); %[MW] 

%% Economic loss because of running the shunt pump, Monthly_Loss
E_pump = 3e-3; %[MW], shunt pump energy demand
%Average electricity price over the data period (21.feb-8.mar)
Monthly_Loss_shunt = (E_heatLoss_shunt + E_pump)*24*31*El_price;



