%=======================================================================
% PROCESSREALTIMEPOSITIONS (Octave Version)
% Processes, plots, and exports real-time position data.
%=======================================================================
function [rtpos] = ProcessRealtimePositions(T0xFileNameOrDir, varargin)

% --- OCTAVE PACKAGE LOADING ---
try
    pkg load statistics;
    pkg load nan;
    % pkg load io; % For xlsx support if available
catch
    % warning('Statistics or NaN packages not found. Some calculations may fail.');
end

% --- DEFAULT OPTIONS (Replacing arguments block) ---
options = struct();
options.titel = 'Position Solution';
options.FixOnly = false;
options.AntennaArpHeight = NaN;
options.stime = 1;
options.etime = 10E10;
options.stimeGPS = []; 
options.etimeGPS = [];
options.ylim = [];
options.RefPosXYZ = [];
options.RefPosLLH = [];
options.RefPosNEU = [];
options.showFig = false;
options.scaledFontSize = 8; % Slightly larger for Octave plots
options.fontName = 'Arial';
options.OutputCSV = true;
options.OutputDir = [];
options.LineWidth = 1;
options.CoordFrame = 'Local';
options.CoordUnits = 'Meters';
options.satMapBuffer = 20;

% Reporting
options.ReportGen = true;
options.ReportFormat = 'HTML'; 
options.GapThreshold = 300;
options.scriptVersion = '1.34-Octave';
options.LogFcn = [];
options.RoundingDigits = 4;

% Semi-Kinematic
options.SemiKinematic = false;
options.MinStaticDuration = 5;
options.StaticClusterRadius = 0.03;
options.StaticVelocityThresh = 0.02;
options.StaticConvergenceThresh = 0.003;
options.StaticSettlingThresh = 0.008;
options.StaticSettlingWindow = 2.0;

% Rotation
options.RotationAnalysis = false;
options.RotGapThreshold = 30;
options.RotMinSpeed = 0.02;
options.RotMaxRadius = 50;
options.RotFitThreshold = 0.15;

options.ThumbRes = 100; % Reduced for Octave performance
options.PlotFixStatusColumn = true;
options.ReturnGNSSInfoOnly = false;
options.OverwriteExistingResults = false;
options.metaData = struct();

% --- INPUT PARSING ---
if nargin < 1 || isempty(T0xFileNameOrDir)
    T0xFileNameOrDir = pwd;
end

% Parse Name-Value pairs from varargin
if nargin > 1
    for i = 1:2:length(varargin)
        if (i + 1) <= length(varargin)
            key = varargin{i};
            val = varargin{i+1};
            % Allow dynamic field creation for struct options
            options.(key) = val;
        end
    end
end

% Ensure char types
T0xFileNameOrDir = char(T0xFileNameOrDir);

% --- VERSION CONTROL ---
scriptVersion = '1.34';
options.scriptVersion = scriptVersion;

% Initialize the metaData struct first.
options.metaData = initMetaData(T0xFileNameOrDir, options);

if options.ReturnGNSSInfoOnly
    try
        infoStruct = getFileInfoOnly(T0xFileNameOrDir, options);
        rtpos = infoStruct; % "highjack" the rtpos struct to return the data to fill UI textfields
    catch ME
        warning('GNSS data could not be read: %s', ME.message);
        rtpos = struct();
    end
    return; % Important: End the script here before actual processing starts
end

if strcmpi(T0xFileNameOrDir, '-v') || strcmpi(T0xFileNameOrDir, '--version')
    rtpos = scriptVersion;
    return;
end

% --- 0. INITIALIZE FIX MANAGER ONCE ---
options.FixStatusManager = FixStatusManager();

if options.SemiKinematic
    options.GapThreshold = 1e10;
    logMsg(sprintf('Semi Kinematic processing selected, deactivating gap search in file.\n'), options);
end

% --- 1. HELP COMMAND CHECK ---
if strcmpi(T0xFileNameOrDir, '-h') || strcmpi(T0xFileNameOrDir, '--help') || strcmpi(T0xFileNameOrDir, '/?')
    printHelp();
    rtpos = [];
    return;
end

tTotal = tic;
logMsg(sprintf('Start: ProcessRealtimePositions v%s', scriptVersion), options);
rtpos = {};

