%% SUMMARY
% This MATLAB script processes and visualizes Wetness Perception (WP) and Skin Temperature data.
%
% Workflow:
% 1. Load raw data from .mat files for Wetness Perception and Skin Temperature.
% 2. Organize data by participant, location, liquid type, and trial repeats.
% 3. Normalize timestamps and ensure consistent duration (120 seconds).
% 4. Downsample data to 1 sample per second.
% 5. Combine underarm locations (UA1 and UA2) into UA and compute mean & standard deviation.
% 6. Export processed data to Excel files for each location and liquid.
% 7. Generate error bar plots for Wetness Perception and Skin Temperature.
% 8. Create combined figures with dual y-axes for WP and Skin Temperature.
%
% Notes:
% - Ensure .mat files (WetnessPerception.mat and SkinTemperature.mat) are available in the working directory.
% - Figures are saved as TIFF files using the custom function printTiff.
% - Code includes detailed inline comments for clarity.

%% Code to plot figures of Wetness Perception (WP) and Skin Temperature data
clc
clear all
close all

%% Load .mat file with raw data
load('WetnessPerception.mat');   % Load wetness perception data
% load('SkinTemperature.mat');   % Uncomment if loading skin temperature separately

%% Define labels for liquids, repeats, and locations
Liquids = {'Water', 'PG', 'MO', 'PMX'};      % Types of liquids tested
Repeats = {'Repeat1', 'Repeat2', 'Repeat3', 'Repeat4', 'Repeat5'}; % Trial repeats
Locations = {'UA1', 'UA2', 'FA'};            % Test locations: Underarm1, Underarm2, Forearm

% Count participants and other dimensions
nparticipants = numel(fieldnames(Mat));      % Number of participants from data structure
nlocations = length(Locations);
nliquids = length(Liquids);

%% Create participant labels (P01, P02, ...)
part = [];
for np = 1:nparticipants
    if np <= 9
        part{np,1} = sprintf('P0%d',np);
    else
        part{np,1} = sprintf('P%d',np);
    end
end

%% Analyze Wetness Perception data
for k = 1:nliquids
    L = Liquids{1,k};                        % Current liquid
    for j = 1:nlocations
        Loc = Locations{1,j};                % Current location
        Data.WetnessPerception.(Loc).(L) = [];   % Initialize storage
        headers.WetnessPerception.(Loc).(L) = [];
        
        for i = 1:nparticipants
            P = part{i,1};                   % Current participant
            nrepeats = numel(Repeats);
            RepeatsFinder = fieldnames(Mat.(P).(Loc).(L)); % Find repeats in data
            
            for n = 1:nrepeats
                R = Repeats{1,n};            % Current repeat
                for j = 1:numel(RepeatsFinder)
                    temp = [];
                    
                    if strcmp(RepeatsFinder(j),R) == 1
                        % Add header for this participant-repeat
                        h = cellstr(strcat(P,'_', R));
                        headers.WetnessPerception.(Loc).(L) = [headers.WetnessPerception.(Loc).(L) h];
                        
                        % Extract start and end contact times
                        startcontact = Mat.(P).(Loc).(L).StartContact.(R)(1,2);
                        endcontact = Mat.(P).(Loc).(L).EndContact.(R)(1,2);
                        
                        % Check if baseline exists
                        if any(strcmp(fieldnames(Mat.(P).(Loc).(L).ToBaseline),R) == 1) == 1
                            tobaseline = Mat.(P).(Loc).(L).ToBaseline.(R)(1,2);
                            contact = Mat.(P).(Loc).(L).(R)(startcontact:endcontact-1,:);
                            aftercontact = Mat.(P).(Loc).(L).(R)(endcontact:tobaseline-1,:);
                        else
                            contact = Mat.(P).(Loc).(L).(R)(startcontact:endcontact-1,:);
                            aftercontact = [];
                        end
                        
                        %% Ensure 60 seconds after contact
                        if ~isempty(aftercontact)
                            count = 1;
                            l = length(aftercontact);
                            if length(aftercontact) < 600
                                % Pad with zeros if shorter than 60s
                                for loop = 1:10:(600-length(aftercontact))
                                    aftercontact = [aftercontact; (aftercontact(l)+count)*ones(10,1) zeros((10),1)];
                                    count = count + 1;
                                end
                            else
                                % Trim if longer than 60s
                                aftercontact = aftercontact(1:600,:);
                            end
                        else
                            % If no aftercontact, create dummy zeros
                            count = 1;
                            for loop = 1:10:600
                                aftercontact = [aftercontact; (contact(end,1)+count)*ones(10,1) zeros((10),1)];
                                count = count + 1;
                            end
                        end
                        
                        % Combine timestamps and normalize
                        timestamp = [contact(:,1); aftercontact(:,1)];
                        timestamp(:,1) = timestamp(:,1) - timestamp(1,1);
                        
                        % Combine wetness perception values
                        temp = [temp contact(:,2); aftercontact(:,2)];
                        
                        %% Downsample to 1 sample per second
                        [C,ia,ic] = unique(timestamp);
                        for z = 1:length(ia)
                            temp_everysec(z,:) = (temp(ia(z)));
                        end
                        
                        % Store first 120 seconds of data
                        Data.WetnessPerception.(Loc).(L) = [Data.WetnessPerception.(Loc).(L) temp_everysec(1:120,:)];
                    end
                    
                    % Reset temporary arrays
                    aftercontact = [];
                    contact = [];
                    temp_everysec = [];
                end
            end
        end
    end
