clear all; 
close all; 
clc

%waste+DH load
availab=[1.4 0.7 2.2 7.5 18.1 23.5 37.6 28.1 9.0 12.6 5.0 0.6]*10^3./days_per_month;%[kW]--> excess heat not used by the DH and so it is used to charge the borehole
days_per_month=[31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31];
result_vector = [];
for i = 1:length(days_per_month)
    monthly_data=repmat(availab(i),[days_per_month(i)*24,1]);
    result_vector = [result_vector;monthly_data];
end
result_vector = repmat(result_vector, 1, ceil(8760 / length(result_vector)));
result_vector = result_vector(1:8760)';

%load
LOAD = load('Real_load_VKBtest.txt'); 
Input = load('input.txt');
HL=LOAD(:,2); %heating load %[kW]

%boreholes
Ny=25; % number of boreholes in Y direction
Nx=10; % number of boreholes in X direction
B=7.5; % Distanza tra i boreholes (m)
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;
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;
ag=1*10^-6;% m^2/s thermal diffusivity of the ground
qg=30;%[W/m]


%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)));

% Time discretization (e.g., hourly for one charging season)
dt = 3600;
days = 365; 
time = 1:1:days*24;
n=length(time);
 
%% 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));

%% INTIALIZATION 
mdot_BTES_array=[]; 
mdot_WASTE_array=[];
Qborehole_array=[];
Rtot_array=[];
Tf_array=[];
Tin_BTES_array=[];
Tout_BTES_array=[];
Tin_evap_array=[];
Tout_evap_array=[];
QHEX1_array=[];
Tin_WASTE_array=[];
Tout_WASTE_array=[];
deltaT_array=[];
T_gr_array=[];
Tb_array=[];
T_supp_array=[];
T_ret_array=[];
mdot_DH_array=[];
T2=zeros(n,1);
T1=zeros(n,1);
mdot_2=zeros(n,1);
mdot_1=zeros(n,1);
mdot_1_dis=zeros(n,1);
mdot_HP=zeros(n,1);
Tapp=5;
Tin_BTES=20;


%% 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 = 10; % *** 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

%% STORAGE OPERATION
for t=1:n
    % CHARGING SEASON
    if  t>=(120*24) && t<=(273*24) % days between may and september, times 24 hours ---> charging
       
        if result_vector(t)/H>=qg
            ground_flux=qg;
            Qborehole=qg*H;
        else
            ground_flux=result_vector(t)/H;
            Qborehole=result_vector(t);
        end
        deltaT= (ground_flux/ (4 * pi * kg) * integral(intfunc, 1/sqrt(4*ag*time(t)),inf, 'ArrayValued', true));
        Rtot=Rbeff+(deltaT/ground_flux);
        T_gr=8;
        Tb=T_gr+deltaT;    
        Tf = T_gr+(Rtot*ground_flux);
        mdot_BTES=(ground_flux*H)./(2*cpw*abs(Tin_BTES-Tf));
        Tout_BTES=Tf-(ground_flux*H)/(2*cpw*mdot_BTES);
        Tin_WASTE=50;
        Tout_WASTE=95;
        mdot_WASTE=result_vector(t).*1000/(cpw*(Tout_WASTE-Tin_WASTE)); 

        % HEX1
        T_supp=40;
        T_ret=25;
        mdot_DH=HL(t)*1000/(cpw*(T_supp-T_ret));
        T1=T_ret+Tapp;
        mdot_1(t)=mdot_DH.*((T_supp-T_ret)/(Tout_WASTE-T1));
        QHEX1=mdot_1(t).*cpw.*(Tout_WASTE-T1);
        T2(t)=Tout_BTES+Tapp;
        mdot_2(t)=mdot_WASTE-mdot_1(t);

        %HEX2
%         Tin_BTES=Tout_BTES+(mdot_2(t)/(mdot_BTES*Ntot)*(Tout_WASTE-T2(t)));
        Tin_WASTE=(mdot_2(t)*T2(t)+mdot_1(t)*T1)./mdot_WASTE;
    else
        ground_flux=-qg;
        Qborehole=-qg*H;
        deltaT= (ground_flux/ (4 * pi * kg) * integral(intfunc, 1/sqrt(4*ag*time(t)),inf, 'ArrayValued', true));
        Rtot=Rbeff+(deltaT/ground_flux);
        T_gr=8;
        Tb=T_gr+deltaT;    
        Tf= T_gr+Rtot*ground_flux;
        To=1;
        delta_Tev = 5; % K assumed and can be changed
        Tin_BTES=To+delta_Tev;
        mdot_BTES=abs((ground_flux*H)/(2*cpw*(Tin_BTES-Tf)));
        Tout_BTES=Tf-(ground_flux*H)/(2*cpw*mdot_BTES);

        % Defining heat pump input data--> SUPPLY PIPE
        Tin_evap=Tout_BTES;
        Tout_evap=Tin_evap+delta_Tev;
        Tin_BTES=Tout_evap;
        delta_Tcd = 8; % K assumed and can be changed
        Td = 95-delta_Tcd; % condensation temperature
        
        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; % evaporator 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;% Condenser capacity 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; % Compresor power 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; % Condensing temperature function 
        
        % HEX1
        T_supp=40;
        T_ret=25;
        mdot_DH=HL(t)*1000/(cpw*(T_supp-T_ret));
        Tout_WASTE=95;
        Tin_WASTE=T_ret+Tapp;
        mdot_WASTE=result_vector(t).*1000/(cpw*(Tout_WASTE-Tin_WASTE));
        QHEX1=mdot_WASTE*cpw*(Tout_WASTE-Tin_WASTE);
        mdot_1_dis(t)=QHEX1/(cpw*(T_supp-T_ret));
        mdot_HP(t)= mdot_DH-mdot_1_dis(t);

    end
    Tf_array=[Tf_array, Tf];
    Qborehole_array=[Qborehole_array, Qborehole];
    Rtot_array=[Rtot_array, Rtot];
    Tin_BTES_array=[Tin_BTES_array, Tin_BTES];
    Tout_BTES_array=[Tout_BTES_array, Tout_BTES];
    mdot_BTES_array=[mdot_BTES_array, mdot_BTES];
    Tin_evap_array=[Tin_evap_array, Tin_evap];
    Tout_evap_array=[Tout_evap_array, Tout_evap];
    mdot_WASTE_array=[mdot_WASTE_array, mdot_WASTE];
    QHEX1_array=[QHEX1_array, QHEX1];
    Tin_WASTE_array=[Tin_WASTE_array, Tin_WASTE];
    Tout_WASTE_array=[Tout_WASTE_array, Tout_WASTE];
    deltaT_array=[deltaT_array, deltaT];
    T_gr_array=[T_gr_array, T_gr];
    Tb_array=[Tb_array, Tb];
    T_supp_array=[T_supp_array, T_supp];
    T_ret_array=[T_ret_array, T_ret];
    mdot_DH_array=[mdot_DH_array, mdot_DH];
end

figure(1)
plot(Tin_BTES_array)
hold on
plot(Tout_BTES_array)
hold on
plot(Tb_array)

