clear all; close all; clc

% Determine timestep (hourly)
load_timestep = 8760/12; %monthly
simulation_timestep = 8760/365; %Daily resolution
dur=[31 28 31 30 31 30 31 31 30 31 30 31]'*24;

%waste+DH load
availability=0.15; %[%] availability of heat production
Pprod=20e6/8760; %[kW] actual excess of waste heat
ratio=[0.1 0.2 0.18 0.38 0.55 0.82 0.9 1 0.48 0.35 0.12 0.2]';
step_waste_heat=Pprod*ratio*availability*1000; %W
grafico_WIP=step_waste_heat.*dur;

%load
LOAD2 = load('Real_load_VKBtest.txt'); 
Input2 = load('input.txt');
HL2=LOAD2(:,2)*10^3; %heating load %[W]

%load
LOAD = load('Real_load.txt'); 
Input = load('input.txt');
HL=0.1*LOAD(:,2)*10^3; %heating load %[W]

% load aggregation algorithm--> I have to transform the hourly load information
% into monthly load information
months_matrix=reshape(HL2,[],8760/load_timestep)';
step_DH_load=mean(months_matrix,2);%[W] load from the DH in a step form
grafico_DH=step_DH_load.*dur;

%boreholes
Ny=7; % number of boreholes in Y direction
Nx=7; % number of boreholes in X direction
B=7.5; % boreholes distance
H=150;%[m]--> borehole depth
D=4;% [m]--> buried depth
rb=0.055;%[m] borehole radius
Ntot=Ny*Nx; % Number of boreholes
shank=70e-3;
kfill=0.6; %thermal conductivity of the backfilling material
pipe_do=32e-3;%[m] U-tube external diameter
ro=pipe_do/2;
tt=3e-3; %[m]--> U-tube thickness
pipe_di=pipe_do-2*tt; %[m]--> U-tube pipe internal diameter
ri=pipe_di/2;
pipe_k=0.42;%[W/mK]

%ground
kg=3;
rhog=3800; %kg/m^3
ag=1*10^-6;% m^2/s thermal diffusivity of the ground
cpg=kg/(rhog*ag);


%heat carrier fluid
cpw=4186;%[J/kgK]
rhow=1000;%[kg/m^3]
%evaluation of the equivalent diameter for U-tube pipes
req=(0.414*ro)+(0.5*shank);
Rbeff=abs(1/(2*pi*kfill)*log((rb/req)*sqrt(req/shank/2)));

%% INTIALIZATION 
n = (8760/simulation_timestep);

T_gr = 8;
Tapp=5;
Tb1=zeros(n,1);
Tf1=zeros(n,1);

%% BOREHOLE CONFIGURATION
% Calcola la coordinata X e Y di ciascun borehole
x = 0:B:(Nx-1)*B;
y = 0:B:(Ny-1)*B;

[X, Y] = meshgrid(x, y);

C1 = [];
C2 = [];
d = 0;

%field model
for i = 1:Ny
    for j = 1:Nx
        for z = 1:Ny
            for k = 1:Nx
                if i == z && j == k
                    d = rb^2;
                else
                    d = (X(i, j) - X(z, k))^2 + (Y(i, j) - Y(z, k))^2;
                end
                index = find(C2 == d);
                if isempty(index)
                    C2 = [C2, d];
                    C1 = [C1, 1];
                else
                    C1(index) =  C1(index) + 1;
                end
            end
        end
    end
end

ierf = @(x) (x .* erf(x)) - (1 / sqrt(pi) * (1 - exp(-x.^2)));
hc = @(s) H*s; 
dc = @(s) D*s; 
Ie = @(s) sum(C1 .* exp(-C2 .* s.^2)) / (Nx * Ny);
ILS = @(hc, dc) 2 * ierf(hc) + 2 * ierf(hc + (2 * dc)) - ierf((2 * hc) + (2 * dc)) - ierf(2 * dc);
intfunc = @(s) Ie(s) .* ILS(hc(s), dc(s)) / (H.*(s.^2));

%% HEAT PUMP MODEL
% Loading heat pump model coefficients
Condenser_Coefficients = load('Condenser_capacity.txt');
Evaporator_Coefficients = load('Evaporator_capacity.txt');
Compressor_Coefficients = load('Compressor_capacity.txt');
CondenserTemperature_Coefficients = load('Condenser_temperature.txt');

% Heat pump selection
index = 6; % *** Enter the index of chosen compressor model from Heatpump_info.txt *** % index = 6 for Load 1, index = 10 for Load 2
q = 1*Evaporator_Coefficients(index,:);% Function coefficients for the evaporator model
C = 1*Condenser_Coefficients(index,:); % Function coefficients for the condenser model 
P = 1*Compressor_Coefficients(index,:); % Function coefficients for the compressor model
g = 1*CondenserTemperature_Coefficients(index,:); %function coefficients for the condensing temperature

%HP model
To=20;
Td = 65; % discharge temperature of the compressor
% evaporator capacity function
Qcl_0 = q(1) + q(2).*To + q(3).*To.^2 + q(4).*To.^3 + q(5).*Td.^3 + q(6).*Td + q(7).*Td.*To + q(8).*Td.*To.^2; 
% Condenser capacity function
Qcd_0 = C(1) + C(2).*Td + C(3).*To + C(4).*Qcl_0 + C(5).*Td.*To + C(6).*Qcl_0.*Td + C(7).*To.^2 + C(8).*Qcl_0.*To;
% Compressor power function
Pcom_0 = P(1) + P(2).*To + P(3).*Td + P(4).*Qcl_0.*Td + P(5).*To.^2 + P(6).*Td.*Qcl_0.^2 + P(7).*Qcl_0 + P(8).*To.*Qcl_0.^2;
% Condensation temperature function
Tc = g(1).*Td + g(2).*To + g(3).*To.^2 + g(4) + g(5).*Td.*To + g(6).*To.^3 + g(7).*Td.^2 + g(8).*Td.*To.^2; 
%COP
COP = Qcd_0/ Qcl_0;