% --- UTC to GPS Time Conversion Logic ---
leap_seconds = 0;
sample_rate = 1;
isTimeRangeSet = false;

if ~isempty(options.stimeGPS)
    try
        if ischar(options.stimeGPS) || isstring(options.stimeGPS)
             date_str = datestr(datevec(options.stimeGPS), 'dd.mm.yyyy');
             time_str = datestr(datevec(options.stimeGPS), 'HH:MM:SS');
             [~, tow_vector] = utc_to_gps_time(date_str, time_str, leap_seconds, sample_rate);
             options.stime = tow_vector(1);
             logMsg(sprintf('stimeGPS provided. Overwriting stime to: %.2f', options.stime), options);
             isTimeRangeSet = true;
        end
    catch ME
        warning(ME.identifier, 'Failed to convert stimeGPS. Using default stime. Error: %s', ME.message);
    end
end

if ~isempty(options.etimeGPS)
    try
        if ischar(options.etimeGPS) || isstring(options.etimeGPS)
             date_str = datestr(datevec(options.etimeGPS), 'dd.mm.yyyy');
             time_str = datestr(datevec(options.etimeGPS), 'HH:MM:SS');
             [~, tow_vector] = utc_to_gps_time(date_str, time_str, leap_seconds, sample_rate);
             options.etime = tow_vector(1);
             logMsg(sprintf('etimeGPS provided. Overwriting etime to: %.2f', options.etime), options);
             isTimeRangeSet = true;
        end
    catch ME
        warning(ME.identifier,'Failed to convert etimeGPS. Using default etime. Error: %s', ME.message);
    end
end

if options.stime ~= 1 || options.etime ~= 10E10
    isTimeRangeSet = true;
end

% --- 2. VALIDATION & FILE LISTING ---
htmlScanDir = "";

if isfolder(T0xFileNameOrDir)
    htmlScanDir = T0xFileNameOrDir;
    files = dir(fullfile(T0xFileNameOrDir, '*.T0?'));
    if isempty(files)
        logMsg('No *.t0? files found in directory.', options);
        return;
    end
    % OPTIMIZATION: Preallocate cell array
    nFiles = length(files);
    T0xFilesToProcess = cell(1, nFiles);
    for i = 1:nFiles
        T0xFilesToProcess{i} = fullfile(files(i).folder, files(i).name);
    end
elseif exist(T0xFileNameOrDir, 'file') == 2
    T0xFilesToProcess = {T0xFileNameOrDir};
    [p, ~, ~] = fileparts(T0xFileNameOrDir);
    if isempty(p)
        p = pwd;
    end
    htmlScanDir = p;
else
    fprintf('\nError: Input "%s" not found.\n', T0xFileNameOrDir);
    rtpos = [];
    return;
end

if ~isempty(options.OutputDir)
    % only create if it is not existing yet
    if ~isfolder(options.OutputDir)
        mkdir(options.OutputDir);
    end
    htmlScanDir = options.OutputDir;
end

% --- 3. AUTO GAP DETECTION MODE ---
if ~isTimeRangeSet && ~isempty(T0xFilesToProcess)
    logMsg('Mode: Auto Gap Detection', options);
    allResults = {};
    for fileIndex = 1:length(T0xFilesToProcess)
        T0xFileName = T0xFilesToProcess{fileIndex};
        logMsg(sprintf('Scanning file %d/%d: %s', fileIndex, length(T0xFilesToProcess), T0xFileName), options);
        tScan = tic;

        allGPS_Sec = getFileGPSData(T0xFileName, options);

        if isempty(allGPS_Sec)
            logMsg(['  No GNSS data found for file: ', T0xFileName], options);
            continue;
        end

        segments = splitDataByGap(allGPS_Sec, options.GapThreshold);

        if length(segments) > 1
            logMsg(sprintf('  Gap Detected! Found %d segments (gaps > %.0fs). Scan time: %.2fs', ...
                length(segments), options.GapThreshold, toc(tScan)), options);
        else
            logMsg(sprintf('  Data is continuous (no gaps > %.0fs). Processing as 1 segment. Scan time: %.2fs', ...
                options.GapThreshold, toc(tScan)), options);
        end

        for i = 1:length(segments)
            seg = segments{i};
            segmentOptions = options;
            segmentOptions.stime = seg.startTime;
            segmentOptions.endTime = seg.endTime; 

            if length(segments) > 1
                logMsg(sprintf('  -- Processing Segment %d: %.2f to %.2f', i, seg.startTime, seg.endTime), options);
            end

            res = processFile(T0xFileName, segmentOptions);
            if ~isempty(res), allResults{end+1} = res; end
        end
    end
    rtpos = allResults;

    if options.ReportGen

        if options.SemiKinematic
            sDir = 'SemiKinematic';
        else
            sDir = 'Static';
        end

        targetReportDir = fullfile(char(htmlScanDir), sDir);

        if isfolder(targetReportDir)
            generateReportsRouter(targetReportDir, options);
        end
    end

    logMsg(sprintf('Total Execution Time: %.2fs', toc(tTotal)), options);
    return;
