%**************************************************************************
%**************************************************************************
% 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 complex
% 
%**************************************************************************
%**************************************************************************

   clc
   clear all;
   L     = 30;      
   %2L is the length along the zeta-axis
   LZ    = 30;      
   %2LZ is the length along the z-axis
   N     = 100;      
   %the number of subintervals from the origin to L (zeta-axis)
   P     = 100;      
   %the number of subintervals from the origin to LZ (z-axis)
   Q     = 2*P;     
   %the total number of subintervals from -LZ to LZ
   M     = 2*N-1;    
   %the total number of subintervals from -LX to LX                
   h     = L/N;	    
   %distance between the mesh points (zeta-axis)
   g     = LZ/P;    
   %distance between the mesh points (z-axis)
   a_1   = 1.0;     
   %interior of unit step function (left of interface)
   a_2   = 2.0;    
   %interior of unit step function (right of interface)
   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   = 4.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_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
   
   c_0 = 1;             
   %vector c coeffecient
%***********************************************************************
%ABSORBING LAYER IN "COMPLEX TRANSFORMATION"
%***********************************************************************
   
   delta  = 50*h;
   %thickness of the layer (ie. 45 layers)
   Delta_0 = 1;
   Delta = -Delta_0/delta;
   %LINEAR ABSORBING LAYER
   
   %Set-up the vector containing the mesh points zeta_n for matrix A_1

%***********************************************************************
%CONSTUCTION OF MATRIX A_1 to the LEFT of the INTERFACE
%***********************************************************************
   
   % NOTE NON-SYMMETIC INTERVAL: zeta = -L+h:h:L-h; OR zeta = -L+h:h:L;
   %Set-up the subinterval zeta_n that discretizes the  matrix A_1 
   zeta = -L+h:h:L-h;
               
   for n = 1:1:M
   
        if zeta(n)<-a_1 | zeta(n)>a_1
        q_n = Q_1;  
        else
        q_n = Q_2;
        end
    
        if (zeta(n)<=(L-delta)) & (zeta(n)>=(-L+delta))
        C_n = 1; 
        else
        C_n = 1-i*Delta;
        end
        %C_n is the discretization of the derivative of the absorbing layer
   
   A_1(n,n)= -2/(C_n*h^2) + q_n;     
   %diagonal elements of A_1
        
        if n == 1
        A_1(1,2) = 1/(C_n*h^2);
        elseif n == M
        A_1(M,M-1) = 1/(C_n*h^2);
        else
        A_1(n,n-1)= 1/(C_n*h^2);
        A_1(n,n+1)= 1/(C_n*h^2);
        %Finish the constuction of A_1
        end
        
    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
        if zeta(t)<-a_2 | zeta(t)>a_2
        q_t=Q_3;  
        else
        q_t=Q_4;
        end

        if (zeta(t)<=(L-delta)) & (zeta(t)>=(-L+delta))
        C_t = 1; 
        else
        C_t = 1-i*Delta;
        end
    
   A_2(t,t)= -2/(C_t*h^2) + q_t;
   %Diagonal elements of A_2
   
        if t == 1
        A_2(1,2) = 1/(C_t*h^2);
        elseif t == M
        A_2(M,M-1) = 1/(C_t*h^2);
        else
        A_2(t,t-1)= 1/(C_t*h^2);
        A_2(t,t+1)= 1/(C_t*h^2);
        %Finish the constuction of the matrix A_2
        end
        
   end

   [V_2,D_2] = eig(A_2);     
   abs_norm = norm((A_2*V_2) - (D_2*V_2)); 
   %V_2 is the FUNDAMENTAL MATRIX of eigenvectors of A_2
   %D_2 is the diagonal matrix of eigenvalues of A_2
   
