close all;
clear all;
clc;
%%
%DATEIIMPORT
% Dialog für die Auswahl der FTIR Datei anzeigen
[Dateiname, Dateipfad] = uigetfile('*.txt', 'FTIR-Datei auswählen');
cd (Dateipfad)

% Überprüfen, ob eine Datei ausgewählt wurde
if isequal(Dateiname, 0) || isequal(Dateipfad, 0)
    disp('Auswahl abgebrochen.');
end
FTIR                = importdata(Dateiname);
FTIR_Zeitstempel    = FTIR.textdata;
FTIR_Daten          = FTIR.data;

% Dialog für die Auswahl der Traversen Datei anzeigen
[Dateiname1, Dateipfad1] = uigetfile('*.txt', 'Traversen-Datei auswählen');

% Überprüfen, ob eine Datei ausgewählt wurde
if isequal(Dateiname1, 0) || isequal(Dateipfad1, 0)
    disp('Auswahl abgebrochen.');
end
Traverse                = importdata(Dateiname1);
Traverse_Zeitstempel    = Traverse.textdata;
Traverse_Daten          = Traverse.data;

% Dialog für die Auswahl der Labview Datei anzeigen
[Dateiname2, Dateipfad2] = uigetfile('*.txt', 'Labview-Datei auswählen');
cd (Dateipfad)

% Überprüfen, ob eine Datei ausgewählt wurde
if isequal(Dateiname2, 0) || isequal(Dateipfad2, 0)
    disp('Auswahl abgebrochen.');
end
LABVIEW                = importdata(Dateiname2,'\t');
LABVIEW_Zeitstempel    = LABVIEW.textdata;
LABVIEW_Daten          = LABVIEW.data;

%%
%GLOBALE VARIABLEN DEFINIEREN
%Zeiten in [s] eintragen
Messdelay   = 10;
Messdauer   = 20;
%% 
%FTIR UND TRAVERSE DATEN KOMBINIEREN
%Zeitstempel von Pos_erreicht wird aus Traversendatei übergeben 
L = length(Traverse_Daten);
%Abfrage ob Labview-Log nur T enthält
s=size(LABVIEW_Daten);
if s(2)==1
Data = zeros(L,6);
else
Data = zeros(L,9);
end

for i = 1:1:L
Pos_erreicht = char(extractAfter(Traverse_Zeitstempel(i+1,1),11)); %+1 wegen Header

%Zeitstempel in absolute Sekunden umrechnen
Messung_start_abs   = Messdelay+str2double(Pos_erreicht(1:2))*3600+str2double(Pos_erreicht(4:5))*60+str2double(Pos_erreicht(7:8));
Messung_ende_abs    = Messung_start_abs+Messdauer;

%Sekunden in Zeitstempel zurückrechnen // Index _T für Temperaturdaten
Messung_start_h     = num2str(floor(Messung_start_abs/3600),'%02.f');
Messung_start_min   = num2str(floor(mod(Messung_start_abs,3600)/60),'%02.f');
Messung_start_sec   = num2str(floor(mod(mod(Messung_start_abs,3600),60)),'%02.f');
Messung_start       = strcat(Messung_start_h,':',Messung_start_min,':',Messung_start_sec,'.00');
Messung_start_T     = strcat(Messung_start_h,':',Messung_start_min,':',Messung_start_sec);

Messung_ende_h      = num2str(floor(Messung_ende_abs/3600),'%02.f');
Messung_ende_min    = num2str(floor(mod(Messung_ende_abs,3600)/60),'%02.f');
Messung_ende_sec    = num2str(floor(mod(mod(Messung_ende_abs,3600),60)),'%02.f');
Messung_ende        = strcat(Messung_ende_h,':',Messung_ende_min,':',Messung_ende_sec,'.00');
Messung_ende_T      = strcat(Messung_ende_h,':',Messung_ende_min,':',Messung_ende_sec);

%Indizes zu den Zeitstempeln herauslesen
Index_start     = find(contains(FTIR_Zeitstempel,Messung_start));
Index_start_T   = find(contains(LABVIEW_Zeitstempel,Messung_start_T));

Index_ende      = find(contains(FTIR_Zeitstempel,Messung_ende));
Index_ende_T    = find(contains(LABVIEW_Zeitstempel,Messung_ende_T));

%Mittelwert der Messdaten im Zeitintervall Start-Ende
NOx_sum     = 0;
NH3_sum     = 0;
HNCO_sum    = 0;
T_sum       = 0;
BL_sum      = 0;
VDL_sum     = 0;
T_BKaus_sum = 0;

    for k = Index_start-1:1:Index_ende-1      %-1 weil Daten-Matrix keinen Header hat
        NOx_sum     = NOx_sum+FTIR_Daten(k,7);
        NH3_sum     = NH3_sum+FTIR_Daten(k,29);
        HNCO_sum    = HNCO_sum+FTIR_Daten(k,22);
    end
    
