function fitResults = fitTwoGaussianCurves(S_theta_all, theta_vals, varargin)
%% fitTwoGaussianCurves
% Fits a two-Gaussian model to multiple spectral curves from θ=0 up to
% the first secondary peak ≥50% of the primary peak, within 0–π/2 radians.
%
% Author: Karthik
% Date: 2025-10-06
% Version: 1.1
%
% Inputs:
%   S_theta_all - 1xN cell array of spectral curves
%   theta_vals  - vector of corresponding theta values (radians)
%
% Optional name-value pairs:
%   'MaxTheta'           - Maximum theta for search (default: pi/2)
%   'DeviationLimit'     - Max relative deviation allowed (default: 0.3)
%   'PositionThreshold'  - Minimum angular separation between main and
%                          secondary peaks (default: 0.1 rad)
%   'AmplitudeThreshold' - Minimum fraction of primary peak amplitude for secondary peak (default: 0.5)
%
% Output:
%   fitResults - struct array with fields:
%       .pFit, .thetaFit, .xFit, .yFit, .thetaFine, .isValid, .fitMaxTheta

    % --- Parse optional inputs ---
    p = inputParser;
    addParameter(p, 'MaxTheta', pi/2, @(x) isnumeric(x) && isscalar(x));
    addParameter(p, 'DeviationLimit', 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));
    addParameter(p, 'RecenterCurves', true, @islogical);
    parse(p, varargin{:});
    opts = p.Results;

    Ncurves = numel(S_theta_all);
    fitResults = struct('pFit',[],'thetaFit',[],'xFit',[],'yFit',[],...
                        'thetaFine',[],'isValid',[],'fitMaxTheta',[]);
    
    if opts.RecenterCurves
        % --- Preprocess curves: shift first peak to zero without wrapping ---
        S_theta_all_shifted = cell(size(S_theta_all));
        for k = 1:Ncurves
            curve = S_theta_all{k};
    
            idx_range = find(theta_vals >= 0 & theta_vals <= opts.MaxTheta);
            if isempty(idx_range)
                [~, peak_idx] = max(curve);
            else
                [~, local_max_idx] = max(curve(idx_range));
                peak_idx = idx_range(local_max_idx);
            end
    
            S_theta_all_shifted{k} = curve(peak_idx:end);
        end
    
        % --- Pad curves to same length ---
        Npoints_shifted = max(cellfun(@numel, S_theta_all_shifted));
        for k = 1:Ncurves
            len = numel(S_theta_all_shifted{k});
            if len < Npoints_shifted
                S_theta_all_shifted{k} = [S_theta_all_shifted{k}, nan(1, Npoints_shifted-len)];
            end
        end
    
        theta_recentered = theta_vals(1:Npoints_shifted) - theta_vals(1);
    else
        S_theta_all_shifted = S_theta_all;
        theta_recentered    = theta_vals;
    end

    for k = 1:Ncurves
        try
            x        = S_theta_all_shifted{k};
            theta    = theta_recentered;

            validIdx = ~isnan(x);
            x        = x(validIdx);
            theta    = theta(validIdx);

            mask  = theta>=0 & theta<=opts.MaxTheta;
            x     = x(mask);
            theta = theta(mask);
            if numel(theta) < 8
                warning('Curve %d too short (<8 points), skipping.', k);
                fitResults(k) = fillNaNStruct();
                continue;
            end

            x        = x / x(1);
            xSmooth  = smooth(x, 5, 'moving');
            curveAmp = max(xSmooth) - min(xSmooth);
            
            % --- Robust peak detection ---
            [pk, locIdx] = findpeaks(xSmooth, 'MinPeakProminence', opts.AmplitudeThreshold * curveAmp);
            thetaPeaks   = theta(locIdx);
            
            if isempty(pk)
                warning('Curve %d has no significant peaks.', k);
                fitResults(k) = fillZeroStruct(theta, opts.MaxTheta);
                continue;
            end
            
            thetaPeaks = thetaPeaks(:);
            pk         = pk(:);
            
            % --- Find secondary peak ---
            secondaryIdx = find((thetaPeaks > opts.PositionThreshold) & ...
                                (pk >= opts.AmplitudeThreshold * curveAmp), 1, 'first');
            
            if isempty(secondaryIdx)
                % No secondary peak → bypass fit, fill with zeros
                fitResults(k) = fillZeroStruct(theta, opts.MaxTheta);
                continue;
            else
                fitMaxTheta = thetaPeaks(secondaryIdx);
                secondaryAmp = pk(secondaryIdx);
            end
            
            % --- Use the full curve for fitting ---
            thetaFit = theta(:);
            xFit = x(:);
            
            if numel(thetaFit) < 8
                warning('Curve %d has too few points for fitting.', k);
                fitResults(k) = fillNaNStruct();
                continue;
            end
            
            % --- Two-Gaussian model ---
            twoGauss = @(p,theta) ...
                1.0 * exp(-0.5*((theta - p(1))/max(p(2),1e-6)).^2) + ...     % main Gaussian
                p(3) * exp(-0.5*((theta - p(4))/max(p(5),1e-6)).^2);         % secondary Gaussian
            
            % --- Parameter guesses based on secondary peak ---
            mu1_guess    = 0;
            sigma1_guess = 0.1 * fitMaxTheta;
            A2_guess     = secondaryAmp;             % use detected secondary amplitude
            mu2_guess    = fitMaxTheta;              % center at detected secondary peak
            sigma2_guess = fitMaxTheta / 3;          % broad guess
            
            p0           = [mu1_guess, sigma1_guess, A2_guess, mu2_guess, sigma2_guess];
            
            % --- Bounds that enforce secondary peak constraints ---
            lb = [0, 1e-6, 0, max(fitMaxTheta - 0.1, 0), 1e-6];
            ub = [fitMaxTheta/4, fitMaxTheta, 2, min(fitMaxTheta + 0.1, opts.MaxTheta), fitMaxTheta];
            
            optsLSQ = optimoptions('lsqcurvefit','Display','off', ...
                                   'MaxFunctionEvaluations',1e4,'MaxIterations',1e4);
            
            % --- Fit ---
            try
                pFit = lsqcurvefit(twoGauss, p0, thetaFit, xFit, lb, ub, optsLSQ);
            catch
                warning('Curve %d fit failed.', k);
                fitResults(k) = fillNaNStruct();
                continue;
            end
            
            % --- Evaluate fitted curve ---
            thetaFine = linspace(0, opts.MaxTheta, 500);
            yFit = twoGauss(pFit, thetaFine);
            
            % --- Assess deviation ---
            relDevAmplitude   = abs(pFit(3) - secondaryAmp) / max(secondaryAmp,1e-6);
            relDevPosition    = abs(pFit(4) - fitMaxTheta) / max(fitMaxTheta, 1e-6);
            combinedDeviation = sqrt(relDevAmplitude^2 + relDevPosition^2);
            isValid           = combinedDeviation <= opts.DeviationLimit;
            
            % --- Store results ---
            fitResults(k).pFit = pFit;
            fitResults(k).thetaFit = thetaFit;
            fitResults(k).xFit = xFit;
            fitResults(k).thetaFine = thetaFine;
            fitResults(k).yFit = yFit;
            fitResults(k).isValid = isValid;
            fitResults(k).fitMaxTheta = fitMaxTheta;
        
        catch
            warning('Curve %d fit failed completely.', k);
            fitResults(k) = fillNaNStruct();
        end
    end
end

% --- Helper: return struct with NaNs (fit failed) ---
function s = fillNaNStruct()
    s = struct('pFit',nan(1,5),'thetaFit',nan,'xFit',nan,'yFit',nan,'thetaFine',nan,...
               'isValid',false,'fitMaxTheta',nan);
end

% --- Helper: return struct with zeros (no secondary peak found) ---
function s = fillZeroStruct(theta, fitMaxTheta)
    s = struct('pFit',zeros(1,5),'thetaFit',theta,'xFit',zeros(size(theta)),...
               'yFit',zeros(size(theta)),'thetaFine',zeros(1,500),...
               'isValid',false,'fitMaxTheta',fitMaxTheta);
end