end

%% Skin Temperature Processing
load('SkinTemperature.mat');  % Load skin temperature measurements

for k = 1:nliquids
    L = Liquids{1,k};                         % Current liquid
    for j = 1:nlocations
        Loc = Locations{1,j};                 % Current location
        Data.SkinTemperature.(Loc).(L) = [];  % Initialize storage
        headers.SkinTemperature.(Loc).(L) = [];
        
        for i = 1:nparticipants
            P = part{i,1};                    % Current participant
            nrepeats = numel(Repeats);
            RepeatsFinder = fieldnames(SkinTemp.(P).(Loc).(L)); % Find repeats
            
            for n = 1:nrepeats
                R = Repeats{1,n};             % Current repeat
                for j = 1:numel(RepeatsFinder)
                    temp = [];
                    
                    if strcmp(RepeatsFinder(j),R) == 1
                        % Add header for this participant-repeat
                        h = cellstr(strcat(P,'_', R));
                        headers.SkinTemperature.(Loc).(L) = [headers.SkinTemperature.(Loc).(L) h];
                        
                        % Extract start and end contact times
                        startcontact = Mat.(P).(Loc).(L).StartContact.(R)(1,2);
                        endcontact = Mat.(P).(Loc).(L).EndContact.(R)(1,2);
                        
                        % Check if baseline exists
                        if any(strcmp(fieldnames(Mat.(P).(Loc).(L).ToBaseline),R) == 1) == 1
                            tobaseline = Mat.(P).(Loc).(L).ToBaseline.(R)(1,2);
                            contact = SkinTemp.(P).(Loc).(L).(R)(startcontact:endcontact-1,:);
                            aftercontact = SkinTemp.(P).(Loc).(L).(R)(endcontact:tobaseline-1,:);
                        else
                            contact = SkinTemp.(P).(Loc).(L).(R)(startcontact:endcontact-1,:);
                            aftercontact = [];
                        end
                        
                        %% Ensure 60 seconds after contact
                        if ~isempty(aftercontact)
                            count = 1;
                            l = length(aftercontact);
                            if length(aftercontact) < 600
                                % Pad with zeros if shorter than 60s
                                for loop = 1:10:(600-length(aftercontact))
                                    aftercontact = [aftercontact; (aftercontact(l)+count)*ones(10,1) zeros((10),1)];
                                    count = count + 1;
                                end
                            else
                                % Trim if longer than 60s
                                aftercontact = aftercontact(1:600,:);
                            end
                        else
                            % If no aftercontact, create dummy zeros
                            count = 1;
                            for loop = 1:10:600
                                aftercontact = [aftercontact; (contact(end,1)+count)*ones(10,1) zeros((10),1)];
                                count = count + 1;
                            end
                        end
                        
                        % Combine timestamps and normalize
                        timestamp = [contact(:,1); aftercontact(:,1)];
                        timestamp(:,1) = timestamp(:,1) - timestamp(1,1);
                        
                        % Combine skin temperature values
                        temp = [temp contact(:,2); aftercontact(:,2)];
                        
                        %% Downsample to 1 sample per second
                        [C,ia,ic] = unique(timestamp);
                        for z = 1:length(ia)
                            temp_everysec(z,:) = (temp(ia(z)));
                        end
                        
                        % Store first 120 seconds of data
                        Data.SkinTemperature.(Loc).(L) = [Data.SkinTemperature.(Loc).(L) temp_everysec(1:120,:)];
                    end
                    
                    % Reset temporary arrays
                    aftercontact = [];
                    contact = [];
                    temp_everysec = [];
                end
            end
        end
    end
