function [mtx_A, mtx_B] = mtx_GEP_Maxwell(n_1, vec_k)
% Generate the matrices mtx_A and mtx_B for the generalized eigenvalue
% problem
%                   mtx_A * x = lambda * mtx_B * x
%
% Input
% n_1   : number of partition for the simple cube in x-direction
% vec_k : wave vector 
%
flag_eps = 1;
n_2      = n_1;
n_3      = n_1;
n        = n_1 * n_2 * n_3;
[diagB, inv_diag_D, delta] = ConstructMtxB(n_1, flag_eps, vec_k);

ii      = 1:3*n;
mtx_B   = sparse(ii,ii,diagB); %diag(diagB);
inv_D   = sparse(ii,ii,inv_diag_D); %diag(inv_diag_D);

[mtx_G] = Sparse_mtx_G(n_1,n_2,n_3,delta,vec_k);

mtx_A  = inv_D * mtx_G;
mtx_A  = mtx_A * mtx_A';
mtx_A  = ( mtx_A' + mtx_A ) / 2;

end %mtx_GEP_Maxwell

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function [mtx_G] = Sparse_mtx_G(n_1,n_2,n_3,delta,vec_k)

     n12  = n_1 * n_2;
     n123 = n12 * n_3;
     
     [i_row_1, j_col_1, val_1] = Gamma_1(n_1, n_2, n_3, delta, vec_k);
     [i_row_2, j_col_2, val_2] = Gamma_2(n_1, n_2, n_3, delta, vec_k);
     [i_row_3, j_col_3, val_3] = Gamma_3(n_1,n_2,n_3, vec_k);
     
     leng_1 = length(i_row_1);
     leng_2 = length(i_row_2);
     leng_3 = length(i_row_3);
     leng   = 2 * ( leng_1  + leng_2 + leng_3 );
     i_row  = zeros(leng,1);
     j_col  = zeros(leng,1);
     val    = zeros(leng,1);
     
     val_1  = val_1 / delta(1);
     val_2  = val_2 / delta(2);
     val_3  = val_3 / delta(3);
%    First block rows     
     i_row(1:leng_3,1) = i_row_3;
     j_col(1:leng_3,1) = n123 + j_col_3;
     val  (1:leng_3,1) = val_3;
     
     ii                   = leng_3 + leng_2;
     i_row(leng_3+1:ii,1) = i_row_2;
     j_col(leng_3+1:ii,1) = 2 * n123 + j_col_2;
     val  (leng_3+1:ii,1) = - val_2;
% Second block rows     
     i2               = ii + leng_3;
     i_row(ii+1:i2,1) = n123 + i_row_3;
     j_col(ii+1:i2,1) = j_col_3;
     val  (ii+1:i2,1) = - val_3;
     
     ii               = i2 + leng_1;
     i_row(i2+1:ii,1) = n123 + i_row_1;
     j_col(i2+1:ii,1) = 2 * n123 + j_col_1;
     val  (i2+1:ii,1) = val_1;
% Third block rows     
     i2               = ii + leng_2;
     i_row(ii+1:i2,1) = 2 * n123 + i_row_2;
     j_col(ii+1:i2,1) = j_col_2;
     val  (ii+1:i2,1) = val_2;
     
     ii               = i2 + leng_1;
     i_row(i2+1:ii,1) = 2 * n123 + i_row_1;
     j_col(i2+1:ii,1) = n123 + j_col_1;
     val  (i2+1:ii,1) = - val_1;
     
     mtx_G            = sparse(i_row, j_col, val);
     mtx_G            = mtx_G';
     
 end %Sparse_mtx_G
 
 %%%%%%%%%%%%%%%%%%%%%%
 
function [i_row, j_col, val] = Gamma_1(n_1, n_2, n_3, delta, vec_k)
 
mtx_C_1 = gen_mtx_C_1(n_1,delta,vec_k);
[i_row_tmp, j_col_tmp, val_tmp] = find(mtx_C_1);
leng                            = length(i_row_tmp);
i_row                           = zeros(leng*n_2*n_3,1);
j_col                           = zeros(leng*n_2*n_3,1);
val                             = zeros(leng*n_2*n_3,1);
for ii = 1:n_2*n_3
    jj                    = ( ii - 1 ) * n_1;
    kk                    = ( ii - 1 ) * leng;
    i_row(kk+1:kk+leng,1) = jj + i_row_tmp;
    j_col(kk+1:kk+leng,1) = jj + j_col_tmp;
    val  (kk+1:kk+leng,1) = val_tmp;
end
end %Gamma_1

%%%%%%%%%%%%%%%%%%%%%%
function [i_row, j_col, val] = Gamma_2(n_1, n_2, n_3, delta, vec_k)
 
[i_row_tmp, j_col_tmp, val_tmp] = gen_mtx_C_2(n_1,n_2,delta,vec_k); 
leng                            = length(i_row_tmp);
i_row                           = zeros(leng*n_3,1);
j_col                           = zeros(leng*n_3,1);
val                             = zeros(leng*n_3,1);
for ii = 1:n_3
    jj                    = ( ii - 1 ) * n_1 * n_2;
    kk                    = ( ii - 1 ) * leng;
    i_row(kk+1:kk+leng,1) = jj + i_row_tmp;
    j_col(kk+1:kk+leng,1) = jj + j_col_tmp;
    val  (kk+1:kk+leng,1) = val_tmp;
end
end %Gamma_2

%%%%%%%%%%%%%%%%%

 function [i_row, j_col, val] = Gamma_3(n_1,n_2,n_3, vec_k)

 global length_z
 
     n12   = n_1 * n_2; 
     i_row = zeros(2*n_1*n_2*n_3,1);
     j_col = zeros(2*n_1*n_2*n_3,1);
     val   = zeros(2*n_1*n_2*n_3,1);
     for jj = 1:n_3-1
         ii                   = ( jj - 1 ) * n12;
         kk                   = 2 * ii;
         i_row(kk+1:kk+n12,1) = ii+1:ii+n12;
         j_col(kk+1:kk+n12,1) = ii+1:ii+n12;
         val  (kk+1:kk+n12,1) = ones(n12,1);
         
         i_row(kk+n12+1:kk+2*n12,1) = ii+1:ii+n12;
         j_col(kk+n12+1:kk+2*n12,1) = ii+n12+1:ii+2*n12;
         val  (kk+n12+1:kk+2*n12,1) = -ones(n12,1);
     end
   
     ii                   = ( n_3 - 1 ) * n12;
     kk                   = 2 * ii;
     i_row(kk+1:kk+n12,1) = ii+1:ii+n12;
     j_col(kk+1:kk+n12,1) = ii+1:ii+n12;
     val  (kk+1:kk+n12,1) = ones(n12,1);
         
     i_row(kk+n12+1:kk+2*n12,1) = ii+1:ii+n12;
     j_col(kk+n12+1:kk+2*n12,1) = 1:n12;
     tmp                        = 2 * pi * length_z * vec_k(3) * 1i;
     val  (kk+n12+1:kk+2*n12,1) = -exp(tmp) * ones(n12,1); 
 end %Gamma_3
 
 %%%%%%%%%%%%%%%%
 

 function [mtx_C_1] = gen_mtx_C_1(n_1,delta,vec_k)
%
% Generate matrix K_1^{*}
%
global length_x length_y length_z

     mtx_C_1 = eye(n_1);

     for ii = 1:n_1-1
        mtx_C_1(ii,ii+1) = - 1;
     end

     tmp            = 2 * pi * length_x * vec_k(1) * 1i;
     mtx_C_1(n_1,1) = -exp(tmp);
 end
 
 % ==========================
 

 function [i_row, j_col, val] = gen_mtx_C_2(n_1,n_2,delta,vec_k)
%
% Generate matrix K_2^{*}
%
global length_x length_y length_z

     i_row = zeros(2*n_1*n_2,1);
     j_col = zeros(2*n_1*n_2,1);
     val   = zeros(2*n_1*n_2,1);
     for jj = 1:n_2-1
         ii                   = ( jj - 1 ) * n_1;
         kk                   = 2 * ii;
         i_row(kk+1:kk+n_1,1) = ii+1:ii+n_1;
         j_col(kk+1:kk+n_1,1) = ii+1:ii+n_1;
         val  (kk+1:kk+n_1,1) = ones(n_1,1);
         
         i_row(kk+n_1+1:kk+2*n_1,1) = ii+1:ii+n_1;
         j_col(kk+n_1+1:kk+2*n_1,1) = ii+n_1+1:ii+2*n_1;
         val  (kk+n_1+1:kk+2*n_1,1) = -ones(n_1,1);
     end
   
     ii                   = ( n_2 - 1 ) * n_1;
     kk                   = 2 * ii;
     i_row(kk+1:kk+n_1,1) = ii+1:ii+n_1;
     j_col(kk+1:kk+n_1,1) = ii+1:ii+n_1;
     val  (kk+1:kk+n_1,1) = ones(n_1,1);
         
     i_row(kk+n_1+1:kk+2*n_1,1) = ii+1:ii+n_1;
     j_col(kk+n_1+1:kk+2*n_1,1) = 1:n_1;
     tmp                        = 2 * pi * length_y * vec_k(2) * 1i;
     val  (kk+n_1+1:kk+2*n_1,1) = -exp(tmp) * ones(n_1,1);
         
 end
 
 % ========================
  
function [diagB, inv_diag_D, delta] = ConstructMtxB(n1, flag_eps, vec_k)
global length_x length_y length_z
                            
      zone                = 1i; 
       
      cubic_length        = 1.0 ;       %2.0d0
      ratio_radius_length = 0.345;
      ratio_cyld_length   = 0.11;
      length_x            = cubic_length;
      length_y            = cubic_length;
      length_z            = cubic_length;
      rd_shpere           = ratio_radius_length * cubic_length;
      rd_cylinder         = ratio_cyld_length * cubic_length; 
      
      if ( flag_eps == 1 ) 
         ratio_epsilon       = 1.3d1;   % silicon-air
      else
         ratio_epsilon       = 1.14d1;  % GaAs-air
      end

      varepsilon_out      =  1.0;   % 0.3d0
      varepsilon_in       = ratio_epsilon * varepsilon_out;
                      
      ii    = 1;
      delta = zeros(3,1);
      if ( ii == 1 ) 
         delta(1)         = length_x / n1;
         delta(2)         = length_y / n1;
         delta(3)         = length_z / n1;
         n2               = n1;
         n3               = n1;
      else
         delta(1)         = length_x / n1;
         tmp              = length_y / delta(1);
         n2               = ( tmp - mod(tmp,1) ) + 1;
         delta(2)         = length_y / n2;
         tmp              = length_z / delta(1);
         n3               = ( tmp - mod(tmp,1) ) + 1;
         delta(3)         = length_z / n3;
      end 
      probdim             = 3 * n1 * n2 * n3;
    
      fprintf('n1 = %3.0f, n2 = %3.0f, n3 = %3.0f, dim = %7.0f \n', n1, n2, n3, 3*n1*n2*n3); 
                 
      diagB          = zeros(probdim,1); 
      inv_diag_D     = zeros(probdim,1);  
   
      idx = 0;

      for kk = 1: n3
         z_k = (kk-1) * delta(3);
         for jj = 1: n2
            y_j = (jj-1) * delta(2);
            for ii = 1: n1
               idx = idx + 1;
               x_i = ( ii - 0.5d0 ) * delta(1);

               tmp             = vec_k(1) * x_i + vec_k(2) * y_j + vec_k(3) * z_k;
               inv_diag_D(idx) = exp(-2.0d0 * pi * tmp * zone);
 
               in_out = locate_new(x_i, y_j, z_k, cubic_length, rd_shpere, rd_cylinder);
               if ( in_out == 'I' )  
                  diagB(idx)          = varepsilon_in;
               else
                  diagB(idx)          = varepsilon_out;
               end 
            end 
         end 
      end 
                  
      for kk = 1: n3
         z_k = (kk-1) * delta(3);
         for jj = 1: n2
            y_j = ( jj - 0.5d0 ) * delta(2);
            for ii = 1: n1
               idx = idx + 1;
               x_i = (ii-1) * delta(1);

               tmp             = vec_k(1) * x_i + vec_k(2) * y_j + vec_k(3) * z_k;
               inv_diag_D(idx) = exp(-2.0d0 * pi * tmp * zone);

               in_out = locate_new(x_i, y_j, z_k, cubic_length, rd_shpere, rd_cylinder);
               if ( in_out == 'I' ) 
                  diagB(idx)          = varepsilon_in;
               else
                  diagB(idx)          = varepsilon_out;
               end 
            end 
         end 
      end 
             
      for kk = 1: n3
         z_k = ( kk - 0.5d0 ) * delta(3);
         for jj = 1: n2
            y_j = (jj-1) * delta(2);
            for ii = 1: n1
               idx = idx + 1;
               x_i = (ii-1) * delta(1);

               tmp = vec_k(1) * x_i + vec_k(2) * y_j + vec_k(3) * z_k;
               inv_diag_D(idx) = exp(-2.0d0 * pi * tmp * zone);

               in_out = locate_new(x_i, y_j, z_k, cubic_length, rd_shpere, rd_cylinder);
               if ( in_out == 'I' )
                  diagB(idx)          = varepsilon_in;
               else
                  diagB(idx)          = varepsilon_out;
               end
            end
         end
      end
                               
end  %ConstructMtxB

%%%%%%%%%%%%%%%

function [ in_out ] = locate_new(x,y,z,cubic_length,rd_sphere,rd_cylinder)

      if ( x > 0.5d0*cubic_length )  
           tx = x - cubic_length;
      else
           tx = x;
      end
      if ( y > 0.5d0*cubic_length )  
           ty = y - cubic_length;
      else
           ty = y;
      end
      if ( z > 0.5d0*cubic_length )  
           tz = z - cubic_length;
      else
           tz = z;
      end
% --------------------------------------------------------------------

      if ( sqrt( tx^2 + ty^2 ) <= rd_cylinder )  
         in_out = 'I';
      elseif ( sqrt( ty^2 + tz^2 ) <= rd_cylinder )  
         in_out = 'I';
      elseif ( sqrt( tz^2 + tx^2 ) <= rd_cylinder )  
         in_out = 'I';
      elseif ( sqrt( tx^2 + ty^2 + tz^2 ) <= rd_sphere ) 
         in_out = 'I';
      else
         in_out = 'O';
      end
 
end %locate_new