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;

% Inizializzazione delle variabili per la media mensile e le deviazioni standard
time_dt = datetime(2024,1,1) + hours(time);
Tb1_monthly_avg = zeros(12, 1);
Tf1_monthly_avg = zeros(12, 1);
Tb1_std = zeros(12, 1);
Tf1_std = zeros(12, 1);

% Calcolo della media mensile e delle deviazioni standard
for i = 1:12
    indices = find(month(time_dt) == i); % Trova gli indici corrispondenti al mese i
    Tb1_monthly_avg(i) = mean(Tb1(indices));
    Tf1_monthly_avg(i) = mean(Tf1(indices));
    Tb1_std(i) = std(Tb1(indices));
    Tf1_std(i) = std(Tf1(indices));
end

% Calcolo dei valori massimi e minimi per le bande di errore
Tb1_max = Tb1_monthly_avg + Tb1_std;
Tb1_min = Tb1_monthly_avg - Tb1_std;
Tf1_max = Tf1_monthly_avg + Tf1_std;
Tf1_min = Tf1_monthly_avg - Tf1_std;

% Plot dei grafici con bande di errore
figure;
hold on;

errorbar(1:12, Tb1_monthly_avg, Tb1_std ,'o-','' 'DisplayName', 'Borehole Wall Temperature');
errorbar(1:12, Tf1_monthly_avg, Tf1_std,'s-', 'DisplayName', 'Average Water Temperature');

legend('show');
xlabel('Time [month]');
ylabel('Temperature [°C]');
set(gca,'fontsize',30)




