function [DMA_new,Shift_factor,datashiftedall,temp_structs]=DMA_processing(DMA,T_ref,filename, filepath)

   if isempty(DMA)
       errordlg('The input DMA data is empty.', 'Error', 'modal');
    return;
   end
   
    % Process and classify the DMA data
    T_ref=T_ref+273;
    data = dataprocess(DMA);
    [unique_temps, temp_structs] = dataclassification(data);
    orgplot(unique_temps, temp_structs);
    [dataref, datashiftedall] = datainitilisation(temp_structs, T_ref);

   if isempty(dataref)
       errorMessage = sprintf('No data found for the reference temperature %d°C.', T_ref - 273);
       errordlg(errorMessage, 'Error', 'modal');
     return;
   end

   fprintf('\n*** Start of Master Curve Generation ***\n\n');
    % Separate temperatures lower and higher than T_ref
    temp_lower = unique_temps(unique_temps < T_ref);
    temp_higher = unique_temps(unique_temps > T_ref); 

    % Initialize SF array
    Shift_factor = [T_ref, 0];

    % Process temperatures lower and higher than T_ref
    [Shift_factor, datashiftedall] = lowshift(temp_lower, Shift_factor, temp_structs, dataref, datashiftedall);
    [Shift_factor, datashiftedall] = highshift(temp_higher, Shift_factor, temp_structs, dataref, datashiftedall);
    SF_1=Shift_factor;
    Shift_factor = sortrows(SF_1, 1, 'descend');
    
    % Apply LOESS data smoothing
    [unique_f, smoothed_storage] = smoothData(datashiftedall);

    % Plot shifted and smoothed data
    plotshifted(datashiftedall);
    plotsmoothed(unique_f, smoothed_storage);

    % Save the shift factor and processed DMA data
    [~, name, ~] = fileparts(filename);
    ext01 ='.csv';
    SF_name = fullfile(filepath, ['Shift_Factor_' name ext01]);
    SF_table = array2table(Shift_factor,'VariableNames', {'Temperature', 'Alpha(Shift factor)'});
    writetable(SF_table, SF_name);

    DMA_new = [unique_f', smoothed_storage']; 
    DMA_name = fullfile(filepath, ['Master_Curve_' name ext01]);
    MS_table = array2table(DMA_new,'VariableNames', {'Angular frequency', 'Storage Modulus'});
    writetable(MS_table, DMA_name);

    fprintf('\n*** Master Curve Generation Completed ***\n\n');
end

function data = dataprocess(DMA)
    j = 1;
    data = struct('angf', {},'logangfreq', {}, 'temp', {}, 'storage', {});
    for i = 1:length(DMA)
        f = DMA(i, 1);  % Frequency
        angf = 2 * pi * f;%Transform to angular frequency 
        temp = DMA(i, 2)+273;  % Temperature in K
        storage = DMA(i, 3);  % Storage modulus

       if storage < 0 || isnan(storage)
            continue;
       end
        data(j).angf = angf;
        data(j).logangfreq = log10(angf);
        data(j).temp = temp;
        data(j).storage = storage;
        j = j + 1;
    end

   
    [~, idx] = sort([data.temp]);
    data = data(idx);
end

function [unique_temps, temp_structs] = dataclassification(data)
    unique_temps = unique([data.temp]);
    temp_structs = struct();

    for k = 1:length(unique_temps)
        temp = unique_temps(k);
        idx = [data.temp] == temp;
        temp_data = data(idx);
        temp_struct_name = sprintf('temp_%d', temp);
        temp_structs.(temp_struct_name) = temp_data;
    end
end


function [dataref, datashiftedall] = datainitilisation(temp_structs, T_ref)
    temp_struct_name = sprintf('temp_%d', T_ref);
    if isfield(temp_structs, temp_struct_name)
        dataref = temp_structs.(temp_struct_name); 
        datashiftedall = dataref; 
    else
        dataref = [];
        datashiftedall = [];
    end