end

% --- 4. STANDARD MODE ---
for i = 1:length(T0xFilesToProcess)
    logMsg(sprintf('Processing File %d/%d', i, length(T0xFilesToProcess)), options);
    res = processFile(T0xFilesToProcess{i}, options);
    if ~isempty(res), rtpos{end+1} = res; end
end

if options.ReportGen
    if options.SemiKinematic, sDir = 'SemiKinematic'; else, sDir = 'Static'; end
    targetReportDir = fullfile(char(htmlScanDir), sDir);

    if isfolder(targetReportDir)
        generateReportsRouter(targetReportDir, options);
    end
end

logMsg(sprintf('Total Execution Time: %.2fs', toc(tTotal)), options);
end



%=======================================================================
%PROCESSFILE Core processing logic for a single T0x file.
%   Decodes, filters, and computes statistics for the given file.
%=======================================================================
function [rtpos] = processFile(T0xFileName, options)

if iscell(T0xFileName)
    T0xFileName = T0xFileName{1};
end
T0xFileName = char(T0xFileName);

% --- LOGGING SETUP ---
options.logState = containers.Map({'Buffer', 'File'}, {{}, ''});

tFile = tic;
[filepath, name, ext] = fileparts(T0xFileName);
DisplayFileName = strcat(name, ext);


if options.SemiKinematic
    subDirName = 'SemiKinematic';
else
    subDirName = 'Static';
end

outputDir = fullfile(filepath, options.OutputDir, subDirName);
% Simple check if OutputDir was absolute. In Octave, simpler to just use fullfile logic
if ~isempty(options.OutputDir)
     if (isunix && options.OutputDir(1) == '/') || (ispc && length(options.OutputDir)>1 && options.OutputDir(2)==':')
         outputDir = fullfile(options.OutputDir, subDirName);
     end
end

if isAbsolute
    % Case A: User provided absolute path (e.g., from App) -> Use directly
    outputDir = fullfile(options.OutputDir, subDirName);
else
    % Case B: Relative path or empty -> Prepend Input File Directory
    outputDir = fullfile(filepath, options.OutputDir, subDirName);
end

% Create folder if not existing
if ~isfolder(outputDir)
    mkdir(outputDir);
elseif options.OverwriteExistingResults
    % delete files in this folder if wanted
    try
        filesToDelete = fullfile(outputDir, '*.*');
        fls = glob(filesToDelete);
        for k=1:length(fls), delete(fls{k}); end
        logMsg(sprintf('    Overwriting enabled: Cleared existing results for "%s".', outputDir), options);
    catch ME
        logMsg(sprintf('    Warning: Could not delete old files: %s', ME.message), options);
    end
end

% --- PATH LENGTH CHECK ---
estimatedMaxLen = length(fullfile(filepath, name)) + 85;
if estimatedMaxLen > 259
    errMsg = sprintf('Skipping File: Path too long! (%d chars). Limit is ~260.\n    -> Move input file to a shorter folder (e.g. C:\\Data)', estimatedMaxLen);
    logMsg(errMsg, options);
    rtpos = [];
    return;
end

if ~exist(T0xFileName, 'file')
    logMsg('  Error: File not found.', options);
    rtpos = [];
    return;
end

utcTitleStr = "";
if ~isempty(options.stimeGPS), utcTitleStr = [utcTitleStr, ' (Start specified)']; end

