Code


%
% TM_Optical_v4
% 20/2/2019
% J. Allam
% Generalised OPTICAL transmission through multiple layers
% with peak detection and calculation of eigenfunctions
% based on Levi Exercise 4.2
% converted to scan wavelength not energy
 
clear
clf
 
%change the following file locations:
input_file='TM_O_layers.dat';
output_transmission='TM_O_trans.txt';
output_peaks='TM_O_peaks.txt';
output_waves='TM_O_waves.txt';
output_DiffPhase='TM_O_DiffPhase.txt';
 
% generalised input from file for arbitrary number of layers
load (input_file);
 
layer = TM_O_layers(:,1);
L = TM_O_layers(:,2)*1e-9;  %distance array (nm)
 
 
 
N=max(layer);           %number of layers
 
lamda_min=100 ;               %minimum wavelength (nm)
lamda_max=4000;               %maximum wavelength (nm)      
npoints=1000;                 %number of points in wavelength scan
 
% constants and parameters
dlamda=(lamda_max-lamda_min)/npoints;       %wavelength increment (nm)
eye=complex(0.,1.);     %square root of -1
c0=2.998e8;             %speed of light
 
%refractive index dependent on lambda
for i=1:N
for jIndex=1:npoints
j = lamda_min+i*dlamda;
if mod(i,2)== 1  
    mu(i,jIndex) = sqrt(1+(0.6961663*j^2)/(j^2 - 0.0684043^2) +(0.4079426*j^2)/(j^2 - 0.1162414^2) + (0.8974794*j^2)/(j^2 - 9.896161^2));      %refractive index
    mumax=max(mu(i,jIndex));
    mumin=min(mu(i,jIndex));
else
    mu(i,jIndex)= sqrt(5.913+(0.2441/((j^2)-0.0803)));
    mumax=max(mu(i,jIndex));
    mumin=min(mu(i,jIndex));
end
end
end       %refractive index
 
%write input data to figure as a table
f = figure(1);
set (f,'Position',[200 200 400 300]);
cnames = {'Layer','L(nm)','Index'};
t = uitable('Parent',f,'Data',TM_O_layers,'ColumnName',cnames,'Position',[20 20 360 250]);
 
%main calculation of transmission parameters against energy
   for j=1:npoints
      lamda(j)=dlamda*j+lamda_min+1e-5;   %add (pi*1.0e-5) to avoid divide by zero
      !energy(j)=1240/lamda(j)
        if j<npoints
          lamdaminusone(j)=lamda(j);    %to make an array with the same dimensions as the derivative
        end
        bigP=[1,0;0,1];         %default value of matrix bigP
       for i=1:N
         k(i)=2*pi*mu(i,j)/lamda(j)/1e-9;
      end
      for n=1:(N-1)
         fac = k(n+1)/k(n);    
         p(1,1)=0.5*(1+fac)*exp(-eye*k(n)*L(n));
         p(1,2)=0.5*(1-fac)*exp(-eye*k(n)*L(n));
         p(2,1)=0.5*(1-fac)*exp(eye*k(n)*L(n));
         p(2,2)=0.5*(1+fac)*exp(eye*k(n)*L(n));
         bigP=bigP*p;
      end
      Trans(j)=(abs(1/bigP(1,1)))^2;                    %transmission probability
      TransPhase(j)=angle(1/bigP(1,1))*180/pi;          %transmitted phase
      Ref(j)=(abs(bigP(2,1)/bigP(1,1)))^2;              %reflection pro ability
      RefPhase(j)=angle(bigP(2,1)/bigP(1,1))*180/pi;    %reflected phase
   end
Tmin=min(Trans);    % min transmission for graph plotting
DPhase=diff(TransPhase);
DiffPhase=(DPhase/dlamda)*lamda(j)/2/pi/c0; %convert to angular frequency CHECK THIS
MaxDiffPhase=max(DiffPhase);
MinDiffPhase=min(abs(DiffPhase));
   
%plot potential, transmission and reflection coefficients
figure(2);              
% generalised generation of potential-distance for
% arbitrary number of input layers
mup=[mu(:,1)';mu(:,1)'];mup=mup(:);
dx=1e-12;               %small distance increment used in potential plot
Lx(1)=L(1);
x1=L;
x2=L;
x1(1)=0;
x2(1)=Lx(1)-dx;
for i=2:N
   for j=2:i
      Lx(i)=L(j)+Lx(j-1);               %distance, x
   end
   x1(j)= Lx(j-1);
   x2(j) = Lx(j)- dx;
end
x3=[x1';x2'];
xp=x3(:)*1e9;
maxL=Lx(N)*1e9;

fprintf('size(A) is [%s]\n', int2str(size(mup)))
fprintf('size(B) is [%s]\n', int2str(size(xp)))

%subplot(3,3,1),plot(xp,mup),axis([0,Lx(N)*1e9,0.9*mumin,1.2*mumax]);
subplot(3,3,1),plot(xp,mup),axis([0,Lx(N)*1e9,0.9*mumin,1.2*mumax]);
xlabel('Position, x (nm)'),ylabel('Refractive index');
 
subplot(3,3,2),plot(Trans,lamda),axis([0,2,lamda_min,lamda_max]);
xlabel('Transmission coefficient'),ylabel('lamda (nm)');
 
subplot(3,3,3),plot(log10(Trans),lamda),axis([log10(Tmin),0,lamda_min,lamda_max]);
%subplot(3,3,3),plot((Trans),lamda),axis([1e-5,1,700,900]);
xlabel('log10(trans. coeff.)'),ylabel('lamda (nm)');
 
subplot(3,3,4),plot(TransPhase,lamda),axis([-180,180,lamda_min,lamda_max]);
xlabel('TransPhase'),ylabel('lamda (nm)');
 
