电磁偶极辐射在平面分层介质中，计算由任意层数的平面多层介质中的电偶极或磁偶极产生的近场电场（Matlab代码实现）

 💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

 ⛳️赠与读者

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

 ⛳️赠与读者

👨‍💻做科研，涉及到一个深在的思想系统，需要科研者逻辑缜密，踏实认真，但是不能只是努力，很多时候借力比努力更重要，然后还要有仰望星空的创新点和启发点。当哲学课上老师问你什么是科学，什么是电的时候，不要觉得这些问题搞笑。哲学是科学之母，哲学就是追究终极问题，寻找那些不言自明只有小孩子会问的但是你却回答不出来的问题。建议读者按目录次序逐一浏览，免得骤然跌入幽暗的迷宫找不到来时的路，它不足为你揭示全部问题的答案，但若能让人胸中升起一朵朵疑云，也未尝不会酿成晚霞斑斓的别一番景致，万一它居然给你带来了一场精神世界的苦雨，那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。

     或许，雨过云收，神驰的天地更清朗.......🔎🔎🔎

💥1 概述

电磁偶极辐射在平面分层介质中，计算由任意层数的平面多层介质中的电偶极或磁偶极产生的近场电场。
该代码计算了在平面分层介质环境中由电偶极或磁偶极（J或M）辐射产生的近场。该代码解决了Sommerfeld积分，并使用传输线模型处理了设置。该代码对偶极位置或方向没有限制。此外，可以使用任意数量的层（请修改SampleConfiguration.m文件以创建您的配置）。介质参数被假定为各向同性，但可以接受复杂值（损耗和等离子体），以及LH介质参数。除此之外，用户可以在界面上添加表面导电片来模拟类似各向同性石墨烯或任何类似的二维金属等场景。此外，用户可以在顶部和/或底部层添加PEC/PMC终端以考虑类似于带有GND平面的基板或平行板波导等场景。

注意：计算Sommerfeld积分很耗时，尽管代码有加速功能，计算大量点数（>100x100）可能需要几分钟或几小时。

在探究平面分层介质结构中的电磁现象时，准确模拟由嵌入其中的电偶极或磁偶极源产生的近场电场是至关重要的研究内容。本研究提供的先进计算代码，专注于在任意配置的多层介质体系中，精确计算电偶极（J）或磁偶极（M）辐射源激发的近场电场分布。这项工作的核心在于高效解决Sommerfeld积分问题，一个在波动理论中极具挑战性的任务，同时巧妙地借助传输线模型来精确描绘和处理这种复杂的介质分层环境。

此代码的灵活性与通用性尤为突出：无论是偶极源的位置、方向，还是介质结构的复杂性，都几乎不受限制。用户仅需调整SampleConfiguration.m文件，即可轻松设定从简单到极端复杂的任意层数介质配置，每一层的介质参数均可独立定义，充分体现了其高度自定义的能力。值得注意的是，虽然默认假设介质参数为各向同性，但代码同样支持包含损耗因子和等离子体效应的复数介质参数，甚至可以纳入左手介质（LH介质）特性，进一步拓宽了应用范围。

在界面处理方面，代码还赋予了用户前所未有的控制权。用户能够自由地在任一介质层界面引入表面导电片层，借此模拟现代纳米技术中常见的各向同性石墨烯层或其他二维材料的电磁响应。此外，为了满足更多研究需求，代码还支持在体系的顶部和/或底部添加理想导电板（PEC）或理想磁导板（PMC），从而能够有效地模拟诸如具有接地平面的印刷电路板（PCB）基板、平行板波导等工程技术中的典型场景。

值得注意的是，尽管代码内置了优化措施以加快计算速度，但由于Sommerfeld积分固有的计算密集性，当需要处理大规模网格点（例如超过100x100的点阵）时，完成计算仍可能耗时较长，从几分钟到数小时不等，这取决于具体计算环境的性能和配置复杂度。因此，用户在进行大规模模拟前，建议合理规划计算资源，以优化整体研究效率和体验。

总而言之，本研究提供的代码代表了电磁场仿真领域的一大进步，它不仅极大地推动了对平面分层介质中电磁偶极辐射行为的深入理解，也为涉及复杂介质结构的工程设计和纳米光子学研究提供了强大的计算工具，其灵活度、准确性和实用性均为相关科研工作者带来了极大的便利。

📚2 运行结果

主函数代码：

