%**************************************************************************
%**************************************************************************
% Clayton Garrett Webster
% M.Sc Project 

%PURPOSE:
%The  purpose of this project was to develop a computational algorithm for 
%solving the stationary Maxwell equations at the interface between two planar 
%waveguides. The algorithm had be general enough to operate with arbitrary 
%refractive index profiles on both sides of the interface and with arbitrary 
%incident electric field from one side. The algorithm computes the reflection 
%and transmission coefficients for trapped waveguide modes, radiation waves 
%and edge states localized along the interface.
%
%ASSUMPTION:  x is real
% 
%**************************************************************************
%**************************************************************************
           
   clc
   clear all;
   counter=0;
  
   for a_2=1:0.1:4 
       counter=counter+1;
       cwd = pwd;
            cd(tempdir);
            pack
            cd(cwd)
   
   L     = 20;      
   %periodic boundary condition (2L is the length along the x-axis)
   N     = 1000;      
   %the number of subintervals from the origin to L (x-axis)
   M     = 2*N;     
   %the total number of subintervals from -L to L   
   h     = L/N;	    
   %distance between the mesh points (x-axis)
   a_1   = 1.0;     
   %interior of unit step function
   %a_2   = 2.0;     
   %interior of unit step function
   Q_1   = 1.0;     
   %value of the step fuunction q(x) where x<=-a_1 and x>=a_1
   Q_2   = 2.0;     
   %value of the step fuunction q(x) where -a_1<=x<=a_1
   Q_3   = 1.0;     
   %value of the step fuunction q(x) where x<=-a_2 and x>=a_2
   Q_4   = 2.0;     
   %value of the step fuunction q(x) where -a_2<=x<=a_2
   tol   = 0;
   
   A_1   = zeros (M,M); 
   %A_1 is a MxM matrix of zeros
   %A_1 becomes the matrix defined using finite-difference to the left 
   %of the interface
   A_1(1,M)= 1/h^2;
   A_1(1,2)= 1/h^2;
   A_1(M,1)= 1/h^2;
   A_1(M,M-1)= 1/h^2;
   
   A_2   = zeros (M,M); 
   %A_2 is a MxM matrix of zeros
   %A_2 becomes the matrix defined using finite-difference to the right 
   %of the interface
   A_2(1,M)= 1/h^2;
   A_2(1,2)= 1/h^2;
   A_2(M,1)= 1/h^2;
   A_2(M,M-1)= 1/h^2;
  
   c_0 = 1;             
   %vector c coeffecient
   
   %***********************************************************************
   %CONSTUCTION OF MATRIX A_1 to the LEFT of the INTERFACE
   %***********************************************************************
   x = -L+h:h:L;

   x_0=-L;
   %Set-up the subinterval x_n that discretizes the  matrix A_1 
   for n = 1:1:M
        x_n=-L+n*h;
        if x_n<-a_1 | x_n>a_1
        q_n=Q_1;  
        else
        q_n=Q_2;
        end
    A_1(n,n)= -2/h^2 + q_n;     
    %The diagonal elements of the matrix A_1
   end

   for n=2:1:M-1
    A_1(n,n-1)= 1/h^2;
    A_1(n,n+1)= 1/h^2;
    %Finish the constuction of the matrix A_1
   end

   [V_1,D_1] = eig(A_1);      
   %V_1 is the FUNDAMENTAL MATRIX of eigenvectors of A_1 
   %D_1 is the diagonal matrix of eigenvalues of A_1

%***********************************************************************
%CONSTUCTION OF MATRIX A_2 to the RIGHT of the INTERFACE
%***********************************************************************

   %Set-up the subinterval x_t that discretizes the  matrix A_2
   for t = 1:1:M
        x_t=-L+t*h;
        if x_t<-a_2 | x_t>a_2
        q_t=Q_3;  
        else
        q_t=Q_4;
        end
   A_2(t,t)= -2/h^2 + q_t;   
   %The diagonal elements of the matrix A_2
    
   end

   for t=2:1:M-1
    A_2(t,t-1)= 1/h^2;
    A_2(t,t+1)= 1/h^2;
    %Finish the constuction of the matrix A_2
   end

