function [fitResults, rawCurves] = fitSingleGaussianCurvesToRadialSpectralDistribution(S_k_all, k_rho_vals, varargin)
%% fitSingleGaussianCurvesToRadialSpectralDistribution
% Fits a single-Gaussian model to multiple radial spectral curves (S(k_rho)),
% truncating each curve to user-specified k_rho range.
%
% Inputs:
%   S_k_all     - 1xN cell array of radial spectral curves
%   k_rho_vals  - vector of corresponding k_rho values
%
% Optional name-value pairs:
%   'KRhoRange'        - [k_min, k_max] (default: [min(k_rho_vals), max(k_rho_vals)])
%   'ResidualThreshold' - threshold for residual validation (default: 0.3)
%
% Outputs:
%   fitResults - struct array with fields:
%       .pFit, .kFit, .xFit, .yFit, .kFine, .isValid, .residual, .maxAbs, .R2, .fitMaxKRho
%   rawCurves  - struct array of truncated raw curves for plotting

    % --- Parse optional inputs ---
    p = inputParser;
    addParameter(p, 'KRhoRange', [min(k_rho_vals), max(k_rho_vals)], @(x) isnumeric(x) && numel(x)==2);
    addParameter(p, 'ResidualThreshold', 0.3, @(x) isnumeric(x) && isscalar(x));
    parse(p, varargin{:});
    opts = p.Results;

    Ncurves = numel(S_k_all);
    fitResults = struct('pFit',[],'kFit',[],'xFit',[],'yFit',[],...
                        'kFine',[],'isValid',[],'residual',[],'maxAbs',[],'R2',[],'fitMaxKRho',[]);
    rawCurves = struct('x', cell(1, Ncurves), 'k', cell(1, Ncurves));
    
    for k = 1:Ncurves
        % --- Truncate curve to k_rho and amplitude ranges ---
        xRaw    = S_k_all{k};
        % --- Normalize and smooth ---
        xNorm   = xRaw / max(xRaw);
        xSmooth = smooth(xNorm, 5, 'moving');
        
        mask  = (k_rho_vals(:) >= opts.KRhoRange(1)) & (k_rho_vals(:) <= opts.KRhoRange(2));
        
        x     = xSmooth(mask);
        kVals = k_rho_vals(mask);

        rawCurves(k).x = x;
        rawCurves(k).k = kVals;
        
        if any(isnan(xRaw))
            warning('Curve %d has NaNs, skipping.', k);
            fitResults(k) = fillNaNStructRadial();
            continue;
        elseif numel(x) < 5
            warning('Curve %d has too few points after truncation, skipping.', k);
            fitResults(k) = fillNaNStructRadial();
            continue;
        elseif isempty(x)
            warning('Curve %d is empty after truncation.', k);
            fitResults(k) = fillZeroStructRadial(kVals, opts.KRhoRange(2));
            continue;
        end

        try
            %% Single Gaussian fit with offset: A*exp(-(k-mu)^2/(2*sigma^2)) + C
            singleGauss = @(p,k) p(1) * exp(-0.5*((k - p(2))/max(p(3),1e-6)).^2) + p(4);
            
            % Initial guess
            [~, idxMax] = max(x);
            A_guess     = x(idxMax) - min(x);          % amplitude above baseline
            mu_guess    = kVals(idxMax);
            sigma_guess = 0.1 * (opts.KRhoRange(2)-opts.KRhoRange(1));
            C_guess     = min(x);                      % baseline offset
            p0          = [A_guess, mu_guess, sigma_guess, C_guess];
            
            lb = [0, opts.KRhoRange(1), 1e-6, -Inf];
            ub = [2*max(x), opts.KRhoRange(2), opts.KRhoRange(2), Inf];
            
            optsLSQ  = optimoptions('lsqcurvefit','Display','off', ...
                                       'MaxFunctionEvaluations',1e4,'MaxIterations',1e4);
            pFit     = lsqcurvefit(singleGauss, p0, kVals(:), x(:), lb, ub, optsLSQ);
            
            % --- Post-fit diagnostics ---
            r = x(:) - singleGauss(pFit, kVals(:));
            residualRMS = sqrt(mean(r.^2));
            maxAbsR = max(abs(r));
            
            kFine = linspace(min(kVals), max(kVals), 500);
            yFit  = singleGauss(pFit, kFine);
            fitMaxKRho = pFit(2);
            
            SS_res = sum(r.^2);
            SS_tot = sum((x(:) - mean(x(:))).^2);
            R2 = 1 - SS_res/SS_tot;
            
            % Store results
            fitResults(k).pFit       = pFit;
            fitResults(k).kFit       = kVals;
            fitResults(k).xFit       = x;
            fitResults(k).kFine      = kFine;
            fitResults(k).yFit       = yFit;
            fitResults(k).isValid    = residualRMS <= opts.ResidualThreshold;
            fitResults(k).residual   = residualRMS;
            fitResults(k).maxAbs     = maxAbsR;
            fitResults(k).R2         = R2;
            fitResults(k).fitMaxKRho = fitMaxKRho;
        catch
            warning('Curve %d fit failed completely.', k);
            fitResults(k) = fillNaNStructRadial();
        end
    end
end

%% --- Helper: NaN struct for failed fits ---
function s = fillNaNStructRadial()
    s = struct( ...
        'pFit',        nan(1,3), ...
        'kFit',        nan, ...
        'xFit',        nan, ...
        'yFit',        nan, ...
        'kFine',       nan, ...
        'isValid',     false, ...
        'residual',    nan, ...
        'maxAbs',      nan, ...
        'R2',          nan, ...
        'fitMaxKRho',  nan ...
     );
end

%% --- Helper: Zero struct for empty curves ---
function s = fillZeroStructRadial(kVals, MaxKRho)
    s = struct( ...
        'pFit',        zeros(1,3), ...
        'kFit',        kVals, ...
        'xFit',        zeros(size(kVals)), ...
        'yFit',        zeros(size(kVals)), ...
        'kFine',       zeros(1,500), ...
        'isValid',     false, ...
        'residual',    zeros(size(kVals)), ...
        'maxAbs',      zeros(size(kVals)), ...
        'R2',          zeros(size(kVals)), ...
        'fitMaxKRho',  MaxKRho ...
    );
end