clear;
clc
warning off MATLAB:divideByZero;
% This program calculates contours in the n,k-plane pertaining to
% experimental input values of transmittance and reflectance. Hence the
% program calculates all n,k values that give those values of R and T. The
% intersections between the R and T contours give the possible solutions to
% of n, k that are compatible with the experimental data. The output is
% given as figures appearing on the screen. It is also stored in videofile
% "optcontour" but the frame rate of that format is far too high.

%Read optical data in form of a space-separated LTR (lambda [m], transmittance [0-1],
%and reflectance [0-1]) - file

[filenameLTR, pathnameLTR] = uigetfile('*.ltr', 'Select a LTR-file containing three columns, lambda, T and R for your sample');
LTR=dlmread(fullfile(pathnameLTR,filenameLTR));

lambda=LTR(:,1);
T_m=LTR(:,2);
R_m=LTR(:,3);

% Set parameters of unknown layer, film thickness (dfilm), coherence(coher),
% incident angle (phi), polarization (pol)

prompt      =   {'film thickness [nm]','incident angle [degrees]', 'coherence (1 for coherent, 0 for non-coherent)'};
dlg_title   =   'Set parameters';
num_lines   =   1;
def         =   {'50','0','1'};
answer      =   inputdlg(prompt,dlg_title,num_lines,def);

d           =   str2num(cell2mat(answer(1)));
phi         =   str2num(cell2mat(answer(2)));
coher(1)    =   str2num(cell2mat(answer(3)));

% Set number of substrate layers.

prompt      =   {'Number of substrate layers'};
dlg_title   =   'Set number of substrate layers';
num_lines   =   1;
def         =   {'1'};
sub_layers  =   inputdlg(prompt,dlg_title,num_lines,def);
num_sub_layers  =   str2num(cell2mat(sub_layers(1)));

% Enter substrate optical constants, thickness (dsub) and coherence.
% Presently the file must have values at the same wavelengths as the
% experimental data. The problem with interpolating is that the wavelengths
% go from large to small.

if (num_sub_layers>0)
    for s=1:num_sub_layers

        [filenameLnksub, pathnameLnksub] = uigetfile('*.nk', ['Select a Lnk-file for substrate layer ' num2str(s)]);
        Lnksub=dlmread(fullfile(pathnameLnksub,filenameLnksub));
       % for i=2:3	% Interpolate to lambda steplength, does not work here.
       % subint(:,i)=interp1(Lnksub(:,1),Lnksub(:,i),lambda(:));  
       % end
       % nsub(:,s)=subint(:,2);
       % ksub(:,s)=subint(:,3);
        nsub(:,s)=Lnksub(:,2);
        ksub(:,s)=Lnksub(:,3);
        Nsub(s,:)=nsub(:,s)'-i*ksub(:,s)';
        prompt      =   {'Layer thickness [nm]','Coherence (1 for coherent, 0 for non-coherent)'};
        dlg_title   =   'Set layer thickness and coherence';
        num_lines   =   1;
        def         =   {'1e6','0'};
        
        dsub_coher  =   inputdlg(prompt,dlg_title,num_lines,def);
        dsub(s)     =   str2num(cell2mat(dsub_coher(1)));
        coher(s+1)    =   str2num(cell2mat(dsub_coher(2)));
        
    end
end

lb_film1=[0 0];
ub_film1=[10 10];

options=optimset('Display','off','TolFun',1e-10);
j=0;

% The calculkations take very long if performed for all wavelengths.
% Therefore it is advisable to select one wavelength or a smaller number.
for L=1:40:max(size(lambda)) % This chooses as an example every 40th wavelength.
% for L=400
    % You may also select a specific wavelength of interest 
    % by choosing a certain index L (in the line above L=400)
    % corresponding to the Lth wavelength in the input file.
    
    lambda(L)
    j=j+1;
    dn1=(ub_film1(1)-lb_film1(1))/100;
    dk1=(ub_film1(2)-lb_film1(2))/100;

    for y=1:100
        for x=1:100
            n1(x,y)=lb_film1(1)+x*dn1;
            k1(x,y)=lb_film1(2)+y*dk1;
            [R_t,T_t,R_p,R_s,T_p,T_s,theta_out, theta,S_p,S_s]=RT_calc(phi,lambda(L),[n1(x,y)-i*k1(x,y);Nsub(:,L)],[d dsub],coher,1,1);
            R_c(x,y)=R_t;
            T_c(x,y)=T_t;
            felyta1(x,y)=log10((R_c(x,y)-R_m(L)).^2+(T_c(x,y)-T_m(L)).^2);
        end
    end
    
    
    [C,h]=contour(n1,k1,T_c,[T_m(L),T_m(L)]);
    hold on
    [C,h]=contour(n1,k1,R_c,[R_m(L),R_m(L)]);
    title(['d=' num2str(d) 'nm   lambda=' num2str(lambda(L)) 'nm']);
    hold off
    F(j) = getframe;

end

% The movie is played regrettably with the last thickness and wavelength as
% a heading to all the images. I do not know how to view the stored file.
     movie(F,2,0.5)
     v=VideoWriter('optcontour','archival');
     open(v)
     writeVideo(v,F)
     close(v)