%% STORAGE OPERATION

surplus = step_waste_heat - step_DH_load; 

surplus_g2=grafico_WIP-grafico_DH;
%grafico Guelpa
deficit = surplus_g2;
for it = 1:12
    if deficit(it)<0
        deficit(it) = -deficit(it);
    else 
        deficit(it) = 0;
    end
end

excess=surplus_g2;
for it = 1:12
    if excess(it)>0
        excess(it) = excess(it);
    else 
        excess(it) = 0;
    end
end


q=surplus/Ntot/H; %flux per length
q_negative = q <0;
q(q_negative) = q(q_negative)*(1-1/COP); % Adjusted flux in discharging regarding heat pump COP
qp = [q(1);diff(q)]; %picewise load generation
tp = linspace(0, 8760, (8760/load_timestep)+1); %time frames according to load time step

for i = simulation_timestep:simulation_timestep:8760
        is = i/simulation_timestep;
        deltaT = 0;
        index = find(tp >= i, 1)-1;
            
        t = i*3600;
        for j = 1:index
            flux = qp(j);
            t_step = (i-tp(j))*3600;
            gfunc = (flux/ (4 * pi * kg) * integral(intfunc, 1 / sqrt(4 * ag * t_step), inf, 'ArrayValued', true));
            deltaT = deltaT + gfunc;
        end
 
        Tb1(is) = T_gr + deltaT;
        Tf1(is) = Tb1(is) + Rbeff*q(index);
        
end
Tb1 = [T_gr;Tb1];
Tf1 = [T_gr;Tf1];
time = 0:simulation_timestep:8760;

%duration curve on DH loads
time_vector=1:1:365*24;
load_sort1=sort(HL,'descend');
load_sort2=sort(HL2,'descend');

save('temperature_data.mat', 'Tb1', 'Tf1');

%% GRAFICO VERDA

% Calcoliamo il valore medio del surplus per ogni mese
mean_surplus = (max_surplus + min_surplus) / 2;

% Calcoliamo gli errori come differenza tra il valore medio e il massimo/minimo
errors = [mean_surplus - min_surplus, max_surplus - mean_surplus];

% Grafico con barre di errore
figure;
bar(1:n_months, mean_surplus, 'FaceColor', [0.8, 0.8, 1]); % Barre con il valore medio del surplus
hold on;
errorbar(1:n_months, mean_surplus, errors(:, 1), errors(:, 2), '.', 'Color', 'black', 'LineWidth', 1); % Barre di errore

xlabel('Mesi');
ylabel('Surplus');
title('Surplus mensile con barre di errore (massimo e minimo)');
grid on;

% figure(1)
% plot(time,Tb1,'*')
% hold on
% plot(time,Tf1,'*')
% xlabel('Time [h]')
% ylabel('Temperature [°C]')
% legend('Wall temperature','Average fluid temperature','location','best')
% set(gca,'fontsize',30)
% 
% 
% figure(3) 
% bar(deficit*10^-3) %not covered by WIP
% hold on
% bar(excess*10^-3) %used for charging season
% xlabel('Time [month]')
% ylabel('Energy [kWh]')
% legend('Energy during winter not covered by WIP','Energy used for charging the BTES from WIP','location','best')
% set(gca,'fontsize',30)
% energy_deficit=sum(deficit)*10^-3;%kWh
% energy_excess=sum(excess)*10^-3;%kWh
% 
% figure(4)
% subplot(2,1,1)
% plot(time_vector,load_sort1*10^-3,'linewidth',2)
% hold on
% xlabel('Time [h]')
% ylabel('Heat flux [kW]')
% set(gca,'fontsize',30)
% title ('Duration curve Load 1')
% 
% subplot(2,1,2)
% plot(time_vector,load_sort2*10^-3,'linewidth',2)
% hold on
% xlabel('Time [h]')
% ylabel('Heat flux [kW]')
% set(gca,'fontsize',30)
% title ('Duration curve Load 2')
% 
% figure(5)
% subplot(2,1,1)
% plot(time_vector,HL*10^-3,'r');
% title('Temperature profile (Fluid)');
% xlabel('Time [h]')
% ylabel('Heating load [kW]')
% set(gca,'fontsize',30)
% title ('Heating demand 1')
% 
% 
% subplot(2,1,2)
% plot(time_vector,HL2*10^-3,'r')
% xlabel('Time [h]')
% ylabel('Heating load [kW]')
% set(gca,'fontsize',30)
% title ('Heating demand 2')

%losses evaluation
area_bhe=B^2;
VV=area_bhe*Ntot*H;% storage volume
mass_g=rhog*VV;
E_ground=1/20*abs(mass_g*cpg*(Tb1(end)-Tb1(1)))*10^-3;%kJ
E_dis = (sum(surplus(1:3,1)) + sum(surplus(end-1:end,1))) * 3624 * 10^-3;%kWh
E_char=0.5*(sum(surplus(4:10,1)))*5136*10^-3;
Net=E_char-E_dis;
Losses=E_ground-Net;
Percentage_L=(E_ground-Losses)/E_ground*100;



save('loss.mat', 'Percentage_L','E_char','E_dis');




