num_params = 3;

f = fits.exp{1};
model_func = @(x)exponential_bscatter(x,f.diel,f.s,f.A);
chi.exp = chi_square(model_func, s, num_params);

f = fits.gauss{1};
model_func = @(x)gaussian_bscatter(x,f.diel,f.s,f.A);
chi.gauss = chi_square(model_func, s, num_params);

f = fits.hagfors{1};
model_func = @(x)hagfors(x,f.diel,f.s,f.A);
chi.hagfors = chi_square(model_func, s, num_params);

num_params = 5;
f = fits.simpleTwoLayer{1};
model_func = @(x)simple_two_layer(x,f.eps1,f.eps2,f.s1,f.s2,f.A);
chi.two_layer = chi_square(model_func, s, num_params);

f = fits.simpleTwoLayer{2};
model_func = @(x)simple_two_layer(x,f.eps1,f.eps2,f.s1,f.s2,f.A);
chi.two_layer_limited = chi_square(model_func, s, num_params);




function [chi_squared_reduced, chi_squared] = chi_square(model_func, s, num_params)
    % CHI_SQUARE Calculate the reduced chi-squared metric for model fitting
    %   Inputs:
    %   - model_func: Anonymous function of the model to be evaluated
    %   - s: Struct containing the data (as described in the problem)
    %   - num_params: Number of free parameters in the model
    %
    %   Outputs:
    %   - chi_squared_reduced: The reduced chi-squared value
    %   - chi_squared: The chi-squared value (before dividing by degrees of freedom)

    % Initialize variables to store combined data
    all_inc = [];
    all_sigma0 = [];
    all_sigma0_std = [];

    % Combine data from all sites
    for idx = 1:numel(s)
        ave_inc = mean(s(idx).inc);
        ave_sigma0 = mean(s(idx).sigma0);
        std_sigma0 = std(s(idx).sigma0);
        
        all_inc = [all_inc  ave_inc];
        all_sigma0 = [all_sigma0  ave_sigma0];
        all_sigma0_std = [all_sigma0_std  std_sigma0];
        
        if ~isempty(s(idx).alt_inc)

        ave_alt_inc = mean(s(idx).alt_inc);
        ave_alt_sigma0 = mean(s(idx).alt_sigma0);
        std_alt_sigma0 = std(s(idx).alt_sigma0);
        
        all_inc = [all_inc  ave_alt_inc];
        all_sigma0 = [all_sigma0  ave_alt_sigma0];
        all_sigma0_std = [all_sigma0_std  std_alt_sigma0];

        end

    end

    all_inc =  deg2rad(all_inc); %s is stored in degrees for some rea

    % Calculate model predictions
    model_predictions = model_func(all_inc);

    % Calculate chi-squared
    chi_squared = sum(((all_sigma0 - model_predictions) ./ all_sigma0_std).^2);

    % Calculate degrees of freedom
    dof = numel(all_inc) - num_params;

    % Calculate reduced chi-squared
    chi_squared_reduced = chi_squared ./ dof;
end