There are some bugs in the code which I can't deal with

[Samson-Ming]

Hi lydia1895 and all,

I have studied lydia1895's PMM code for a long time, though still I'm not quite understand Edee's paper on PMM in matched coordinate. I tried to look into the code step by step assisted with Edee's theory (Matched coordinates in the framework of polynomial modal methods for complex metasurface modeling, Vol. 35, No. 4 / April 2018 / Journal of the Optical Society of America A).

I made some trival modifications to the orignal code to accelerate the simulation, and also correct some bugs. But still I find some bugs that I can't deal with.

The code works quite well with stratified multi-layers film structures (still has a column in the unit cell, but its permittivity is the same with environment, i.e. epsilon(1)=epsilon(2), it is in fact a slab rather than grating), the output results are the same with those from the analytical transfer matrix method, no matter for lossless material or lossy material.

But when it comes to the patterned grating sturctue (the permitivity of the column in the unit cell is different from its environment, i.e. epsilon(1)≠epsilon(2)), there are some errors. For lossless grating, the results are inaccurate, only close to those from FMM or FEM methods, but never converge and never energy conservative. For lossy grating, the situation is much worse, for small N_b (numbers of polynominals used) and N_p (order of the polynominals, I modified the code to permit it to be different from 5 defined in the original code), the results is inaccurate, only close to those from FMM or FEM methods, but for large N_b and N_p, the results diverges (become very weird, either approximate to 0 or to infinity) due to the too large conditon number of variable "pminus_p1" and "Z" in the S method code (new_recursion.m) in matrix inversion.

I checked the code carefully, the apparent problem is that when we compute the eigenvalues in the PMM_gamma.m，some eigenvalues (especially for lossy grating, and the large the N_b, the more the "bad" eigenvalues) satisfy the Re(gamma)+Im(gamma)>0 but Im(gamma)<<0 (or Re(gamma)+Im(gamma)<0 but Im(gamma)>>0), thus makes the corresponding pplusv(m) = exp(1i*gammaplus(m)h)>>0 and diverges (say, 10^28) (or pminusv(m) = exp(-1igammaminus(m)*h)>>0 and diverges). The too large elements together with most orther elements (<=1) make the pplus and pminus matrix close to singular.

At first I guess there might be some problem in the eigvalue computation process, but because it works correctly for multi-layer films, I think the problem might not be here. So I looked into the part considering the material (PMM_epsilon_ellipse_matrices.m), and I bielve the problem can result from several lines in PMM_epsilon_ellipse_matrices.m below: . I tried for a while, but still can't get it right.

After a long time's struggle, I finally dicided to post my confusions and understandings here, I would hope lydia or anyone else can help me understand the principle and make it work, thank you so much!

Samson

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
.......
k1 = 2;
k2 = 2;

if (eps_in~=eps_out)

for m1=Nx(k1)+1:Nx(k1)+nx(k1)
for m2=Ny(k2)+1:Ny(k2)+ny(k2)
row = m2 + (m1-1)N_total_y3;
for i1=1:N_total_x3
for i2=1:N_total_y3
%%%%%%%%%%??????????????????????
%for i1=Nx(k1)+1:Nx(k1)+nx(k1)
%for i2=Ny(k2)+1:Ny(k2)+ny(k2)
col = i2 + (i1-1)N_total_y3;
eps_int_Ex_g_down22(row,col)=eps_inint_Ex_g_down22_full(row,col);
eps_int_Ey_g_down12(row,col)=eps_inint_Ey_g_down12_full(row,col);
eps_int_Ex_g_down21(row,col)=eps_inint_Ex_g_down21_full(row,col);
eps_int_Ey_g_down11(row,col)=eps_inint_Ey_g_down11_full(row,col);

    mu_int_Ex_g_down22(row,col) =mu_in*int_Ex_g_down22_full(row,col);
    mu_int_Ey_g_down12(row,col) =mu_in*int_Ey_g_down12_full(row,col);
    mu_int_Ex_g_down21(row,col) =mu_in*int_Ex_g_down21_full(row,col);
    mu_int_Ey_g_down11(row,col) =mu_in*int_Ey_g_down11_full(row,col);
    
    eps_int_Ez_g_sqrt(row,col)=eps_in*int_Ez_sqrt_g_full(row,col);
    mu_int_Ez_g_sqrt(row,col)=mu_in*int_Ez_sqrt_g_full(row,col);

end
end

end
end

end
.............
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Samson-Ming]

Also I read some PMM for a 1D grating with R method rather than S method, which works quite good, but I don't know if it matters. Neverthernless, whethe R method or S method, you need first get the eigenvalues, if several eigenvalues are abnormal, it must impact the R/S method and the final results.

[Samson-Ming]

Or there might be some problem in the definition of int_Ex_g_down22_full, int_Ex_g_down12_full,nt_Ey_g_down21_full,int_Ey_g_down11_full, and int_Ez_sqrt_g_full in PMM_metric_integral_polyfit_matrices.m file.