% Data processing of 2021.07
close all
clc
clear
warning off
format long

%

% Input time frame
start_test = '2021/07/17  03:00:00'; % start of the test
start_heating = '2021/07/17  08:00:00'; % start of the heating
end_heating = '2021/07/18  20:00:00'; % end of the heating
end_test = '2021/07/19  10:00:00'; % end of the test

%

% Read XML data
[fname1,pathn]=uigetfile('*.xls;*.xlsx','Select the files');
[num1,txt1,raw1]=xlsread(fname1,1);
[num2,txt2,raw2]=xlsread(fname1,2);
data_time_Ultima = [txt1 txt2];
data_temperature_Ultima = [num1(:,2:end) num2(:,2:end)];
data_depth_Ultima = num1(:,1);

%

% Time processing

% the test
p = 0;
for i = 1:length(data_time_Ultima)
    t_processing = data_time_Ultima(1,i);
    t_processing = char(t_processing);%
    t_processing = datetime(t_processing,'InputFormat','yyyy/MM/dd HH:mm:ss');
    
    if t_processing>=start_test && t_processing<=end_test
        p = p+1;
        txt_processed(1,p) = data_time_Ultima(1,i);
        temperature_Ultima(:,p) = data_temperature_Ultima(1:end,i);
        t_processed = txt_processed(1,p);
        t_processed = char(t_processed);
        t_processed = datetime(t_processing,'InputFormat','yyyy/MM/dd HH:mm:ss');
        t_seconds_Ultima(1,p) = seconds(t_processed-start_test); 
    end

end

clear t_processing
clear t_processed
clear txt_processed

% the heating stage
p = 0;
for i = 1:length(data_time_Ultima)
    t_processing = data_time_Ultima(1,i);
    t_processing = char(t_processing);
    t_processing = datetime(t_processing,'InputFormat','yyyy/MM/dd HH:mm:ss');
    
    if t_processing>=start_heating && t_processing<=end_heating
        p = p+1;
        txt_processed(1,p) = data_time_Ultima(1,i);
        temperature_Ultima_heating(:,p) = data_temperature_Ultima(1:end,i);
        t_processed = txt_processed(1,p);
        t_processed = char(t_processed);
        t_processed = datetime(t_processing,'InputFormat','yyyy/MM/dd HH:mm:ss');
        t_seconds_Ultima_heating(1,p) = seconds(t_processed-start_heating); 
    end

end

clear t_processing
clear t_processed
clear txt_processed

length = data_depth_Ultima;

%

% Save the origin data
save insitu_measurements.mat

%

% Temperature-length figures
%   Borehole location (32.078-190.19,999-1621,bottom 111.388-1311)
figure(1)
plot(length,temperature_Ultima(:,1:200:end))

%

% View the thermal response test
figure(2)
length_borehole = length(999:1621);
temperature_Ultima_brehole= temperature_Ultima(999:1621,:);
temperature_rise_Ultima_brehole = temperature_Ultima(999:1621,:)-temperature_Ultima(999:1621,1);
clims = [0 8];
imagesc(t_seconds_Ultima,length_borehole,temperature_rise_Ultima_brehole,clims)
xlabel('time(s)');
ylabel('depth(m)');
colorbar

%

% Divide temperature at different depths

% In-stage
q = 0;
for k = 1:size(length)
    
    if length(k,1)>=32.078 && length(k,1)<=111.388
        q = q+1;
        depth_in(q,1) = length(k,1)-32.078;
        temperature_Ultima_borehole_in(q,:) = temperature_Ultima(k,:);
        temperature_Ultima_heating_borehole_in(q,:) = temperature_Ultima_heating(k,:);
    end

end

% Out-stage
q = 0;
for k = 1:size(length)

    if length(k,1)>=111.388 && length(k,1)<=190.19
        q = q+1;
        depth_out(q,1) = 190.19-length(k,1);
        temperature_Ultima_borehole_out(q,:) = temperature_Ultima(k,:);
        temperature_Ultima_heating_borehole_out(q,:) = temperature_Ultima_heating(k,:);
    end

end
depth_out = flip(depth_out);
temperature_Ultima_borehole_out = flip(temperature_Ultima_borehole_out,1);
temperature_Ultima_heating_borehole_out = flip(temperature_Ultima_heating_borehole_out,1);

%

% Initial temperature-depth profile
figure(3)
plot(temperature_Ultima_borehole_in(:,1),depth_in)
set(gca,'YDir','reverse');
ylim([-5 80])
xlim([15 25])
ylabel('Depth (m)');
xlabel('Temperature (^oC)');
grid on
set(gcf,'unit','centimeters','position',[5 5 5 12]);

%

% The temperature response in one depth
figure(4)
plot(t_seconds_Ultima,temperature_Ultima_borehole_in(5,:))
hold on
plot([18000 18000], get(gca, 'YLim'), '-o', 'LineWidth', 1)
hold on
plot([147600 147600], get(gca, 'YLim'), '-o', 'LineWidth', 1)
hold off
legend('1m','Heating-start','Heating-end','Location','Southeast')
xlabel('time(s)');
ylabel('temperature(^oC)');
grid on

