clear all; close all; clc

% Determine timestep (hourly)
load_timestep = 8760 / 12; % Monthly
simulation_timestep = 8760 / 365; % Daily resolution

% Heat load and waste data
availability = 0.15; % [%] Heat production availability
Pprod = 20e6 / 8760; % [kW] Actual excess of residual 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]'; % Data ratio
step_waste_heat = Pprod * ratio * availability * 1000; % W

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

% Aggregation algorithm to transform hourly load information into monthly load information
months_matrix = reshape(HL,[],8760/load_timestep)';
step_DH_load = mean(months_matrix,2); %[W] Load from the district in step form
step_DH_load(1:3) = 0;
step_DH_load(end-1:end) = 0;

% Boreholes configuration
Ny = 7; % Number of boreholes in the Y direction
Nx = 7; % Number of boreholes in the X direction
B = 7.5; % Borehole distance
H = 150; %[m] Borehole depth
D = 4; %[m] Buried depth
rb = 0.055; %[m] Borehole radius
Ntot = Ny * Nx; % Total number of boreholes
shank = 70e-3;
kfill = 0.6; % Thermal conductivity of the filling 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 internal diameter
ri = pipe_di / 2;
pipe_k = 0.42; %[W/mK]

% Ground
kg = 3; % Ground thermal conductivity
ag = 1 * 10^-6; %[m^2/s] Ground thermal diffusivity

% Heat carrier fluid
cpw = 4186; %[J/kgK]
rhow = 1000; %[kg/m^3]
req = (0.414 * ro) + (0.5 * shank);
Rbeff = abs(1 / (2 * pi * kfill) * log((rb / req) * sqrt(req / shank / 2)));

% Heat pump model
% Model parameters loading

% Initialization
n = 8760 / simulation_timestep;
T_gr = 8;
Tb1 = zeros(n+1, 1);
Tf1 = zeros(n+1, 1);

% Boreholes configuration
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;


time = 0:simulation_timestep:8760; % Time vector
% Operazione di stoccaggio
years = 5; % Numero di anni
Tb = cell(years, 1); % Inizializza celle per contenere temperature del borehole
Tf = cell(years, 1); % Inizializza celle per contenere temperature del fluido

% Ciclo per ogni anno
for year = 1:years
    % Prepara vettori per memorizzare le temperature
    Tb{year} = zeros(n, 1);
    Tf{year} = zeros(n, 1);
    
    % Calcolo del flusso termico
    surplus = step_waste_heat - step_DH_load; 
    q = surplus / Ntot / H; % flusso per lunghezza
    q_negative = q < 0;
    q(q_negative) = q(q_negative) * (1 - 1 / COP); 
    qp = [q(1); diff(q)]; 
    tp = linspace(0, 8760, (8760/load_timestep)+1);
    
    % Ciclo per ogni passo temporale
    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
    
        % Aggiorna le temperature del borehole e del fluido
        Tb{year}(is) = T_gr + deltaT;
        Tf{year}(is) = Tb{year}(is) + Rbeff * q(index);
    end  

    % Aggiungi la temperatura iniziale
    Tb{year} = [T_gr; Tb{year}];
    Tf{year} = [T_gr; Tf{year}];
    T_gr = Tb{year}(end);
end

time = 0:simulation_timestep:8760; % Vettore temporale

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

combined_time = []; % Inizializza un vettore per contenere il tempo combinato
combined_Tb = []; % Inizializza un vettore per contenere la temperatura del borehole combinata
combined_Tf = [];
for year = 1:years
    % Concatena i vettori di tempo e temperatura del borehole
    combined_time = [combined_time, time + (year - 1) * 8760]; % Aggiungi un offset di tempo per ogni anno
    combined_Tb = [combined_Tb; Tb{year}];
    combined_Tf = [combined_Tf; Tf{year}];
end

% Inizializzazione delle variabili per la media annuale e le deviazioni standard
years_dt = datetime(2024:2028, 1, 1);
Tb_yearly_avg = zeros(5, 1);
Tf_yearly_avg = zeros(5, 1);
Tb_std = zeros(5, 1);
Tf_std = zeros(5, 1);

% Calcolo della media annuale e delle deviazioni standard per 5 anni
for i = 1:5
    indices = find(year(years_dt) == (2024 + i - 1)); % Trova gli indici corrispondenti all'anno i
    Tb_yearly_avg(i) = mean(combined_Tb(indices));
    Tf_yearly_avg(i) = mean(combined_Tf(indices));
    Tb_std(i) = std(combined_Tb(indices));
    Tf_std(i) = std(combined_Tf(indices));
end

% Calcolo dei valori massimi e minimi per le bande di errore
Tb_max = Tb_yearly_avg + Tb_std;
Tb_min = Tb_yearly_avg - Tb_std;
Tf_max = Tf_yearly_avg + Tf_std;
Tf_min = Tf_yearly_avg - Tf_std;

% Plot dei grafici con bande di errore
errorbar(2024:2028, Tb_yearly_avg, Tb_std, 'o-', 'LineWidth', 1.5, 'DisplayName', 'Borehole Wall Temperature');
errorbar(2024:2028, Tf_yearly_avg, Tf_std, 's-', 'LineWidth', 1.5, 'DisplayName', 'Average Fluid Temperature');
xlabel('Year');
ylabel('Temperature (°C)');
title('Temperature Profile with Error Bars for 5 Years');
legend('Location', 'best');
grid on;