%%This file generates the function MC_simulation which is called in the main code MC_smm_main.m (lines 50 - 56)
%The function uses number of runs, number of basefunctions N, wavelength, and input data (E_meas,E_ref,S_ref, and S_DUT) as inputs from the main code file
%and returns the standard uncertainty /% and the spectral mismatch factor. 
%First the wavelength is called as 'wl' at line 14, then input parameter to be analysed is selected as 'spec' in lines 19 - 21. 'unc' calls the input standard uncertainty from lines 26 - 32
%in the main code. An orthogonal function is generated; lines 28 - 29. A loop starts on line 38 and runs for 50000 times. Normally distributed r is generated on lines 40 - 47. A secondary loop
%starts at 56.   
%Weights are then calculated as Equation 5 in the article. These weights are then used to calculate the error function (for i = 1,2,3,... N.) given as Equation 3 in the article. 
%The error function delta is then used in Equation 2 and line 70 to form
%a wavelength dependent error, irre. This wavelength depenedent error is then called in Equation 1 lines 75 - 84 to calculate the spectral mismatch error and to analyse the 
%uncertainties. The primary loop ends after this at line 85. The standard uncertainty is obtained by simply calculating the uncertainty between 50000 spectral mismatch values. 
function [std_uncertainty, mismatch] = MC_simulation(runs, n_basefunc, wavelength, ... 
     E_meas, E_ref,  S_ref,S_DUT,  standard_uncertainty)
    % Monte Carlo function for estimating the standard uncertainties 
    wl = wavelength;      % Wavelength   
    %Selecting which parameter to analyse. This is made one parameter at a time. Options are: (spec = S_ref or spec = S_DUT or spec = E_meas) 
    
    %Uncomment one of the following three lines 
    %depending on the parameter you want to evaluate:
    spec =  S_DUT; %uncomment to analyse uncertainty caused by S_DUT
    %spec = S_ref %uncomment to analyse uncertainty caused by S_ref
    %spec = E_meas %uncomment to analyse uncertainty caused by SNR, radiocalibration, stability, bandwidth or wavelength. These are the subcomponents of the uncertainty due to the 
		   %spectral irradiance of solar cell (E_meas). Therefore, they are only considered in the distorted spectral irradiance. The parameter (subcomponent) being analysed here 
		   %depends on the parameter selected in the main code file in lines 26 - 30. 		   

    unc = standard_uncertainty; % Standard uncertainty 
    
    % ----- Function for creating the base functions (Equation 4 in the article) ----- 
    create_base_function = @(wl, phi, order) ... 
        (sqrt(2)*sin(order*2*pi*((wl - wl(1))./(wl(end) - wl(1))) + phi));
    % ----- 

    % Initialize vectors 
    irre = zeros(runs, 1); 
    smm = zeros(runs,1);

    % Display the base functions used in the current simulation in the command window
    disp(['Order of the base function: ',  num2str(n_basefunc)]);

    % ----- Monte Carlo starts ----- 
    for idx_runs = 1:runs  %runs for 50000 times

        % Normal distributed r
        r = [random('Normal', 0, 1, n_basefunc, 1)]; 
        
        R = sqrt(sum(r.^2));

        a = r./R; 
        
        phase = [random('Uniform', 0, 2*pi, n_basefunc, 1)];
        
        delta_smm = zeros(length(wl), 1);
        function_smm = zeros(length(wl), 1);
     

        order = 1; 
        for idx = 1:n_basefunc
            phi = phase(n_basefunc); 
            if idx == 1
                delta_smm = delta_smm + a(n_basefunc)*ones(length(wl), 1);           
            else 
                function_smm = create_base_function(wl, phi, order);
                delta_smm = delta_smm + a(n_basefunc)*function_smm;
                order = order + 1; 
            end
  
        end
        delta = delta_smm;
		
           % Distorted spectrum
           irre = (1 + delta.*unc).*spec;  %Equation 2 in article
           
		   %Selecting which parameter to analyse.
           %Spectral mismatch (the parameter under test has to be replaced by distorted irre function in the equation). This corresponds to Equation 1 in article
           
		   %uncomment to analyse uncertainty due to spectral responsivity of cell (device) under test
           smm(idx_runs) = ((trapz(E_ref.*S_ref)).*(trapz(E_meas.*irre)))./((trapz(E_meas.*S_ref)).*trapz(E_ref.*irre)); %distorted S_DUT
           
           %uncomment to analyse uncertainty due to SNR, Radiometric calibration,stability, bandwidth or wavelength. These are the subcomponents of the uncertainty due to the 
		   %spectral irradiance of solar cell (E_meas). Therefore, they are only considered in the distorted spectral irradiance. The parameter being analysed here depends on the 
           %parameter selected in the main code file in lines 26 - 30. 		   
           %smm(idx_runs) = ((trapz(E_ref.*S_ref)).*(trapz(irre.*S_DUT)))./((trapz(irre.*S_ref)).*trapz(E_ref.*S_DUT));
           
           %uncomment to calcalute uncertainty due to spectral responsivity of reference cell
           %smm(idx_runs) = ((trapz(E_ref.*irre)).*(trapz(E_meas.*S_DUT)))./((trapz(E_meas.*irre)).*trapz(E_ref.*S_DUT));
    end 
    % ----- 
    
    % ----- Calculate standard uncertainty ----- 
    std_uncertainty = std(smm)*100;  %Standard uncertainty 
    mismatch = mean(smm); %Spectral mismatch value
    % ----- 
    
    
    % Display the standard uncertainty in the command window
    disp(['Standard uncertainty: ',  num2str(std_uncertainty)]);
    disp([' ']); 
    
%%End of function file