% ----- Main Script - Monte Carlo Method to Estimate Uncertainties in SMM ----- 
% ----- 
clear all; close all; clc;
set(0, 'DefaultTextFontname', 'Times New Roman', ... 
    'DefaultAxesFontsize', 14, 'DefaultTextFontsize', 14);
set(0, 'DefaultAxesFontName', 'Times New Roman', ... 
    'DefaultAxesFontsize', 14, 'DefaultTextFontsize', 14);
% ----- 
% Import files
E_sun_data = importdata('E_sun.txt', '\t', 1); %spectral irradiance AM1.5G 
E_meas_data = importdata('E_meas.txt', '\t', 1); %measured spectral irradiance
S_ref_data = importdata('S_ref.txt', '\t', 1); %spectral responsivity of reference cell
S_DUT_data = importdata('S_DUT.txt','\t', 1); %spectral responsivity of device (cell) under test

%Interpolation of data
lambda = [290:1:1200]'; % Wavelength (nm)
wl_E_sun = E_sun_data.data(:, 1); 
E_ref = spline(wl_E_sun, E_sun_data.data(:, 2),lambda); %AM1.5G
wl_E_meas = E_meas_data.data(:, 1);
E_meas = spline(wl_E_meas, E_meas_data.data(:, 2),lambda); %irradiance
wl_S_ref = S_ref_data.data(:, 1);
S_ref = spline(wl_S_ref,S_ref_data.data(:, 2),lambda); %spectral responsivity of reference cell
wl_S_DUT = S_DUT_data.data(:, 1);
S_DUT = S_DUT_data.data(:,2); %extrapolated spectral responsivity of detector under test
%Selecting which parameter to analyse.
%uncertainty = spline(wl_E_meas, (E_meas_data.data(:, 3)./200),lambda) %uncomment to to analyse uncertainties caused by SNR
%uncertainty = spline(wl_E_meas, (E_meas_data.data(:, 4)./200),lambda) %uncomment to analyse uncertainties caused by Radiometric Calibration
%uncertainty = spline(wl_E_meas, (E_meas_data.data(:, 5)./200),lambda)%uncomment to analyse uncertainties caused by Stability
%uncertainty = spline(wl_E_meas, (E_meas_data.data(:, 6)./200),lambda) %uncomment to analyse uncertainties caused by Bandwidth
%uncertainty = spline(wl_E_meas, (E_meas_data.data(:, 7)./200),lambda) %uncomment to analyse uncertainties caused by Wavelength
%uncertainty = spline(wl_S_ref, (S_ref_data.data(:, 4)./200),lambda); %uncomment to analyse uncertainties caused by spectral responsivity of reference cell
uncertainty = spline(wl_S_DUT, (S_DUT_data.data(:, 4)./200),lambda); %uncomment to analyse uncertainties caused by spectral responsivity of the cell under test

% Number of Monte Carlo rounds 
runs = 50000;

%Making a list of orders of complexity N to analyse. Highest order possible is 456 due to Nyquist criterion
basefunc = [0:10 16 32 50 64 75 100 128 150 175 200 256 280 300 375 456]; 

% Standard uncertainty of input parameter
standard_uncertainty = uncertainty;

% Select the step size, e.g. step = 1 is 1 nm, step = 5 is 5 nm
step = 1; 
wavelength = lambda(1:step:end);

%The MC_simulation.smm function is called here. This loop runs the function
%file and returns the spectral mismatch error and standard uncertainty / % for
%basefunctions 0 to 456
idx = 1;
for n_basefunc = basefunc
    [unc(idx),mismatch(idx)] ... 
        = MC_simulation(runs, n_basefunc, wavelength,   E_meas,  E_ref, S_ref, S_DUT,  ... 
        standard_uncertainty);
    idx = idx + 1;
end 

%The file is saved as a .txt file where the first column is the number of basefunctions N and the second column is the standard uncertainty in / %
MT = [(basefunc)' unc']; 
save uncertainty_table.txt MT -ascii

% ----- 
%A plot is generated with number of base functions N on the x-axis and the standard uncertainty / % on the y-axis for the input parameter under study.
figure;
set(gcf, 'color', 'white');
plot(basefunc, unc);
title('Spectral responsivity of DUT') %title corresponds to the parameter being analysed 
ylabel('Standard uncertainty / %');
xlabel('Number of base functions {\it N}');

%%End of code