% --- INITIALIZE METADATA (Early Init for D16 storage) ---
metaData = struct();
metaData.Filename = DisplayFileName;
metaData.ScriptVersion = options.scriptVersion;
metaData.CoordFrame = options.CoordFrame;
% Initialize standard D16 fields to empty in case file is missing
metaData.ReceiverId = '';
metaData.RtkVersion = '';
metaData.IonoMitigationEnabled = '';
metaData.IonoMitigationFallBackmode = '';
metaData.AntennaApcHeight = NaN;
metaData.AntennaType = NaN;
metaData.AntennaArpHeight = options.AntennaArpHeight;
metaData.ThrottleEngineHighLoadCount = 0;

% decoding the d16 records first
[~, name, ~] = fileparts(T0xFileName);
ffpname = fullfile(filepath,[name,'.d16']);
tStep = tic;

if ~exist(ffpname,"file")
    logMsg(sprintf('  > Decoding/Reading Records D16...'), options);

    pathToViewDat = getPathToViewDat(); % Use helper function

    flag = '-d16';
    cmd = sprintf('"%s" %s "%s" > "%s"', pathToViewDat, flag, T0xFileName, ffpname);
    [status, ~] = system(cmd);

    if status ~= 0
        logMsg('    Decoding Record D16 failed. Processing without antenna information.', options);
        DataD16 = [];
    else
            DataD16 = ReadD16(ffpname);
    end

else
    logMsg('    Found existing Record D16 file. Reparsing instead of decoding again.', options);
        DataD16 = ReadD16(ffpname);
end

% --- SAVE D16 FIELDS TO METADATA ---
if ~isempty(DataD16)
    if isfield(DataD16, 'ReceiverId'), metaData.ReceiverId = DataD16.ReceiverId; end
    if isfield(DataD16, 'RtkVersion'), metaData.RtkVersion = DataD16.RtkVersion; end
    if isfield(DataD16, 'IonoMitigationEnabled'), metaData.IonoMitigationEnabled = DataD16.IonoMitigationEnabled; end
    if isfield(DataD16, 'IonoMitigationFallBackmode'), metaData.IonoMitigationFallBackmode = DataD16.IonoMitigationFallBackmode; end
    if isfield(DataD16, 'AntennaApcHeight'), metaData.AntennaApcHeight = DataD16.AntennaApcHeight; end
    if isfield(DataD16, 'AntennaType'), metaData.AntennaType = DataD16.AntennaType; end
    if isfield(DataD16, 'AntennaArpHeight'), metaData.AntennaArpHeight = DataD16.AntennaArpHeight; end
    if isfield(DataD16, 'ThrottleEngineHighLoadCount'), metaData.ThrottleEngineHighLoadCount = DataD16.ThrottleEngineHighLoadCount; end

    % Log found metadata
    logMsg(sprintf('    Found Antenna Info: (Antenna Type: %d, Antenna ARP Height: %.3fm)', DataD16.AntennaType, DataD16.AntennaArpHeight), options);
end

logMsg(sprintf('    Done (%.2fs). Record D16 parsing.', toc(tStep)), options);

% Loading the Antenna.ini
tStep = tic;
logMsg(sprintf('  > Loading antenna.ini'), options);
try
    % Attempt to create the database instance
    db = TrimbleAntennaDB();
    logMsg(sprintf('    Done (%.2fs). Database loaded successfully from: %s', toc(tStep), db.LoadedFilePath),options);

catch ME
    % Handle the error gracefully
    % logMsg('    Warning: TrimbleAntennaDB not found or failed.', options);
    logMsg(sprintf('An unexpected error occurred getting TrimbleAntennaDB(): %s\n', ME.message),options);
end

% --- DECODING POSITIONS ---
tStep = tic;
logMsg(sprintf('  > Decoding/Reading %s...', name), options);
isT04 = contains(upper(ext), 'T04');
cleanExt = ".d35_2";
if ~isT04, cleanExt = ".d29"; end
finalCleanPath = fullfile(filepath, strcat(name, cleanExt));