NOx_value   = NOx_sum/(Index_ende-Index_start+1);
NH3_value   = NH3_sum/(Index_ende-Index_start+1);
HNCO_value  = HNCO_sum/(Index_ende-Index_start+1);

if s(2)==1
    for k = Index_start_T-1:1:Index_ende_T-1      %-1 weil Daten-Matrix keinen Header hat
        T_sum       = T_sum+LABVIEW_Daten(k,1); % 1 oder 3 wählen, je nach dem ob im Labview-Log mehr Daten aus die Temperatur gespeichert ist
    end
    T_value = T_sum/(Index_ende_T-Index_start_T+1);
else
    for k = Index_start_T-1:1:Index_ende_T-1      %-1 weil Daten-Matrix keinen Header hat
        BL_sum      = BL_sum+LABVIEW_Daten(k,2);
        VDL_sum     = VDL_sum+LABVIEW_Daten(k,1);
        T_sum       = T_sum+LABVIEW_Daten(k,3); 
        T_BKaus_sum = T_BKaus_sum+LABVIEW_Daten(k,6);
    end
    T_value = T_sum/(Index_ende_T-Index_start_T+1);
    BL_value = BL_sum/(Index_ende_T-Index_start_T+1);
    VDL_value = VDL_sum/(Index_ende_T-Index_start_T+1);
    T_BKaus_value = T_BKaus_sum/(Index_ende_T-Index_start_T+1);
end

%Koordinaten mit Werten verknüpfen
     Data(i,1)= Traverse_Daten(i,1); %x-Werte
     Data(i,2)= Traverse_Daten(i,2); %y-Werte
     Data(i,3)= NOx_value;
     Data(i,4)= NH3_value;
     Data(i,5)= HNCO_value;
     Data(i,6)= T_value;
     if s(2)>1
     Data(i,7)= BL_value;
     Data(i,8)= VDL_value;
     Data(i,9)= T_BKaus_value;
     end
     
end

%%
% NOx WERTE IN MATRIX KONVERTIEREN UND PLOTTEN
nx = length(unique(Data(:,1))) ; 
ny = length(unique(Data(:,2))) ; 
x_wertebereich = unique(Data(:,1));
y_wertebereich = flip(unique(Data(:,2)));
NOx_Map     = nan(ny,nx);
NH3_Map     = nan(ny,nx);
HNCO_Map    = nan(ny,nx);
T_Map       = nan(ny,nx);

if s(2)>1
BL_Map      = nan(ny,nx);
VDL_Map     = nan(ny,nx);
T_BKaus_Map = nan(ny,nx);
end

for i = 1:1:nx
    for k = 1:1:ny
        I = find(Data(:,1)==x_wertebereich(i) & Data(:,2)==y_wertebereich(k));
        if isempty(I)==1
   
        else 
        NOx_Map(k,i)    = Data(I,3);
        NH3_Map(k,i)    = Data(I,4);
        HNCO_Map(k,i)   = Data(I,5);
        T_Map(k,i)      = Data(I,6);
            if s(2)>1
                BL_Map(k,i)      = Data(I,7);
                VDL_Map(k,i)     = Data(I,8);
                T_BKaus_Map(k,i) = Data(I,9);
            end
        end
    end
end
%Map-Mittelwerte berechnen 
NOx_avg  = round(mean(NOx_Map,'all','omitnan'),2);
NH3_avg  = round(mean(NH3_Map,'all','omitnan'),2);
HNCO_avg = round(mean(HNCO_Map,'all','omitnan'),2);
T_avg    = round(mean(T_Map,'all','omitnan'),1);

if s(2)>1
    BL_avg      = round(mean(BL_Map,'all','omitnan'),1);
    VDL_avg     = round(mean(VDL_Map,'all','omitnan'),1);
    T_BKaus_avg = round(mean(T_BKaus_Map,'all','omitnan'),1);
end

%Betriebsbedigungen auslesen
velo = extractAfter(extractBefore(Dateiname1,4),1);
Temp = extractAfter(extractBefore(Dateiname1,9),5);
Air  = extractAfter(extractBefore(Dateiname1,13),10);
Urea = strcat(extractAfter(extractBefore(Dateiname1,16),14),',',extractAfter(extractBefore(Dateiname1,17),15));

%%
%GLEICHVERTEILUNGS MAP ERSTELLEN UND INDEX BERECHNEN
%Lokaler Ungleichförmigkeitsindex
UI_value = zeros(L,1);
for i=1:1:L
    UI_value(i)=sqrt((Data(i,3)-NOx_avg)^2)/NOx_avg;