end 

function [bestshifts, datashiftedall] = lowshift(temp_lower, bestshifts, temp_structs, datashifted, datashiftedall)
    for i = 1:length(temp_lower)
        temp_shift_name = sprintf('temp_%d', temp_lower(end + 1 - i));
        datatoshift = temp_structs.(temp_shift_name);
        position=2;
        bestshift = calculateShift(datashifted, datatoshift,bestshifts,position);
        
%         if bestshift < bestshifts(end, 2)
%             disp('Negative shift from temperature below reference temperature occurred, skipping iteration');
%             continue;
%         end
        
        bestshifts = [bestshifts; temp_lower(end + 1 - i), bestshift];  
        [datashifted, datashiftedall] = datashift(datatoshift, datashiftedall, bestshift);
    end
end

function [bestshifts, datashiftedall] = highshift(temp_higher, bestshifts, temp_structs, datashifted, datashiftedall)
    for i = 1:length(temp_higher)
        temp_shift_name = sprintf('temp_%d', temp_higher(i));
        datatoshift = temp_structs.(temp_shift_name);
        position=1;
        bestshift = calculateShift(datashifted, datatoshift,bestshifts,position);

%         if bestshift > bestshifts(end, 2)
%             disp('Positive shift from higher temperature occurred, skipping iteration');
%             continue;
%         end
        
        bestshifts = [bestshifts; temp_higher(i), bestshift];  
        [datashifted, datashiftedall] = datashift(datatoshift, datashiftedall, bestshift);
    end
end

function shift = calculateShift(dataRef, dataShift,bestshifts,position)
    baselineX = [dataRef.logangfreq];
    baselineY = [dataRef.storage];
    curveX = [dataShift.logangfreq];
    curveY = [dataShift.storage];
    lastshift=bestshifts(end,2);
    baseTemp=[dataRef(1).temp]-273;
    curveTemp=[dataShift(1).temp]-273;
    

    overlapStart = max(min(baselineY), min(curveY));
    overlapEnd = min(max(baselineY), max(curveY));

    if overlapStart >= overlapEnd
        newrefY=(baselineY+curveY)/2;
        newrefX=baselineX;
        shift1=shiftfunc(baselineX,baselineY,newrefX,newrefY, position);
        shift2=shiftfunc(newrefX,newrefY, curveX,curveY, position);
        shift=shift1+shift2;
        message = sprintf('No overlap, unreliable data possibly produced at %d and %d\n', baseTemp, curveTemp);
        fprintf(message);
        
    else
        baselineIdx = (baselineY >= overlapStart) & (baselineY <= overlapEnd);
        curveIdx = (curveY >= overlapStart) & (curveY <= overlapEnd);

        overlapBaselineX = baselineX(baselineIdx);
        overlapBaselineY = baselineY(baselineIdx);
        overlapCurveX = curveX(curveIdx);
        overlapCurveY = curveY(curveIdx);

        if length(overlapBaselineX) == 1
            if baselineIdx(end) == 1
                baselineIdx(end-1) = 1;
            end
            if baselineIdx(1) == 1
                baselineIdx(2) = 1;
            end
            overlapBaselineX = baselineX(baselineIdx);
            overlapBaselineY = baselineY(baselineIdx);
        end

        if length(overlapCurveX) == 1
            if curveIdx(end) == 1
                curveIdx(end-1) = 1;
            end
            if curveIdx(1) == 1
                curveIdx(2) = 1;
            end
            overlapCurveX = curveX(curveIdx);
            overlapCurveY = curveY(curveIdx);
        end
    costFunction = @(s) sum(((interp1(overlapCurveY,overlapCurveX + s,  overlapBaselineY, 'linear', 'extrap') - overlapBaselineX) .^ 2) ./ (overlapBaselineX .^ 2));

    if position==2
        shift = fminbnd(costFunction, 0,100);%lastshift
    elseif position==1
        shift = fminbnd(costFunction, -100,0); %This is to make sure good points
    end
    end