if ~exist(finalCleanPath, 'file')
    if isT04, flag = "-d35:2"; else, flag = "-d29"; end
    pathToViewDat = getPathToViewDat();
    rawTempPath = fullfile(filepath, strcat(name, '_RAW_TEMP.txt'));
    cmd = sprintf('"%s" %s -h -m "%s" > "%s"', pathToViewDat, flag, T0xFileName, rawTempPath);
    [status, ~] = system(cmd);
    if status ~= 0
        logMsg('    Decoding failed.', options);
        if exist(rawTempPath, 'file'), delete(rawTempPath); end
        rtpos = []; return;
    end
    if exist(rawTempPath, 'file'), movefile(rawTempPath, finalCleanPath); end
end

if ~exist(finalCleanPath, 'file'), rtpos = []; return; end

% 2. Reading
if isT04
    try
        rtpos = import_d35_2(finalCleanPath, options.FixStatusManager);
        % Manually subset by time because import_d35_2 usually reads everything
        if ~isempty(rtpos) && isfield(rtpos, 'GPS_Sec')
            rtpos = subsetStruct(rtpos, rtpos.GPS_Sec >= options.stime & rtpos.GPS_Sec <= options.etime);
        end
    catch ME
        logMsg(sprintf('    Error using import_d35_2: %s', ME.message), options);
        rtpos = []; return;
    end
else
    rtpos = ReadD29(finalCleanPath);
    rtpos = subsetStruct(rtpos, rtpos.GPS_Sec >= options.stime & rtpos.GPS_Sec <= options.etime);
end

if isempty(rtpos) || isempty(rtpos.GPS_Sec)
    logMsg('    No data loaded.', options); rtpos = []; return;
end

% Preprocessing to compute ground point heights in case of GNSS only
rtpos = PreprocessPositionsData(rtpos, metaData, options);

logMsg(sprintf('    Done (%.2fs). Points: %d', toc(tStep), length(rtpos.GPS_Sec)), options);

% 3. Filter Fix Status
if options.FixOnly
    fixStatusesToPlot = options.FixStatusManager.FixStatusesToPlot;
    filteredIndices = ismember(rtpos.FixStatus, fixStatusesToPlot);
    if ~any(filteredIndices)
        logMsg('    No data with specified fix status.', options); rtpos = []; return;
    end
    rtpos = subsetStruct(rtpos, filteredIndices);
end

% 4. Coordinate Calculations
try
    ref = [median(deg2rad(rtpos.Pos_lat)), median(deg2rad(rtpos.Pos_lon)), median(rtpos.Pos_hgt)];
    enu = LLH2ENU([deg2rad(rtpos.Pos_lat), deg2rad(rtpos.Pos_lon), rtpos.Pos_hgt], ref);
    rtpos.Pos_e = enu(:,1); rtpos.Pos_n = enu(:,2); rtpos.Pos_u = enu(:,3);
catch
    % Fallback if LLH2ENU missing: Simple deg to meters approx
    lat0 = median(rtpos.Pos_lat);
    lon0 = median(rtpos.Pos_lon);
    rtpos.Pos_n = (rtpos.Pos_lat - lat0) * 111132;
    rtpos.Pos_e = (rtpos.Pos_lon - lon0) * 111132 * cosd(lat0);
    rtpos.Pos_u = rtpos.Pos_hgt - median(rtpos.Pos_hgt);
end

% Time array
rtpos.caltime = gps_tow_to_utc_octave(rtpos.GPS_Sec, rtpos.GPS_Cycle);

dHz = sqrt(rtpos.Pos_n.^2 + rtpos.Pos_e.^2);
dHzSig = sqrt(rtpos.Sig_n.^2 + rtpos.Sig_e.^2);

% --- CSV EXPORT ---
cOutputDir = char(outputDir);
cBaseName  = char(strcat(name, '_RT_Pos'));
csvFileName  = fullfile(cOutputDir, [cBaseName, '.csv']);
if options.OutputCSV
     writePositionCSV(rtpos, csvFileName, options, metaData);
end

% --- PLOTTING (Octave Compatible) ---
if options.showFig == false
    set(0, 'DefaultFigureVisible', 'off');
end

% 1. Sat Map
plotSatMap(rtpos, name, utcTitleStr, outputDir, options);

% 2. Scatter
plotScatter_Octave(rtpos, dHz, name, utcTitleStr, outputDir, options);

% 3. Stats
plotHorizontalError_Octave(rtpos, dHz, dHzSig, name, utcTitleStr, outputDir, options);
plotVerticalError_Octave(rtpos, name, utcTitleStr, outputDir, options);