end

%In Matrix konvertieren
UI_map = nan(ny,nx);
for i = 1:1:nx
    for k = 1:1:ny
        I = find(Data(:,1)==x_wertebereich(i) & Data(:,2)==y_wertebereich(k));
        if isempty(I)==1
   
        else
        UI_map(k,i) = UI_value(I);
        end
    end
end

%Globalen Gleichförmigkeits Index berechnen 
UI=1-sum(UI_value)/L;

%%
%2D-Maps erstellen
set(gca, 'FontName', 'Arial','FontSize',13)
figure
hold on
surf(x_wertebereich,y_wertebereich,NOx_Map)
shading interp
view(2)
clim([50, 250])
xlabel('x-Position [mm]')
ylabel('y-Position [mm]')
colorbar
ylabel(colorbar,'NO_{x}-Konzentration [ppm]','FontSize',11)
%title ({strcat('T Gas=',Temp,'°C ',' v Gas=',velo,'m/s');strcat('Air=',Air,'l/min ',' Urea=',Urea,'kg/h')})
txt={'mean=',strcat(num2str(NOx_avg),' [ppm]')};
text(-145,135,txt)
txt={'UI=',strcat(num2str(UI))};
text(-145,-130,txt)
hold off
exportgraphics(gcf, 'NOx.emf', 'ContentType', 'vector');


figure
hold on
surf(x_wertebereich,y_wertebereich,NH3_Map)
shading interp
view(2)
clim([0, 25])
xlabel('x-Position in mm')
ylabel('y-Position in mm')
colorbar
ylabel(colorbar,'NH3 Concentration in ppm','FontSize',11)
title ({strcat('T Gas=',Temp,'°C ',' v Gas=',velo,'m/s');strcat('Air=',Air,'l/min ',' Urea=',Urea,'kg/h')})
txt={'mean=',strcat(num2str(NH3_avg),' ppm')};
text(-145,135,txt)
hold off
exportgraphics(gcf, 'NH3.emf', 'ContentType', 'vector');



figure
hold on
surf(x_wertebereich,y_wertebereich,HNCO_Map)
shading interp
view(2)
clim([0, 8])
xlabel('x-Position in mm')
ylabel('y-Position in mm')
colorbar
ylabel(colorbar,'HNCO Concentration in ppm','FontSize',11)
title ({strcat('T Gas=',Temp,'°C ',' v Gas=',velo,'m/s');strcat('Air=',Air,'l/min ',' Urea=',Urea,'kg/h')})
txt={'mean=',strcat(num2str(HNCO_avg),' ppm')};
text(-145,135,txt)
hold off
exportgraphics(gcf, 'HNCO.emf', 'ContentType', 'vector');


figure
hold on
surf(x_wertebereich,y_wertebereich,T_Map)
shading interp
view(2)
clim([350, 500])
xlabel('x-Position in mm')
ylabel('y-Position in mm')
colorbar
ylabel(colorbar,'T in °C','FontSize',11)
title ({strcat('T Gas=',Temp,'°C ',' v Gas=',velo,'m/s');strcat('Air=',Air,'l/min ',' Urea=',Urea,'kg/h')})
txt={'mean=',strcat(num2str(T_avg),' °C')};
text(-145,135,txt)  
hold off
exportgraphics(gcf, 'T.emf', 'ContentType', 'vector');



%%
%Tabellen zur Datenausgabe
%Export T- und NOx-Map
NOx_Map_export  = nan(ny+1,nx+1);
T_Map_export    = nan(ny+1,nx+1);
for i=1:1:nx
    NOx_Map_export(1,i+1)   = x_wertebereich(i);
    NOx_Map_export(i+1,1)   = y_wertebereich(i);
    T_Map_export(1,i+1)     = x_wertebereich(i);
    T_Map_export(i+1,1)     = y_wertebereich(i);
    for k=1:1:ny
        NOx_Map_export(i+1,k+1) = NOx_Map(i,k);
        T_Map_export(i+1,k+1)   = T_Map(i,k);
    end
end

dlmwrite('NOx_Map.txt',NOx_Map_export,'delimiter', '\t');
dlmwrite('T_Map.txt',T_Map_export,'delimiter', '\t');
%export stats
if s(2)>1
stats=table(VDL_avg,BL_avg,T_avg,T_BKaus_avg,NOx_avg);
writetable(stats,'stats.txt','Delimiter','\t')  
end    
%%
%UI 2D-Map erstellen
figure 
hold on
surf(x_wertebereich,y_wertebereich,UI_map)
shading interp
view(2)
clim([0, 0.4])
colorbar
title (num2str(UI))
hold off