end

function shift=shiftfunc(Xref,Yref,Xshift,Yshift,position)
    overlapStart = max(min(Yref), min(Yshift));
    overlapEnd = min(max(Yref), max(Yshift));
    baselineIdx = (Yref >= overlapStart) & (Yref <= overlapEnd);
    curveIdx = (Yshift>= overlapStart) & (Yshift<= overlapEnd);
        overlapBaselineX = Xref(baselineIdx);
        overlapBaselineY = Yref(baselineIdx);
        overlapCurveX = Xshift(curveIdx);
        overlapCurveY = Yshift(curveIdx);
 costFunction = @(s) sum(((interp1(overlapCurveY, overlapCurveX + s,  overlapBaselineY, 'linear', 'extrap') - overlapBaselineX) .^ 2) ./ (overlapBaselineX .^ 2));
    if position==2
        shift = fminbnd(costFunction, 0, 100);
    elseif position==1
        shift = fminbnd(costFunction, -100,0);
    end
end

function [datashifted, datashiftedall] = datashift(datatoshift, datashiftedall, bestshift)
    for j = 1:length(datatoshift)
        datatoshift(j).logangfreq = datatoshift(j).logangfreq + bestshift;
    end
    datashifted = datatoshift;
    datashiftedall(end+1:end+length(datatoshift)) = datatoshift;  
    [~, idx] = sort([datashiftedall.logangfreq]);
    datashiftedall = datashiftedall(idx);
end

function [f, storage] = smoothData(datashiftedall)
    [unique_f, ia] = unique([datashiftedall.logangfreq]);
    unique_storage = [datashiftedall(ia).storage];
    f=unique_f;
    span = 0.1;
    storage = smooth(unique_f, unique_storage, span, 'loess');
    storage=storage';
end

function orgplot(unique_temps, temp_structs)
    figure;
    hold on;
    legendEntries = {};
    scatterHandles = [];

    for k = 1:length(unique_temps)
        temp = unique_temps(k);
        temp_struct_name = sprintf('temp_%d', temp);
        finstan = [temp_structs.(temp_struct_name).angf];
        storinstan = [temp_structs.(temp_struct_name).storage];

        hScatter = scatter(finstan, storinstan, 'filled');
        scatterHandles(end+1) = hScatter;
        hLine = plot(finstan, storinstan, 'LineWidth', 2);
        hLine.Color = hScatter.CData;
        temp=temp-273; 
        legendEntries{end+1} = sprintf('%d^oC', temp);
    end

    legend(scatterHandles, legendEntries);
    xlabel('[Angular frequency (rad/s)]');
    ylabel('Storage modulus (MPa)');
    set(gca, 'XScale', 'log');
    hold off;
end


function plotshifted(data)

    figure;
    uniqueTemps = unique([data.temp]);
    colors = lines(length(uniqueTemps)); 
    
    hold on;
    
    
    for i = 1:length(uniqueTemps)
        idx = [data.temp] == uniqueTemps(i);
        fValues = [data(idx).logangfreq];
        storageValues = [data(idx).storage];
        scatter(fValues, storageValues, [], colors(i, :), 'DisplayName', [num2str(uniqueTemps(i)-273),'^oC']);
    end
    
    xlabel('Log10 [Angular frequency (rad/s)]');
    ylabel('Storage Modulus (MPa)');
    
    legend show; 
    hold off;
end

function plotsmoothed(unique_f, smoothed_storage)
    figure;
    plot(10.^unique_f, smoothed_storage);
    set(gca, 'XScale', 'log');
    xlabel('Log10 [Angular frequency (rad/s)]');
    ylabel('Storage Modulus (MPa)');
    title('Plot of Smoothed Data');
    grid on;
end