%**************************************************************************
%INTERMEDIATE MATRIX CALCULATIONS:  SORT EIGENVALUES and CORRESPONDING
%MATRICES ACCORDING TO THE EIGENVALUE SPECTRUM
%**************************************************************************

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

   eigenvalue1_real = real(eigenvalue1);         
   %vector containing the real part of the vector of eigenvalues for A_1
   eigenvalue1_imaginary = imag(eigenvalue1);    
   %vector containing the imaginary part of the vector of eigenvalues for A_1

   eigenvalue2_real = real(eigenvalue2);         
   %vector containing the real part of the vector of eigenvalues for A_2
   eigenvalue2_imaginary = imag(eigenvalue2);    
   %vector containing the imaginary part of the vector of eigenvalues for A_2

   [x1_real, index_1] = sort(eigenvalue1_real);
   %sort the real part of the eigenvalues of A_1 in ascending order
   %index_1 holds the order of the eigenvalues of A_1
   [x2_real, index_2] = sort(eigenvalue2_real);
   %sort the real part of the eigenvalues of A_2 in ascending order
   %index_2 holds the order of the eigenvalues of A_2

   x1_imaginary = eigenvalue1_imaginary(index_1);
   x2_imaginary = eigenvalue2_imaginary(index_2);
   
   
   V_1_sort = V_1(:,index_1);
   %Re-sorted Fundamental matrix of A_1 according to the sorted eigenvalues
   V_2_sort = V_2(:,index_2);
   %Re-sorted Fundamental matrix of A_2 according to the sorted eigenvalues

   D_1_sort = D_1(:,index_1);
   %Re-sorted diagnonal matrix of the eigenvalues A_1 
   %according to the sorted eigenvalues
   D_2_sort = D_2(:,index_2);
   %Re-sorted diagnonal matrix of the eigenvalues A_2 
   %according to the sorted eigenvalues

   for p=1:1:M
        eigenvalue1_sort(p) = D_1_sort(find(D_1_sort(:,p)), p);
        eigenvalue2_sort(p) = D_2_sort(find(D_2_sort(:,p)), p);
   end
   %short loop to construct row vectors of sorted eigenvalues

   EIGENVALUES1 = transpose(eigenvalue1_sort);
   EIGENVALUES2 = transpose(eigenvalue2_sort);
   %these column vectors are now the vectors containing the sorted
   %eigenvalues described above

   D_1_root = zeros(M,M);
   D_2_root = zeros(M,M);
   %Construction of two matrices whose diagonal elements are the square
   %root of the eigenvalues of A_1 and A_2 resp.

   for s=1:1:M
        D_1_diag(s,s)=EIGENVALUES1(s);
        D_2_diag(s,s)=EIGENVALUES2(s);
   end

   D_1_root = sqrt(D_1_diag);       
   D_2_root = sqrt(D_2_diag);       
   %Complete the construction of the two matrices whose diagonal elements 
   %are the square root of the eigenvalues of A_1 and A_2 resp.

   R_EIG1 = real(EIGENVALUES1);
   R_EIG2 = real(EIGENVALUES2);
   
   index1 = find(R_EIG1>0);
   for k=1:1:length(index1)
       r_eigenvalue1(k) = R_EIG1(index1(k));
   end
   R_Eigenvalue1=transpose(r_eigenvalue1);
   %Contains all eigenvalue1 that are real and positive (sorted)
  
   index2 = find(R_EIG2>0);
   for k=1:1:length(index2)
       r_eigenvalue2(k) = R_EIG2(index2(k));
   end
   R_Eigenvalue2=transpose(r_eigenvalue2);
   %Contains all eigenvalue1 that are real and positive (sorted)
  
   
%**************************************************************************
%CONSTRUCT THE MATRIX A THAT SOLVES FOR THE COEFFICIENT VECTORS a and b 
%SIMULTANEOUSLY, using the equation Ax=y
%**************************************************************************

   V1 = V_1_sort*D_1_root;  
   V2 = V_2_sort*D_2_root;
   %please refer to thesis (Chapter 3) for the reasoning behind the
   %construction of such matrices

   A = [-V_1_sort, V_2_sort; V1, V2];
   %the matrix A stores all info needed to solve for vectors a and b
   y = c_0*[V_1_sort(:,M); D_1_root(M,M)*V_1_sort(:,M)];
   %right-side of the matrix equation

   x = A\y;                      
   %x=inv(A)*y
   %the vector x contains all reflected and transmitted coefficients 

   for r=1:1:M
   a(r) = x(r);
   end
   a_ref = transpose(a);
   %extract the vector a of reflected coefficients
   
   for r=M+1:1:2*M
   b(r-M) = x(r);
   end
   b_trans = transpose(b);
   %extract the vector b of transmitted coefficients

   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 
   
   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);
   

%**************************************************************************
%FIGURES
%**************************************************************************

   %Energy Spectrum for Tranmitted Waves
   figure(1); plot (R_EIG2, (2*sqrt(abs(EIGENVALUES2))).*(abs(b_trans).^2),  '*'); 
   xlabel('Re(\lambda^+)', 'FontSize', 12);
   ylabel('Energy Density P_{ j }^{  + }', 'FontSize', 12);
   title('\it{Energy Spectrum for Transmitted Radiation and Edge Modes}', 'FontSize', 12);

   %Energy Spectrum for Reflected Waves
   figure(2); plot (R_EIG1, (2*sqrt(abs(EIGENVALUES1))).*(abs(a_ref).^2), '*');
   xlabel('Re(\lambda^-)', 'FontSize', 12);
   ylabel('Energy Density P_{ j }^{  - }', 'FontSize', 12);
   title('\it{Energy Spectrum for Reflected Radiation and Edge Modes}', 'FontSize', 12);

   figure(3); plot (x1_real, x1_imaginary, '*', 'Color', 'r');
   xlabel('Re(\lambda^-)', 'FontSize', 12);
   ylabel('Im(\lambda^-)', 'FontSize', 12);
   title('\it{Location of Complex Eigenvalues \lambda^-}', 'FontSize', 12);

   figure(4); plot (x2_real, x2_imaginary, '*', 'Color', 'r');
   xlabel('Re(\lambda^+)', 'FontSize', 12);
   ylabel('Im(\lambda^+)', 'FontSize', 12);
   title('\it{Location of Complex Eigenvalues \lambda^+}', 'FontSize', 12);
   
   
   figure(7);plot(zeta, real(V_1_sort(:, 1)));
   xlabel('\xi', 'FontSize', 12);
   figure(8);plot(zeta, real(V_1_sort(:, 3)));
   xlabel('\xi', 'FontSize', 12);
   figure(9);plot(zeta, real(V_1_sort(:, 6)));
   xlabel('\xi', 'FontSize', 12);
   figure(10);plot(zeta, real(V_1_sort(:, 7)));
   xlabel('\xi', 'FontSize', 12);
   figure(11);plot(zeta, real(V_1_sort(:, 9)));
   xlabel('\xi', 'FontSize', 12);
   figure(12);plot(zeta, real(V_1_sort(:, 11)));
   
