function loss_curve_digitized(filename)
% PLOT_WPD_DATASETS Reads and plots scattering data from a CSV file.
%   filename: String containing the path to the CSV file (e.g., 'wpd_datasets.csv')
%
%   This function expects a CSV with the following structure:
%   Row 1: Dataset Names (every 2nd column)
%   Row 2: Variable Names (X, Y, X, Y...)
%   Row 3+: Numeric Data (potentially with NaNs for unequal lengths)

    % --- 1. Import Data ---
    if nargin < 1
        filename = 'wpd_datasets.csv'; % Default filename
    end

    % Read the numeric data, skipping the first 2 header lines
    % 'TreatAsMissing' ensures empty cells become NaNs
    raw_data = readmatrix(filename, 'NumHeaderLines', 2);

    % Read the first line separately to parse Dataset Names
    fid = fopen(filename, 'r');
    if fid == -1
        error('Could not open file: %s', filename);
    end
    header_line = fgetl(fid);
    fclose(fid);
    
    % Split header by comma to get names
    raw_names = split(header_line, ',');
    
    % Extract non-empty names (assuming names are in col 1, 3, 5...)
    dataset_names = raw_names(~cellfun('isempty', raw_names));
    
    % --- 2. Setup Plot ---
    figure('Color', 'w', 'Position', [100, 100, 800, 600]);
    ax = gca;
    hold(ax, 'on');
   
    
    line_styles = {'-', '-', '-', '-'};
    
    % --- 3. Iterate and Plot Each Dataset ---
    num_datasets = length(dataset_names);
    
    for i = 1:num_datasets
        % Calculate column indices for X and Y
        % Dataset 1: Cols 1,2 | Dataset 2: Cols 3,4 | etc.
        col_x = (i-1)*2 + 1;
        col_y = (i-1)*2 + 2;
        
        % Extract data
        if col_y > size(raw_data, 2)
            warning('Data columns missing for dataset %d', i);
            break;
        end
        
        X = raw_data(:, col_x);
        Y = raw_data(:, col_y);

        
        
        % Remove NaNs (missing data due to unequal lengths)
        valid_mask = ~isnan(X) & ~isnan(Y);
        X = X(valid_mask);
        Y = Y(valid_mask);

        [X, sortIdx] = sort(X);
        Y = Y(sortIdx);
        
        % 3. Smooth the Data

        if numel(X) > 20
        if 1
            Y = smoothdata(Y, 'sgolay', 15); 
        else
            Y = movmean(Y, 15); 
        end
        end

        
        % Plotting
        % Using semilogy because scattering data often spans orders of magnitude
        % (Adjust to 'plot' if linear scale is preferred)
        p = plot(ax, X, Y, ...
            'LineStyle', line_styles{mod(i-1, length(line_styles)) + 1}, ...
            'Marker', 'none', ...
            'Color', 'black', ...
            'LineWidth', 1.5, ...
            'MarkerSize', 6, ...
            'DisplayName', dataset_names{i});
            
        % Optional: Fill marker faces for better visibility
        % p.MarkerFaceColor = p.Color;
        % p.MarkerFaceAlpha = 0.3; % Semi-transparent fill
    end
    
    % --- 4. Styling and Formatting ---
    
    % Axis Labels (Inferred from typical scattering plots)
    xlabel(ax, 'Wavelength (\mu m)', 'FontSize', 12, 'FontWeight', 'bold');
    ylabel(ax, 'Intensity / Cross-Section (a.u.)', 'FontSize', 12, 'FontWeight', 'bold');
    
    % Title
    title(ax, 'Dataset Comparison', 'FontSize', 14);
    
    % Legend
    legend(ax, 'Location', 'best', 'Interpreter', 'none', 'Box', 'on');
    
    % Grid
    grid(ax, 'on');
    ax.GridAlpha = 0.3;
    ax.MinorGridAlpha = 0.1;
    
    % Set Log Scale for Y (likely required for this data type)
    set(ax, 'YScale', 'log'); 
    
    % Enhance axis appearance
    set(ax, 'Box', 'on', 'LineWidth', 1.2, 'FontSize', 10);
    
    hold(ax, 'off');
end