% 编码计算的黄金相对介电常数
function w = cef(z,N);
%可视化参数（用户输入）
lmin = 1e-7; %可视化中的最小波长
lmax = 1e-6; %可视化中的最大波长
dl = 5e-9; %可视化波长增量
N = 1000; %计算复杂误差函数时的迭代次数
%固有物理常数
c = 2.997e8;
h = 6.626e-34;
hbar = h/(2*pi);
q = 1.602e-19;
i = complex(0,1);
for lambda = lmin:dl:lmax;
omega = 2*pi*c./lambda;
%金属的数值常数
wp = 9.03;
f0 = 0.77;
fj = [0.054 0.050 0.312 0.719 1.648];
GAMMA0 = 0.050;
GAMMAj = [0.074 0.035 0.083 0.125 0.179];
wj = [0.218 2.885 4.069 6.137 27.97];
sigmaj = [0.742 0.349 0.830 1.246 1.795];
%常数和计算的转换
wp = q*wp/hbar;
GAMMA0 = (q/hbar)*GAMMA0;
GAMMAj = (q/hbar).*GAMMAj;
wj = (q/hbar).*wj;
sigmaj = (q/hbar).*sigmaj;
OMEGA_P = sqrt(f0).*wp;
%相对介电常数的第一部分
epsf = 1 - OMEGA_P^2./(omega.*(omega + i*GAMMA0));
%相对介电常数的第二（总和）分量
aj = (omega^2 + i*omega.*GAMMAj).^(0.5);
zplus = (aj+wj)./(sqrt(2).*sigmaj);
zminus = (aj-wj)./(sqrt(2).*sigmaj);
epsb = (i*sqrt(pi).*fj.*wp^2./(2*sqrt(2).*aj.*sigmaj)).*(cef(zplus,N)+cef(zminus,N));
%总介电常数计算
eps = epsf + sum(epsb);
%光波长上的可视化
figure(1)
loglog(lambda,abs(real(eps)),'o',lambda,abs(imag(eps)),'+')
hold on
axis([lmin lmax 0.1 2000])
xlabel('波长[m]')
ylabel('|\epsilon_{1r}| & \epsilon_{2r}')
title('黄金的相对介电常数\epsilon_{r} = \epsilon_{1r} + j\epsilon_{2r}')
legend
end
hold off
评论1