% 4. CDF
plotCDF_Octave(dHz, 'Horizontal Error [m]', 'Horizontal Error CDF', [name '_RT_CDF_Hz'], outputDir, name, utcTitleStr, options);
plotCDF_Octave(abs(rtpos.Pos_u), 'Vertical Error [m]', 'Vertical Error CDF', [name '_RT_CDF_Vert'], outputDir, name, utcTitleStr, options);

if options.showFig == false
    set(0, 'DefaultFigureVisible', 'on');
end

logMsg(sprintf('  Completed %s in %.2fs', name, toc(tFile)), options);
end


%=======================================================================
% PLOTTING FUNCTIONS
%=======================================================================

function plotSatMap(rtpos, name, utcTitleStr, outputDir, options)
    fig = figure('Visible', 'off');
    ax = axes(fig);
    plot(ax, rtpos.Pos_lon, rtpos.Pos_lat, 'b.');
    grid(ax, 'on');
    xlabel(ax, 'Longitude [deg]');
    ylabel(ax, 'Latitude [deg]');
    title(ax, ['Satellite View (Simple) - ' name], 'Interpreter', 'none');
    
    saveAndExport(fig, fullfile(outputDir, [name, '_RT_SatMap']), options);
end

function plotScatter_Octave(rtpos, dHz, name, utcTitleStr, outputDir, options)
    fig = figure('Visible', 'off');
    ax = axes(fig);
    plot(ax, rtpos.Pos_e, rtpos.Pos_n, 'b.');
    hold(ax, 'on');
    plot(ax, 0, 0, 'r+', 'MarkerSize', 10);
    grid(ax, 'on');
    axis(ax, 'equal');
    xlabel(ax, 'East [m]');
    ylabel(ax, 'North [m]');
    title(ax, ['Scatter Plot - ' name], 'Interpreter', 'none');
    
    saveAndExport(fig, fullfile(outputDir, [name, '_RT_scatter']), options);
end

function plotHorizontalError_Octave(rtpos, dHz, dHzSig, name, utcTitleStr, outputDir, options)
    fig = figure('Visible', 'off');
    
    ax1 = subplot(2,1,1);
    [t, y] = FillTimseriesGapsWithNaNs(rtpos.caltime, dHz);
    plot(ax1, t, y, 'b');
    ylabel(ax1, 'Hz Error [m]'); grid(ax1, 'on');
    title(ax1, ['Horizontal Error - ' name], 'Interpreter', 'none');
    
    ax2 = subplot(2,1,2);
    [t, y] = FillTimseriesGapsWithNaNs(rtpos.caltime, dHzSig);
    plot(ax2, t, y, 'm');
    ylabel(ax2, 'Hz Sigma [m]'); grid(ax2, 'on');
    xlabel(ax2, 'Time');
    
    saveAndExport(fig, fullfile(outputDir, [name, '_RT_Hz']), options);
end

function plotVerticalError_Octave(rtpos, name, utcTitleStr, outputDir, options)
    fig = figure('Visible', 'off');
    
    ax1 = subplot(2,1,1);
    [t, y] = FillTimseriesGapsWithNaNs(rtpos.caltime, rtpos.Pos_u);
    plot(ax1, t, y, 'b');
    ylabel(ax1, 'Vert Error [m]'); grid(ax1, 'on');
    title(ax1, ['Vertical Error - ' name], 'Interpreter', 'none');
    
    ax2 = subplot(2,1,2);
    [t, y] = FillTimseriesGapsWithNaNs(rtpos.caltime, rtpos.Sig_u);
    plot(ax2, t, y, 'm');
    ylabel(ax2, 'Vert Sigma [m]'); grid(ax2, 'on');
    xlabel(ax2, 'Time');
    
    saveAndExport(fig, fullfile(outputDir, [name, '_RT_Vert']), options);
end