end

%% Combine UA1 and UA2 into UA, compute mean and standard deviation
Locations = {'UA', 'FA'};                     % New combined locations
datatype = {'WetnessPerception', 'SkinTemperature'};

for d = 1:length(datatype)
    data = str2mat(datatype(d));
    for k = 1:length(Liquids)
        L = Liquids{1,k};
        
        % Merge UA1 and UA2
        Data.(data).UA.(L) = [Data.(data).UA1.(L) Data.(data).UA2.(L)];
        [r, c] = size(Data.(data).UA.(L));
        
        % Update headers
        headers.(data).UA.(L) = [headers.(data).UA1.(L) headers.(data).UA2.(L)];
        headers.(data).UA.(L) = [headers.(data).UA.(L) 'Mean' 'Standard Deviation'];
        headers.(data).FA.(L) = [headers.(data).FA.(L) 'Mean' 'Standard Deviation'];
        
        % Compute mean and std for UA
        for i = 1:r
            if sum(Data.(data).UA.(L)(i,:)) == 0
                Mean.(data).UA.(L)(i,1) = mean(Data.(data).UA.(L)(i,:));
                StandardDev.(data).UA.(L)(i,1) = std(Data.(data).UA.(L)(i,:));
            else
                removezeros = find(Data.(data).UA.(L)(i,:) == 0);
                temp = Data.(data).UA.(L)(i,:);
                temp(:,removezeros) = [];
                Mean.(data).UA.(L)(i,1) = mean(temp(:,:));
                StandardDev.(data).UA.(L)(i,1) = std(temp(:,:));
            end
        end
        
        % Compute mean and std for FA
        for i = 1:r
            if sum(Data.(data).FA.(L)(i,:)) == 0
                Mean.(data).FA.(L)(i,1) = mean(Data.(data).FA.(L)(i,:));
                StandardDev.(data).FA.(L)(i,1) = std(Data.(data).FA.(L)(i,:));
            else
                removezeros = find(Data.(data).FA.(L)(i,:) == 0);
                temp = Data.(data).FA.(L)(i,:);
                temp(:,removezeros) = [];
                Mean.(data).FA.(L)(i,1) = mean(temp(:,:));
                StandardDev.(data).FA.(L)(i,1) = std(temp(:,:));
            end
        end
    end
end

%% Save processed data to Excel and plot figures
colours = {[0 0.4470 0.7410], [0.8500 0.3250 0.0980], [0.4660 0.6740 0.1880], [0,0,0]+0.6}; % Color scheme

for d = 1:length(datatype)
    data = str2mat(datatype(d));
    for l = 1:length(Locations)
        loc = Locations{l};
        file2save = strcat(loc, '.xls');
        
        for k = 1:length(Liquids)
            L = Liquids{1,k};
            
            % Prepare output for Excel
            output = num2cell([Data.(data).(loc).(L) Mean.(data).(loc).(L) StandardDev.(data).(loc).(L)]);
            out = [headers.(data).(loc).(L); output];
            sheet = strcat(data, '-', L);
            writecell(out, file2save, 'Sheet', sheet);
            
            % Plot error bars
            figure
            errorbar(Mean.(data).(loc).(L), StandardDev.(data).(loc).(L), 'color', colours{k});
            xlabel('Time [s]');
            ylabel(data);
            if d == 2
                axis([0 120 25 35]);  % Skin Temperature range
            else
                axis([0 120 0 100]);  % Wetness Perception range
            end
            figurename = strcat(data, '-', loc, '-', L);
            title(figurename);
            printTiff(gcf, figurename, [15 15], 500); % Save figure as TIFF
        end
    end
end

%% Plot combined figures with dual y-axis
for l = 1:length(Locations)
    for k = 1:length(Liquids)
        L = Liquids{1,k};
        loc = Locations{l};
        
        fig = figure;
        set(fig,'defaultAxesColorOrder',[colours{k}; colours{k}]);
        hold on
        
        yyaxis left
        plot(Mean.WetnessPerception.(loc).(L),'--', 'color', colours{k});
        xlabel('Time [s]');
        ylabel('Wetness Perception');
        axis([0 120 0 100]);
        
        yyaxis right
        plot(Mean.SkinTemperature.(loc).(L), ':', 'LineWidth',1.5, 'color', colours{k});
        ylabel('Skin Temperature');
        axis([0 120 25 35]);
        
        figname = strcat(loc, '-', L);
        title(figname);
        legend('Wetness Perception', 'Skin Temperature');
        printTiff(gcf, figname, [15 15], 500);
    end
end