close all; clear; clc;
addpath(strcat(pwd,'\Engine'));
%%
Ns          =   500;        % Number of data points
[~,~,~,~,nm,~,~,~,~]=Units();
%% Definitions
Type        =   'J';        % Dipole type 'J' or 'M'
theta0      =   0;            % Tilt theta angle of dipole 
phi0        =   0;          % Tilt phi angle of dipole
J           =   1;          % Dipole strength
x_          =   0*nm;       % Dipole x' position
y_          =   0*nm;       % Dipole y' position
z_          =   20*nm;      % Dipole z' position
%% Results definitions
x_min        =   -2000*nm;   % Results minimum x value
x_max        =   +2000*nm;    % Results maximum x value
y           =   0*nm;       % Results y-plane
z_min         =   -2000*nm;   % Results minimum z value
z_max         =   +2000*nm;    % Results maximum z value
%% Corrections
[x,z]       =   meshgrid(linspace(x_min,x_max,Ns),linspace(z_min,z_max,Ns));
[phi,rho]   =   cart2pol(x-x_,y-y_);
%% Dipole Components
[Jx,Jy,Jz]=DipoleJ(theta0,phi0,J);
%% Compute Fields
Ex          =   zeros(Ns,Ns);
Ey          =   zeros(Ns,Ns);
Ez          =   zeros(Ns,Ns);
tic;
count       =   0;
for m=1:Ns
    for n=1:Ns
        if Type=='J'
            Ex(m,n)     =   GEJxx(rho(m,n),phi(m,n),z(m,n),z_)*Jx+GEJxy(rho(m,n),phi(m,n),z(m,n),z_)*Jy+GEJxz(rho(m,n),phi(m,n),z(m,n),z_)*Jz;
            Ey(m,n)     =   GEJyx(rho(m,n),phi(m,n),z(m,n),z_)*Jx+GEJyy(rho(m,n),phi(m,n),z(m,n),z_)*Jy+GEJyz(rho(m,n),phi(m,n),z(m,n),z_)*Jz;
            Ez(m,n)     =   GEJzx(rho(m,n),phi(m,n),z(m,n),z_)*Jx+GEJzy(rho(m,n),phi(m,n),z(m,n),z_)*Jy+GEJzz(rho(m,n),phi(m,n),z(m,n),z_)*Jz;
        end
        if Type=='M'
            Ex(m,n)     =   GEMxx(rho(m,n),phi(m,n),z(m,n),z_)*Jx+GEMxy(rho(m,n),phi(m,n),z(m,n),z_)*Jy+GEMxz(rho(m,n),phi(m,n),z(m,n),z_)*Jz;
            Ey(m,n)     =   GEMyx(rho(m,n),phi(m,n),z(m,n),z_)*Jx+GEMyy(rho(m,n),phi(m,n),z(m,n),z_)*Jy+GEMyz(rho(m,n),phi(m,n),z(m,n),z_)*Jz;
            Ez(m,n)     =   GEMzx(rho(m,n),phi(m,n),z(m,n),z_)*Jx+GEMzy(rho(m,n),phi(m,n),z(m,n),z_)*Jy+GEMzz(rho(m,n),phi(m,n),z(m,n),z_)*Jz;
        end
        count       =   count+1;
        clc;
        fprintf('Step:\t%0.0f/100\n',count*100/Ns^2);
    end
end
toc; clc;
fprintf('Time Elapsed %0.2f minutes\n',toc/60);
%% E Field
E_abs        =   sqrt(abs(Ex).^2+abs(Ey).^2+abs(Ez).^2);
E_real        =   sqrt(real(Ex).^2+real(Ey).^2+real(Ez).^2);
%% Plot E
figure()
pcolor(x/nm,z/nm,20*log10(E_abs))
shading interp
colormap jet
colorbar
hold on
PlotLayersH(x_min,x_max,nm);
hold off
xlabel('$x$ [nm]','Interpret','Latex')
ylabel('$z$ [nm]','Interpret','Latex')
title('$20\log_{10}|E|$','Interpret','Latex')
set(gca,'TickLabel','Latex')
set(colorbar,'TickLabelInterpreter','Latex')
axis equal
axis([x_min x_max z_min z_max]/nm)
%% Plot Ez
figure()
pcolor(x/nm,z/nm,20*log10(E_real))
shading interp
colormap jet
colorbar
hold on
PlotLayersH(x_min,x_max,nm);
hold off
xlabel('$x$ [nm]','Interpret','Latex')
ylabel('$z$ [nm]','Interpret','Latex')
title('$20\log_{10}|\textrm{Re}[E]|$','Interpret','Latex')
set(gca,'TickLabel','Latex')
set(colorbar,'TickLabelInterpreter','Latex')
axis equal
axis([x_min x_max z_min z_max]/nm)
caxis([200 340])
%%
function[Jx,Jy,Jz]=DipoleJ(theta0,phi0,J)
theta0      =   deg2rad(theta0);
phi0        =   deg2rad(phi0);
Jx          =   J*sin(theta0)*cos(phi0);
Jy          =   J*sin(theta0)*sin(phi0);
Jz          =   J*cos(theta0);
end
%%
 

🎉3 参考文献

文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。

[1]汪智.电磁理论中Sommerfeld积分的新解法[D].电子科技大学[2024-07-03].

[2]徐建华,朱德怀.多层介质中偶极子场的系数推关系[J].石油地球物理勘探, 1994.

[3]周金中.点偶极子在多层同心介质中的热力学函数[J].自然杂志, 1982(08):76-77.DOI:CNKI:SUN:ZRZZ.0.1982-08-022.

🌈4 Matlab代码实现

资料获取，更多粉丝福利，MATLAB|Simulink|Python资源获取