%**************************************************************************
%INTERMEDIATE MATRIX CALCULATIONS
%**************************************************************************
   
   [V_2,D_2] = eig(A_2);       
   %V_2 is the FUNDAMENTAL MATRIX of eigenvectors of A_2
   %D_2 is the diagonal matrix of eigenvalues of A_2

   D_1_root = sqrt(D_1);       
   %take square root of eigenvalues of A_1
   D_2_root = sqrt(D_2);       
   %take square root of eigenvalues of A_2

   eigenvalue1 = diag(D_1);    
   %column vector containing the eigenvalues of A_1
   eigenvalue2 = diag(D_2);    
   %column vector containing the eigenvalues of A_2

   index1 = find(eigenvalue1>0);
   for k=1:1:length(index1)
       eigenvalue_1(k) = eigenvalue1(index1(k));
   end
   Eigenvalue1=transpose(eigenvalue_1);
   %Contains all eigenvalue1 that are positive
  
   index2 = find(eigenvalue2>0);
   for k=1:1:length(index2)
       eigenvalue_2(k) = eigenvalue2(index2(k));
   end
   Eigenvalue2=transpose(eigenvalue_2);
   %Contains all eigenvalue2 that are positive
   
   Grammian_21 = V_2'*V_1;        
   %Grammian matrix of inner products of eigenvectors of both systems A_1 and A_2
   Grammian_11 = V_1'*V_1;        
   %Grammian matrix of inner products of eigenvectors of the system A_1 
   Grammian_22 = V_2'*V_2;       
   %Grammian matrix of inner products of eigenvectors of the system A_2 
   
   %NOTE:  V_2' and V_1' are the transpose of V_2 and V_1 respectively

%**************************************************************************
%SOLVEING THE MATRIX SYSTEM Ba=c for the vectors a and b
%**************************************************************************

   B=D_2_root*Grammian_21 + Grammian_21*D_1_root;
   C = -c_0*(D_2_root*Grammian_21(:,M)-D_1_root(M,M)*Grammian_21(:,M));
   
   %BOTH matrices B and C are found using Orthogonal Projection Operators
   a = B\C;                      
   % a=inv(B)*C

   c = zeros(M,1); c(M) = c_0;     
   %M dimensional vector whose Mth component is c_0 
   
   %NOTE:  Since the eigenvalues are sorted in ascending order the vector 
   %c whose mth component corresponds to the largest eigenvalue
   
   b = Grammian_21*a  + Grammian_21*c;
   %construction of transmitted coefficients
   
   for k=1:1:length(index1)
       a_1(k) = a(index1(k));
   end
   a1=transpose(a_1);
   
   for k=1:1:length(index1)
       c_1(k) = c(index1(k));
   end
   c1=transpose(c_1);
   
   for k=1:1:length(index2)
       b_1(k) = b(index2(k));
   end
   b1=transpose(b_1);
   
   
   %a1 and b1 are the reflected and transmitted wave coefficients 
   %corresponding to the nonnegative eigenvalues


%**************************************************************************
%CONSERVATION OF ENERGY RELATION
%**************************************************************************

   I_reflected= 2*sum(sqrt(diag(Eigenvalue1))*abs(a1).^2);    
   % total reflected intensity (energy)
   I_incident = 2*sum(sqrt(diag(Eigenvalue1))*abs(c1).^2);    
   % total incident intensity (energy)
   I_transmitted = 2*sum(sqrt(diag(Eigenvalue2))*abs(b1).^2);  
   % total transmitted intensity (energy)
   
   R=I_reflected/I_incident;
   T=I_transmitted/I_incident;
   answer = R+T;
   
   R_vector(counter)=R;
   T_vector(counter)=T;
   tau(counter)=a_2;
   
end
figure(1);   
fig1=plot(tau,R_vector)
set(fig1,'Linewidth',2)
set(fig1,{'Color'},{'r'})

xlabel('\rho^{ +}', 'Fontsize', 12)
ylabel('Reflectance (R)', 'Fontsize', 12)
%legend(h,'R','T')
title('Reflectance versus \rho^{+}', 'Fontsize', 12)

figure(2);   
fig2=plot(tau, T_vector)
set(fig2,'Linewidth',2)
set(fig2,{'Color'},{'b'})

xlabel('\rho^{ +}', 'Fontsize', 12)
ylabel('Transmittance (T)', 'Fontsize', 12)
%legend(h,'R','T')
title('Transmittance versus \rho^{+}', 'Fontsize', 12)










   
