全波傅氏算法matlab代码

以下是一种可能的全波傅氏算法的 MATLAB 代码实现：

function [E, H] = fdfd_fullwave(Nx, Ny, dx, dy, eps_r, mu_r, omega, src_pos, src_val)

% Nx, Ny: 网格的尺寸
% dx, dy: 网格的间距
% eps_r, mu_r: 网格中介质的相对电容率和相对磁导率
% omega: 光源的角频率
% src_pos: 光源在网格中的位置
% src_val: 光源的值

% 计算磁场
Hx = zeros(Nx+1, Ny);
Hy = zeros(Nx, Ny+1);

for j = 1:Ny
    for i = 1:Nx+1
        if i == src_pos(1) && j == src_pos(2)
            Hx(i,j) = src_val(1);
        elseif i == src_pos(1)+1 && j == src_pos(2)
            Hx(i,j) = src_val(2);
        else
            Hx(i,j) = 0;
        end
    end
end

for j = 1:Ny+1
    for i = 1:Nx
        if i == src_pos(1) && j == src_pos(2)
            Hy(i,j) = src_val(3);
        elseif i == src_pos(1) && j == src_pos(2)+1
            Hy(i,j) = src_val(4);
        else
            Hy(i,j) = 0;
        end
    end
end

% 计算电场
Ex = zeros(Nx, Ny+1);
Ey = zeros(Nx+1, Ny);

for j = 1:Ny+1
    for i = 1:Nx
        if j == src_pos(2) && i == src_pos(1)
            Ey(i,j) = src_val(5);
        elseif j == src_pos(2)+1 && i == src_pos(1)
            Ey(i,j) = src_val(6);
        else
            Ey(i,j) = 0;
        end
    end
end

for j = 1:Ny
    for i = 1:Nx+1
        if j == src_pos(2) && i == src_pos(1)
            Ex(i,j) = src_val(7);
        elseif j == src_pos(2) && i == src_pos(1)+1
            Ex(i,j) = src_val(8);
        else
            Ex(i,j) = 0;
        end
    end
end

% 傅里叶变换
FEx = fft2(Ex);
FEy = fft2(Ey);
FHx = fft2(Hx);
FHy = fft2(Hy);

% 计算傅里叶系数
kx_vec = (-Nx/2:Nx/2-1) * (2*pi/(Nx*dx));
ky_vec = (-Ny/2:Ny/2-1) * (2*pi/(Ny*dy));

[Kx, Ky] = meshgrid(kx_vec, ky_vec);
K = sqrt(Kx.^2 + Ky.^2);

epsilon = eps_r * ones(Nx, Ny);
mu = mu_r * ones(Nx, Ny);

Z = sqrt(mu ./ epsilon);

Zx = Z .* (Kx ./ K);
Zy = Z .* (Ky ./ K);

% 计算傅里叶系数的倒数
Zx_inv = 1 ./ Zx;
Zy_inv = 1 ./ Zy;

% 计算电场和磁场的傅里叶系数
FEx_new = (Zx_inv .* FHy - Zy_inv .* FEy) ./ (Zx_inv .* Zy_inv - K.^2);
FEy_new = -(Zy_inv .* FHx - Zx_inv .* FEx) ./ (Zx_inv .* Zy_inv - K.^2);
FHx_new = -(Zy_inv .* FEx - Zx_inv .* FEy) ./ (Zx_inv .* Zy_inv - K.^2);
FHy_new = (Zx_inv .* FEy - Zy_inv .* FHy) ./ (Zx_inv .* Zy_inv - K.^2);

% 计算电场和磁场
E = real(ifft2(FEx_new));
E = E(1:Nx, 1:Ny);
E = E / max(abs(E(:)));

E = E .* exp(-1i * omega * dx * src_pos(1));

H = real(ifft2(FHx_new));
H = H(1:Nx+1, 1:Ny);
H = H / max(abs(H(:)));

这个代码的输入是网格的大小 `Nx` 和 `Ny`，网格间距 `dx` 和 `dy`，介质的相对电容率和相对磁导率 `eps_r` 和 `mu_r`，光源的角频率 `omega`，光源在网格中的位置 `src_pos` 和光源的值 `src_val`。其中，`src_val` 是一个长度为 8 的向量，分别表示 x 和 y 方向的电场和 z 方向的磁场在两个网格点上的值。函数的输出是电场 `E` 和磁场 `H`。

这个代码实现的全波傅氏算法可以处理一般的线性均匀各向同性介质中的 Maxwell 方程，但是对于非线性或非均匀介质，或者介质表面有物体时，可能需要使用其他方法。