function plotCDF_Octave(data, xLabelStr, titlePrefix, nameSuffix, fpath, dispName, utcTitleStr, options)
    data = data(~isnan(data));
    if isempty(data), return; end
    
    fig = figure('Visible', 'off');
    ax = axes(fig);
    hold(ax, 'on');
    
    % Use manual KDE smoothing for vertical data, standard for horizontal
    if contains(xLabelStr, 'Vertical', 'IgnoreCase', true)
         [f_smooth, x_smooth] = smoothCDF_KDE(data, 0.001, 1000);
         plot(ax, x_smooth, f_smooth * 100, 'b', 'LineWidth', 2);
    else
         sorted_data = sort(data);
         n = length(sorted_data);
         p = 100 * (1:n) / n;
         plot(ax, sorted_data, p, 'b', 'LineWidth', 2);
    end
    
    grid(ax, 'on');
    ylabel(ax, 'Probability [%]');
    xlabel(ax, xLabelStr);
    title(ax, [titlePrefix ' - ' dispName], 'Interpreter', 'none');
    
    saveAndExport(fig, fullfile(fpath, nameSuffix), options);
end

%=======================================================================
% HELPER FUNCTIONS
%=======================================================================

function saveAndExport(fig, baseFileName, options)
    % Octave save logic
    try
        fName = [baseFileName '.png'];
        print(fig, fName, '-dpng', sprintf('-r%d', options.ThumbRes));
    catch
        warning('Failed to save plot: %s', baseFileName);
    end
    close(fig);
end

function [t_out, y_out] = FillTimseriesGapsWithNaNs(timeinp, valinp)
    % Replaces timetable/retime with vector logic
    if isempty(timeinp), t_out=[]; y_out=[]; return; end
    
    dt = diff(timeinp);
    UR = median(dt(~isnan(dt)));
    if isnan(UR) || UR == 0, t_out = timeinp; y_out = valinp; return; end
    
    gap_threshold = 1.5 * UR; % define gap
    gaps = find(dt > gap_threshold);
    
    if isempty(gaps)
        t_out = timeinp; y_out = valinp;
        return;
    end
    
    t_out = [];
    y_out = [];
    last_idx = 1;
    
    for i = 1:length(gaps)
        idx = gaps(i);
        % Append chunk
        t_out = [t_out; timeinp(last_idx:idx); timeinp(idx)+UR]; % Insert time for gap
        y_out = [y_out; valinp(last_idx:idx); NaN];             % Insert NaN for gap
        last_idx = idx + 1;
    end
    
    t_out = [t_out; timeinp(last_idx:end)];
    y_out = [y_out; valinp(last_idx:end)];
end

