//Lab 1: Odrazivost tenskej vrstvy
//Kolmy dopad

clear; clf;

//TO BE CHANGED

//Working directory must be edited as needed!
WorkDir = "C:\Users\juro\Dropbox\FEI\Prednasky\FKMM\Laboratorne ulohy\01 reflektometria"

//Directory of refractive index of the material used
MaterialDir = WorkDir;

//File names of the material refractive indices
SubIndexFile = "Si.csv";
LayIndexFile = "SiO2.csv";

//Directory of the measured data
DataDir = WorkDir+"\Measurement2";

//File names of the measured data
ReferenceFile = "reference2.txt";
SampleFile = "w1s2.txt";

h = 2; //Thickness of the reference SiO2 layer in nm

//Estimated values of the parameters for fitting
ThicknessEst = 70; //Estimation of the SiO2 layer thickness
I0Est = 0; //Estimation of the constant intensity

//Spectral range:
lambdamin = 350; // nm
lambdamax = 950; // nm

//Spectral resolution
lambdastep = 2; // nm

//END OF CHANGE

//Auxiliary functions

exec (WorkDir + "\Functions.sci");


//Local functions
function R = Rcalc(h, w, eOut, e1)
    //Function calculates reflectivity R of dielectric layer on substrate in normal incidence
    //h - thickness of the layer
    //w - wavelengths - array
    //eOut - complex permittivity of the superstrate (air) - array of the same length as w
    //e1 - complex permittivity of the layer - array of the same length as w
    for jj = 1:length(w)
      R(jj) = OneLayPerpReflectivity(1, 1, eOut(jj), 1, e1(jj), 1, h, w(jj));
    end
endfunction

function e = RefRes(abc, m)
    
    //function RefRes calculates the residui for the least-square fitting
    //abc -- array of the parameters to be fitted
    
    thickness = abc(1);
    I0 = abc(2); //I0 is the constant intensity shift caused by the experimental layout
    
    e = real(Rmeas) - real(I0 + Rcalc(thickness, lambda, epsSub, eps1));
    
endfunction
//

//Wavelength definition
lambda = (lambdamin:lambdastep:lambdamax)';
N = length(lambda);


//MATERIAL PARAMETERS
nSub = RefIndexInterp(WorkDir, SubIndexFile, lambda);
epsSub = nSub.^2;

n1 = RefIndexInterp(WorkDir, LayIndexFile, lambda);
eps1 = n1.^2;
h1 = h; //nm
//

//MEASURED DATA PREPROCESSING
Iref = DataInterp(DataDir, ReferenceFile, lambda);
Imeas = DataInterp(DataDir, SampleFile, lambda);
//

//Calculation of the reference sample reflectivity
Rref = Rcalc(h1, lambda, epsSub, eps1);
//

//Measured reflectivity of the SiO2/Si layer
Rmeas = (Imeas.*Rref)./Iref;
//

//least-square Lavenberg-Marquardt algorithm
[abc, v] = lsqrsolve([ThicknessEst; I0Est], RefRes, N);
//

//Results print to the terminal
disp("I0 = "+string(abc(2)), "h = "+string(abc(1))+" nm");
disp("chi2 = "+string(norm(v)));
//

//Plotting the results
plot2d(lambda,[Rmeas (abc(2)+Rcalc(abc(1),lambda,epsSub,eps1))],axesflag=1);

    //Titles
    xtitle('h = '+string(round(abc(1)))+' nm,   '+'chi2 = '+string(norm(abc(2))), 'Wavelength nm', 'Reflectivity');
    //

    //Legend
    legend("Measurement", "Non-linear fit",boxed=%f);
    //

    //Ploted lines properties
    GrphA = gca();
    GrphA.children(2).children(1).thickness = 1.5;
    GrphA.children(2).children(1).foreground = 5;
    GrphA.children(2).children(1).line_style = 1;
    GrphA.children(2).children(2).thickness = 4;
    GrphA.children(2).children(2).foreground = 2;
    GrphA.children(2).children(2).line_style = 7;
    GrphA.children(1).font_size = 2;
    GrphA.font_size = 2;
    //
//