%% Simple test script for Passler & Paarmann transfer matrix code
        %   calculates Otto geometry reflectivity spectra and field distributions 
        %   for 4 material systems considered in the manuscript
        %   materials limited to the database defined in the function 
        %   "passler_transfer_matrix.m"
        %
        %   Copyright Nikolai Christian Passler and Alexander Paarmann 2017
        %       email: alexander.paarmann@fhi-berlin.mpg.de
        %   Please cite the associated publication is you use this code

clear 
close all

% set up frequency axis w, air gaps d, incidence angle theta
w = 350:1400;        % unit: 1/cm
d = [2,3.5,5.5,7.5]; % unit: um
theta = 30/180*pi;   
% set up parameters for field distributions
dz   = 2e-8; % set to 0 for no calculations of the field distribution
Emax = 5;    % plot range maximum for field distributions
kpick = 2;   % gap index for field distributions

% set up material systems, leave air gap on NaN here
mats{1} = {'KRS5','vac','SiC6Hc'};
ds{1}   = {0     ,NaN , 0};
mats{2} = {'KRS5','vac','GaNc','SiC6Hc'};
ds{2}   = {0     ,NaN , 2e-6  ,0};
mats{3} = {'KRS5','vac','SiC6Hc','GaNc','SiC6Hc'};
ds{3}   = {0     ,NaN , .7e-6    ,2e-6  ,0};
mats{4} = {'KRS5','vac','SiO2a'};
ds{4}   = {0     ,NaN , 0};

for ks = 1:length(mats), Nds{ks} = length(ds{ks}); end

titles  = {'SiC','GaN/SiC','SiC/GaN/SiC','a-SiO_2'};
legends = {['d_{GAP} =',num2str(d(1)),'\mum'],['d_{GAP} =',num2str(d(2)),'\mum'],['d_{GAP} =',num2str(d(3)),'\mum'],['d_{GAP} =',num2str(d(4)),'\mum']};
f1 = figure(1);
f2 = figure(2);

% loop through air gaps and material systems, and plot spectra and field distributions
for kd = 1:length(d)
    for ks = 1:length(mats)
        ds{ks}{2} = d(kd)*1e-6; % set specific air gap here
        if kd == kpick && dz > 0
            [R,T,Ed,zd] = passler_transfer_matrix(w,theta,ds{ks},mats{ks},dz);
        else
            [R,T,Ed,zd] = passler_transfer_matrix(w,theta,ds{ks},mats{ks},0);
        end
        figure(1), p1(ks)=subplot(2,2,ks);
        if kd == 1, hold off, else hold on, end 
        plot(w,R(1,:));
        xlim([min(w),max(w)]); xlabel('Frequency (1/cm)'); ylabel('Reflectivity');
        ztot = max(zd{end});
        if kd == kpick && dz > 0  % Plot field distributions for a specific air gap
            figure(2+ks),
            for kl = 2:length(ds{ks})
                p2{ks}(kl) = subplot(Nds{ks}-1,1,kl-1);
                imagesc(w,zd{kl},permute(abs(Ed{kl}(1,:,:)-Ed{kl}(4,:,:)),[3,2,1]),[0,Emax])                
            end
            % make it look pretty
            for kl = 2:length(ds{ks}) 
                z0 = min(zd{kl}); z1 = max(zd{kl});
                set(p2{ks}(kl),'Position',[0.15,0.95-z1/ztot*0.8,0.8,(z1-z0)/ztot*0.8]);
                p2{ks}(kl); p2{ks}(kl); text(min(w)+(max(w)-min(w))*.05,z0+(ztot)*.05,mats{ks}{kl},'Color','w');
                if kl<length(ds{ks})
                    set(p2{ks}(kl),'XTickLabel',''); 
                else
                    p2{ks}(kl); xlabel('Frequency (1/cm)'); ylabel('z (m)'); 
                end
            end                
            colormap('Jet');
        end
    end
end
% make it look pretty 
for ks = 1:length(mats)
    p1(ks)=subplot(2,2,ks);
    title(titles{ks})
    legend(legends,'Location','southeast');
    set(p1(ks),'Box','on');
end
set(f1,'Position',[50,50,1200,800]);