%**************************************************************************
%CONSERVATION OF ENERGY RELATION
%**************************************************************************

   I_reflected= 2*sum(sqrt(diag(R_Eigenvalue1))*abs(a1).^2)    
   % total reflected intensity (energy)
   I_incident = 2*sum(sqrt(diag(R_Eigenvalue1))*abs(c1).^2)    
   % total incident intensity (energy)
   I_transmitted = 2*sum(sqrt(diag(R_Eigenvalue2))*abs(b1).^2)  
   % total transmitted intensity (energy)
   
   R=I_reflected/I_incident
   T=I_transmitted/I_incident
   answer = R+T

   
%**************************************************************************
%CONSTRUCT THE 3-D PLOT OF THE WAVE-GUIDE
%**************************************************************************


   z = -LZ : g : LZ;
   %Set-up the subinterval z_t that discretizes the z-axis

   %Set-up the solution matrices to the left and right of the interface
   Psi1 = zeros(M,P+1);
   %Solution matrix to the left of the interface
   Psi2 = zeros(M,P);
   %Solution matrix to the right of the interface

   %Solve for the matrix Psi1 (LEFT of the INTERFACE) that solves our 
   %problem containing the reflected and incident waves
   for k = 1 : (P+1)
        E_1 = diag(exp(-i*diag(D_1_root)*z(k)));                        
        %matrix of exponentials 
        Sol_Matrix_1 = V_1_sort*(diag(a)*E_1 + diag(c)*inv(E_1));       
        %matrix containing the coefficients and the basis vectors  
        Psi1(:,k) = sum(Sol_Matrix_1')';            
        %Sums up the columns of transpose(Sol_Matrix_2) and saves them 
        %as a row vector    
   end

   %Solve for the matrix Psi1 (RIGHT of the INTERFACE) that solves our 
   %problem containing the transmitted waves
   for k = (P+2) : Q+1
        E_2 = diag(exp(i*diag(D_2_root)*z(k)));         
        %matrix of exponentials 
        Sol_Matrix_2 = V_2_sort*diag(b)*E_2;                 
        %matrix containing the coefficients and the basis vectors
        Psi2(:,k-P-1) = sum(Sol_Matrix_2')';
        %Sums up the columns of transpose(Sol_Matrix_2) and saves them 
        %as a row vector
   end

%**************************************************************************
%Now define the complete solution from -M to M as a mesh of grid points
%**************************************************************************

   [Z,X] = meshgrid(z,zeta);
   Psi = [Psi1, Psi2];

   %set(figure(3),'Renderer','ZBuffer');
   %figure(3); mesh(Z,X,abs(Psi).^2);
   %xlabel('zeta', 'FontSize', 12);
   %ylabel('z', 'FontSize', 12);
   %zlabel('Wave Solution', 'FontSize', 12);
   %title('\it{Numerical Representation of Waveguide}', 'FontSize', 12);
   
   figure(5); surf(Z,X,abs(Psi).^2);
   %figure(3); surf(Z,X,abs(Psi).^2);
   colorbar;
   view(0,75);
   shading interp;
   lightangle(-45,30);
   set(findobj(gca, 'type', 'surface'),...
        'FaceLighting', 'phong',...
        'AmbientStrength',.3,'DiffuseStrength',.8,...
        'SpecularStrength',.9,'SpecularExponent',25,...
        'BackFaceLighting','unlit');
   xlabel('z', 'FontSize', 12);
   ylabel('\xi', 'FontSize', 12);
   zlabel(' Amplitude of Wave Solution', 'FontSize', 12);
   title('\it{Numerical Representation of Waveguide}', 'FontSize', 12);
   
   figure(6); surf(Z,X,abs(Psi).^2);
   %figure(3); surf(Z,X,abs(Psi).^2);
   %colorbar;
   view(0,75);
   shading interp;
   lightangle(-45,30);
   set(findobj(gca, 'type', 'surface'),...
        'FaceLighting', 'phong',...
        'AmbientStrength',.3,'DiffuseStrength',.8,...
        'SpecularStrength',.9,'SpecularExponent',25,...
        'BackFaceLighting', 'unlit');
   xlabel('z', 'FontSize', 12);
   ylabel('\xi', 'FontSize', 12);
   zlabel(' Amplitude of Wave Solution', 'FontSize', 12);
   title('\it{Numerical Representation of Radiation from Waveguide}', 'FontSize', 12);

   
%**************************************************************************
%THE END
%**************************************************************************