%subplot(3,3,5);
%h=plot(log10(DiffPhase),lamdaminusone,'ro');
%set(h,'markersize',2);
%axis([log10(MinDiffPhase),log10(MaxDiffPhase),lamda_min,lamda_max]);
%xlabel('log10(DiffPhase)'),ylabel('Energy, E (eV)');
 
subplot(3,3,6),plot(Ref,lamda),axis([0,1,lamda_min,lamda_max]);
xlabel('ref. coeff.'),ylabel('lamda(nm)');
 
subplot(3,3,7),plot(RefPhase,lamda),axis([-180,180,lamda_min,lamda_max]);
xlabel('RefPhase'),ylabel('lamda(nm)');
 
figure (4);
plot(lamda,Ref),axis([lamda_min,lamda_max,-0.1,1.2]),xlabel('lamda (nm)'),ylabel('ref. coeff.');
 
 
%output trans to file
fileid=fopen(output_transmission,'w');
fprintf(fileid, '%s\t %s\t %s\t %s\t %s\n', 'lamda','Trans','TransPhase','Ref','RefPhase');
fprintf(fileid, '%f %e %f %f %f \r\n', [lamda; Trans; TransPhase; Ref; RefPhase]);
fclose(fileid);
 
 
%output differentiated phase to file
fileid=fopen(output_DiffPhase,'w');
fprintf(fileid, '%s\t %s\n', 'lamda','DiffPhase');
fprintf(fileid, '%f %e \r\n', [lamdaminusone; DiffPhase]);
fclose(fileid);
 
%find peaks IN REFLECTANCE NOT TRANSMISSION
npeak=0;
for j=2:npoints-1
    Ref(j),lamda(j)
    if Ref(j)>Ref(j-1)
        if Ref(j)>Ref(j+1)
            npeak=npeak+1;
            lamdapeak(npeak)=lamda(j);
            Rpeak(npeak)=Ref(j);
        end
    end
end
 
%output peaks to file
fileid=fopen(output_peaks,'w');
fprintf(fileid, '%s\t %s\n', 'lamdapeak','Rpeak');
fprintf(fileid, '%f %f \r\n', [lamdapeak; Rpeak]);
fclose(fileid);
 
%find wavefunctions at peak energies VVV
%prepared for subplots depending on number of peaks detected;
figure(3);
peakplotrows=1;
if npeak>4
        peakplotrows=2;
end
if npeak>6
        peakplotrows=3;
end
if npeak>9
        peakplotrows=4;
end
peakplotcols=round(npeak/peakplotrows+0.5);
 
%first find complex transmission at peaks for initial waveform on RHS
 
%output waves to file
fileid=fopen(output_waves,'w');
fprintf(fileid, '%s\t %s\n', 'position (nm)','wave');
 
for j=1:npeak
    otherlamda(j)=lamdapeak(j);
        bigP=[1,0;0,1]; %default value of matrix bigP
       for i=1:N
         k(i)=2*pi*mu(i,j)/otherlamda(j)/1e-9;
       end
      for n=1:(N-1)
         fac = k(n+1)/k(n);
         p(1,1)=0.5*(1+fac)*exp(-eye*k(n)*L(n));
         p(1,2)=0.5*(1-fac)*exp(-eye*k(n)*L(n));
         p(2,1)=0.5*(1-fac)*exp(eye*k(n)*L(n));
         p(2,2)=0.5*(1+fac)*exp(eye*k(n)*L(n));
         bigP=bigP*p;
      end
 
   % initial waveform on RHS layer N;
        A(N)=bigP(1,1);
        B(N)=0;
        % subsequent waveforms from sucessive matrices;
       for nn=1:(N-1)
           n=N-nn;
           fac = k(n+1)/k(n);
           p(1,1)=0.5*(1+fac)*exp(-eye*k(n)*L(n));
           p(1,2)=0.5*(1-fac)*exp(-eye*k(n)*L(n));
           p(2,1)=0.5*(1-fac)*exp(eye*k(n)*L(n));
           p(2,2)=0.5*(1+fac)*exp(eye*k(n)*L(n));
           A(n)= p(1,1)*A(n+1)+p(1,2)*B(n+1);
           B(n)= p(2,1)*A(n+1)+p(2,2)*B(n+1);
       end
     
layer=1;
xw(1)=0;
wf(1)=0;
xx=0;
for i = 2:(N*100)
    xw(i) = xw(i-1)+maxL/100/N;
    xx=xx+maxL/N/100;
    wf(i) = real(A(layer)*exp(eye*k(layer)*xx*1e-9)+B(layer)*exp(-eye*k(layer)*xx*1e-9));
    if xx > L(layer)*1e9
        layer=layer+1;
        xx=0;
    end
end
 
 
%output successive waves to same file
fprintf(fileid, '%f %f \r\n', [xw; wf]);
 
%plot each wave as subplot
wfmin=min(wf);
wfmax=max(wf);
subplot(peakplotrows,peakplotcols,j),plot(xw,wf),axis([0,Lx(N)*1e9,wfmin,wfmax]);
xlabel('position (nm)'),ylabel('wavefunction');
ttl = sprintf('peak %3.0f of %3.0f, otherlamda=%3.1f nm, Ref=%1.5f',j,npeak,lamdapeak(j),Rpeak(j));
title (ttl);
end  
 
figure (4);
plot(lamda,Ref),axis([lamda_min,lamda_max,-0.1,1.2]),xlabel('lamda (nm)'),ylabel('ref. coeff.');
 
fclose(fileid);
 
!figure (5);
!plot(energy,Ref),axis([1240/lamda_max,1240/lamda_min,-0.1,1.2]),xlabel('E (eV)'),ylabel('ref. coeff.');