function [f, x_grid] = smoothCDF_KDE(data, bandwidth, numPoints)
    % Manual Gaussian KDE
    data = data(:);
    n = length(data);
    x_min = min(data) - 3*bandwidth;
    x_max = max(data) + 3*bandwidth;
    x_grid = linspace(x_min, x_max, numPoints)';
    
    f = zeros(numPoints, 1);
    % Loop is slow in Octave, but reliable without toolboxes
    % Vectorized approach:
    % CDF of Gaussian: 0.5 * (1 + erf((x - mu) / (sigma * sqrt(2))))
    
    % To vectorize, we create a matrix (NumPoints x N), might be large
    if n * numPoints < 1e7
         X = repmat(x_grid, 1, n);
         MU = repmat(data', numPoints, 1);
         F_mat = 0.5 * (1 + erf((X - MU) / (bandwidth * sqrt(2))));
         f = mean(F_mat, 2);
    else
         for i = 1:n
            f = f + 0.5 * (1 + erf((x_grid - data(i)) / (bandwidth * sqrt(2))));
         end
         f = f / n;
    end
end

function allGPS_Sec = getFileGPSData(T0xFileName, options)
    % Extracts GPS seconds from file via viewdat
    [filepath, name, ext] = fileparts(T0xFileName);
    isT04 = contains(upper(ext), 'T04');
    if isT04, extOut = ".d35_2"; flag = "-d35:2"; else, extOut = ".d29"; flag = "-d29"; end

    outP = fullfile(filepath, [name, extOut]);
    if ~exist(outP, 'file')
        pathToViewDat = getPathToViewDat(); 
        cmd = sprintf('"%s" %s -h -m "%s" > "%s"', pathToViewDat, flag, T0xFileName, outP);
        [status, ~] = system(cmd);
        if status ~= 0, warning(['Decoding failed for ' T0xFileName]); allGPS_Sec = []; return; end
    end

    if exist(outP, 'file')
        if isT04
            try
                % TODO: Ensure import_d35_2 is available
                rtpos_temp = import_d35_2(outP, options.FixStatusManager);
            catch
                rtpos_temp = [];
            end
        else
            % TODO: Ensure ReadD29 is available
            rtpos_temp = ReadD29(outP);
        end
        if ~isempty(rtpos_temp) && isfield(rtpos_temp, 'GPS_Sec')
            allGPS_Sec = rtpos_temp.GPS_Sec; 
        else
            allGPS_Sec = []; 
        end
    else
        allGPS_Sec = []; 
    end
end

function s = subsetStruct(s, idx)
    fields = fieldnames(s);
    for k = 1:numel(fields)
        val = s.(fields{k});
        if isnumeric(val) && numel(val) == numel(idx)
             s.(fields{k}) = val(idx);
        elseif iscell(val) && numel(val) == numel(idx)
             s.(fields{k}) = val(idx);
        end
    end
end

function segments = splitDataByGap(allGPS_Sec, gapThreshold)
    if isempty(allGPS_Sec), segments = {}; return; end
    allGPS_Sec = allGPS_Sec(:); 
    dt = diff(allGPS_Sec); 
    gapIndices = find(dt > gapThreshold);
    startIndices = [1; gapIndices + 1]; 
    endIndices = [gapIndices; length(allGPS_Sec)];
    segments = cell(length(startIndices), 1);
    for i = 1:length(startIndices)
        segments{i} = struct('startTime', allGPS_Sec(startIndices(i)), 'endTime', allGPS_Sec(endIndices(i))); 
    end
end

function t_out = gps_tow_to_utc_octave(gps_tow, gps_week)
    % Simplified conversion: GPS epoch + weeks + seconds
    epoch = datenum(1980, 1, 6);
    leap = 18 / 86400; 
    t_out = epoch + (gps_week * 7) + (gps_tow / 86400) - leap;
end

function writePositionCSV(rtpos, fileName, options, meta)
    % Simplified CSV export for Octave
    fileName = char(fileName);
    fid = fopen(fileName, 'w');
    if fid == -1, return; end
    
    fprintf(fid, 'Filename:,%s\n', meta.Filename);
    fprintf(fid, 'Start Time:,%s\n', meta.Start.CalendarTime);
    fprintf(fid, 'GPS_Week,GPS_Sec,Lat_deg,Lon_deg,Hgt_m\n');
    
    n = length(rtpos.GPS_Sec);
    for i=1:n
        fprintf(fid, '%d,%.3f,%.9f,%.9f,%.3f\n', ...
            rtpos.GPS_Cycle(i), rtpos.GPS_Sec(i), rtpos.Pos_lat(i), rtpos.Pos_lon(i), rtpos.Pos_hgt(i));
    end
    fclose(fid);
end

function path = getPathToViewDat()
    if ispc
        path = 'ViewDAT.exe';
    else
        path = 'viewdat';
    end
    % Fallback checks can be added here
end

function logMsg(msg, options)
    fprintf('[%s] %s\n', datestr(now, 13), msg);
end

function printHelp()
    fprintf('ProcessRealtimePositions (Octave)\nUsage: function(file, "Key", Val...)\n');
end

% --- STUBS / TODOs ---
function generateReportsRouter(scanDir, options)
    % TODO: Implement HTML string generation for reports
    logMsg('Report generation (HTML) not fully implemented in Octave version.', options);
end

function md = initMetaData(filename, options)
    md = options.metaData;
    md.Filename = filename;
    md.Start.CalendarTime = '';
    md.Start.GPS_Sec = 0;
end

function info = getFileInfoOnly(filename, options)
    % TODO: Implement full header parsing
    info = struct();
end

function useIni = askUserWithTimeout(messageText, timeoutSecs)
    % Octave UI fallback: print to console
    fprintf('\n%s\n', messageText);
    fprintf('Defaulting to Antenna.ini in %d seconds...\n', timeoutSecs);
    useIni = true; 
    % pause(timeoutSecs); % Blocking pause, or just proceed
end

function [rtpos] = PreprocessPositionsData(rtpos, metaData, options)
    if isfield(rtpos, 'Pos_hgt') && ~isnan(metaData.AntennaApcHeight)
         % Logic placeholder
    end
end
