function [fitResults, rawCurves] = fitTwoGaussianCurvesToRadialSpectralDistribution(S_k_all, k_rho_vals, varargin)
%% fitTwoGaussianCurvesToRadialSpectralDistribution
% Fits a two-Gaussian model to multiple radial spectral curves (S(k_rho)),
% truncating each curve to user-specified k_rho and amplitude ranges.
%
% Author:       Karthik
% Date: 2025-10-15
% Version: 1.0
%
% 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)])
%   'AmplitudeRange'   - [y_min, y_max] (default: [0, Inf])
%   'ResidualThreshold', 'PositionThreshold', 'AmplitudeThreshold' as before
%
% Outputs:
%   fitResults - struct array with fields:
%       .pFit, .kFit, .xFit, .yFit, .kFine, .isValid, .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, 'AmplitudeRange', [0, Inf], @(x) isnumeric(x) && numel(x)==2);
    addParameter(p, 'ResidualThreshold', 0.3, @(x) isnumeric(x) && isscalar(x));
    addParameter(p, 'PositionThreshold', 0.1, @(x) isnumeric(x) && isscalar(x));
    addParameter(p, 'AmplitudeThreshold', 0.5, @(x) isnumeric(x) && isscalar(x));
    parse(p, varargin{:});
    opts = p.Results;

    Ncurves = numel(S_k_all);
    fitResults = struct('pFit',[],'kFit',[],'xFit',[],'yFit',[],...
                        'kFine',[],'isValid',[],'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)) & ...
               (xSmooth(:) >= opts.AmplitudeRange(1)) & (xSmooth(:) <= opts.AmplitudeRange(2));

        x     = xSmooth(mask);
        kVals = k_rho_vals(mask);

        rawCurves(k).x = x;
        rawCurves(k).k = kVals;

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

        try
            % --- Detect peaks ---
            curveAmp           = max(x) - min(x);
            minProm            = opts.AmplitudeThreshold * curveAmp;
            [pk, locIdx]       = findpeaks(x, 'MinPeakProminence', minProm);
            kPeaks             = kVals(locIdx);

            %% Generic two-Gaussian
            twoGauss = @(p,k) ...
                p(1)*exp(-0.5*((k - p(2))/max(p(3),1e-6)).^2) + ...
                p(4)*exp(-0.5*((k - p(5))/max(p(6),1e-6)).^2);
            
            secondaryIdx = find(kPeaks > opts.PositionThreshold, 1, 'first');
            if ~isempty(secondaryIdx)
                kSecondary     = kPeaks(secondaryIdx);
                secondaryAmp   = pk(secondaryIdx);
                A1_guess     = opts.AmplitudeRange(2);
                mu1_guess    = 0.0;
                sigma1_guess = 0.1 * kSecondary;
                A2_guess     = secondaryAmp;
                mu2_guess    = kSecondary;
                sigma2_guess = 0.1 * kSecondary;
                p0           = [A1_guess, mu1_guess, sigma1_guess, A2_guess, mu2_guess, sigma2_guess];
            else
                p0           = [opts.AmplitudeRange(2), 0, 0.1, 0.01, opts.KRhoRange(2)/3, 0.1];
            end

            lb       = [0, 0, 1e-6, 0, opts.PositionThreshold, 1e-6];
            ub       = [2, opts.KRhoRange(2)/2, opts.KRhoRange(2), 2, opts.KRhoRange(2), opts.KRhoRange(2)];
            optsLSQ  = optimoptions('lsqcurvefit','Display','off', ...
                                       'MaxFunctionEvaluations',1e4,'MaxIterations',1e4);
            pFit     = lsqcurvefit(twoGauss, p0, kVals(:), x(:), lb, ub, optsLSQ);

            if pFit(4) <= 1e-3
                pFit(4)   = 0;
                pFit(5:6) = NaN;
            end
            
            fitMaxKRho               = max(pFit(5), pFit(2));
            kFine                    = linspace(min(kVals), max(kVals), 500);
            yFit                     = twoGauss(pFit, kFine);

            fitResults(k).pFit       = pFit;
            fitResults(k).kFit       = kVals;
            fitResults(k).xFit       = x;
            fitResults(k).kFine      = kFine;
            fitResults(k).yFit       = yFit;
            fitResults(k).isValid    = true;
            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,6), ...
        'kFit',        nan, ...
        'xFit',        nan, ...
        'yFit',        nan, ...
        'kFine',       nan, ...
        'isValid',     false, ...
        'fitMaxKRho',  nan ...
     );
end

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