%

% Check and calibrate outliers(Optional) 
%   Choose whether to calibrate or not according to the results of the inspection

% Calculate the difference - point by point
for i = 1:size(temperature_Ultima_brehole_heating,2)-1
    difference1_temperature_Ultima_borehole_heating(:,i) = temperature_Ultima_brehole_heating(:,i+1)-temperature_Ultima_brehole_heating(:,i);
end

% Calculate the difference -5min is a calculation cycle
for i = 1:size(temperature_Ultima_brehole_heating,2)-30
    difference2_temperature_Ultima_borehole_heating(:,i) = temperature_Ultima_brehole_heating(:,i+30)-temperature_Ultima_brehole_heating(:,i);
end

% Determine the location
q = 0;
for i= 1:size(difference2_temperature_Ultima_borehole_heating,2)
    if difference2_temperature_Ultima_borehole_heating(end,i)>=0.4 | difference2_temperature_Ultima_borehole_heating(end,i)<=-0.4
        q = q+1;
        discontinuity_points_difference(q) = difference2_temperature_Ultima_borehole_heating(end,i);
        discontinuity_points_location(q) = t_seconds_Ultima_heating(i);
    end
end

% Filter heating early data
for i= 1:size(discontinuity_points_location,2)
    if discontinuity_points_location(1,i)<=20000
        discontinuity_points_difference(1,i) = nan;
        discontinuity_points_location(1,i) = nan;
    end
end

% Determine the compensated difference and delete the data at the original mutation junction
q = 0;
for i= 1:size(temperature_Ultima_brehole_heating,2)
    if t_seconds_Ultima_heating(1,i)<24001-30
        q = q+1;
        temperature_Ultima_brehole_heating_part1(:,q) = temperature_Ultima_brehole_heating(:,i);
        t_seconds_Ultima_heating_part1(1,q) = t_seconds_Ultima_heating(1,i);
    end
end
for j = 1:size(temperature_Ultima_brehole_heating,1)
    compensate_part1_r(j) = mean(temperature_Ultima_brehole_heating_part1(j,end-30:end));
end

q = 0;
for i= 1:size(temperature_Ultima_brehole_heating,2)
    if t_seconds_Ultima_heating(1,i)>24652+30 && t_seconds_Ultima_heating(1,i)<96566-30
        q = q+1;
        temperature_Ultima_brehole_heating_part2(:,q) = temperature_Ultima_brehole_heating(:,i);
        t_seconds_Ultima_heating_part2(1,q) = t_seconds_Ultima_heating(1,i);
    end
end
for j = 1:size(temperature_Ultima_brehole_heating,1)
    compensate_part2_l(j) = mean(temperature_Ultima_brehole_heating_part2(j,1:30));
    compensate_part2_r(j) = mean(temperature_Ultima_brehole_heating_part2(j,end-30:end));
end

q = 0;
for i= 1:size(temperature_Ultima_brehole_heating,2)
    if t_seconds_Ultima_heating(1,i)>96810+30
        q = q+1;
        temperature_Ultima_brehole_heating_part3(:,q) = temperature_Ultima_brehole_heating(:,i);
        t_seconds_Ultima_heating_part3(1,q) = t_seconds_Ultima_heating(1,i);
    end
end
for j = 1:size(temperature_Ultima_brehole_heating,1)
    compensate_part3_r(j) = mean(temperature_Ultima_brehole_heating_part3(j,1:30));
end

t_seconds_Ultima_heating_correction = [t_seconds_Ultima_heating_part1 t_seconds_Ultima_heating_part2 t_seconds_Ultima_heating_part3];
for j = 1:size(temperature_Ultima_brehole_heating,1)
    temperature_Ultima_brehole_heating_part2(j,:) = temperature_Ultima_brehole_heating_part2(j,:)+compensate_part1_r(j)-compensate_part2_l(j);
    temperature_Ultima_brehole_heating_part3(j,:) = temperature_Ultima_brehole_heating_part3(j,:)+compensate_part1_r(j)-compensate_part2_l(j);
    temperature_Ultima_brehole_heating_part3(j,:) = temperature_Ultima_brehole_heating_part3(j,:)-compensate_part3_r(j)+compensate_part2_r(j);
    temperature_Ultima_brehole_heating_correction = [temperature_Ultima_brehole_heating_part1 temperature_Ultima_brehole_heating_part2 temperature_Ultima_brehole_heating_part3];
end

%

% Data thinning(Optional) 
t_seconds_Ultima_heating_correction_thin = downsample(t_seconds_Ultima_heating_correction,6); % Imput thinning parameters
for j = 1:size(temperature_Ultima_brehole_heating_correction,1)
    temperature_Ultima_brehole_heating_correction_thin(j,:) = downsample(temperature_Ultima_brehole_heating_correction(j,:),6);
end

%

% Save the data for flow rate estimation
save estimation_data.mat