Calculations/Data-Analyzer/+Analyzer/fitTwoGaussianCurves.m

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.0
%
% 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)
%
% Output:
%   fitResults - struct array with fields:
%       .pFit        - fitted parameters
%       .thetaFit    - theta values used for fit
%       .xFit        - curve values used for fit
%       .yFit        - fitted curve evaluated on thetaFine
%       .thetaFine   - fine theta for plotting
%       .isValid     - whether fit passed deviation threshold
%       .fitMaxTheta - theta used as fit limit

    % --- 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));
    parse(p, varargin{:});
    opts = p.Results;

    Ncurves = numel(S_theta_all);
    fitResults = struct('pFit',[],'thetaFit',[],'xFit',[],'yFit',[],...
                        'thetaFine',[],'isValid',[],'fitMaxTheta',[]);

    % --- 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};

        % --- Find first peak in 0–MaxTheta range ---
        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

        % Take only the part from the peak to the end
        S_theta_all_shifted{k} = curve(peak_idx:end);
    end

    % --- Pad curves to the 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

    % --- Define recentered theta values ---
    theta_recentered = theta_vals(1:Npoints_shifted) - theta_vals(1);

    for k = 1:Ncurves
        x        = S_theta_all_shifted{k};  % use recentered curve
        theta    = theta_recentered;     % recentered x-values

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

        % --- Restrict to 0–MaxTheta ---
        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);
            continue;
        end

        % --- Normalize so S(0)=1 ---
        x            = x / x(1);

        % --- Smooth for stable peak detection ---
        xSmooth      = smooth(x, 5, 'moving');

        % --- Find local maxima ---
        [pk, locIdx] = findpeaks(xSmooth);
        thetaPeaks   = theta(locIdx);

        if isempty(pk)
            warning('Curve %d has no significant peaks, restricting to within pi/4.', k);
            fitMaxTheta = pi/4;
        else
            % --- Primary peak ---
            [mainAmp, mainIdx] = max(pk);
            mainTheta          = thetaPeaks(mainIdx);
    
            thetaPeaks         = thetaPeaks(:);  % column vector
            pk                 = pk(:);
    
            % --- Secondary peak ≥50% of primary after main ---
            secondaryIdx = find((thetaPeaks>mainTheta) & (pk>=0.50*mainAmp), 1, 'first');
            if isempty(secondaryIdx)
                fitMaxTheta = opts.MaxTheta;
            else
                fitMaxTheta = thetaPeaks(secondaryIdx);
            end
        end

        % --- Extract data up to secondary peak ---
        fitMask = theta>=0 & theta<=fitMaxTheta;
        thetaFit = theta(fitMask);
        xFit = x(fitMask);
        if numel(thetaFit) < 8
            warning('Curve %d has too few points for fitting, skipping.', k);
            continue;
        end

        % --- Two-Gaussian model ---
        twoGauss = @(p,theta) ...
            1.0*exp(-0.5*((theta-p(1))/max(p(2),1e-6)).^2) + ...
            p(3)*exp(-0.5*((theta-p(4))/max(p(5),1e-6)).^2);
        % p = [mu1, sigma1, A2, mu2, sigma2]

        % --- Initial guesses ---
        mu1_guess = 0;
        sigma1_guess = 0.1*fitMaxTheta;
        A2_guess = max(xFit)*0.8;
        mu2_guess = fitMaxTheta/2;
        sigma2_guess = fitMaxTheta/5;
        p0 = [mu1_guess, sigma1_guess, A2_guess, mu2_guess, sigma2_guess];

        % --- Bounds ---
        lb = [0, 1e-6, 0, 0, 1e-6];
        ub = [fitMaxTheta/4, fitMaxTheta, 2, fitMaxTheta, 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, using initial guess.', k);
            pFit = p0;
        end

        % --- Evaluate fitted curve ---
        thetaFine = linspace(0, opts.MaxTheta, 500);
        yFit = twoGauss(pFit, thetaFine);

        % --- Compute relative deviation ---
        yFitInterp = twoGauss(pFit, thetaFit);
        relDeviation = abs(yFitInterp - xFit) ./ max(xFit, 1e-6);
        maxRelDev = max(relDeviation);

        % --- Goodness-of-fit ---
        isValid = maxRelDev <= 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;
    end
end