zozokim
Newbie
I've been studying the UPML portion of Computational electrodynamics by ALLEN TAFLOVE and writing matlab code, but the electromagnetic field is not updating.
Please help me
I can't remove things like emoticon :/), so I'm attaching it as a notepad.
Thank you.
I apologize for my poor English.
This is my MATLAB code.
Please help me
I can't remove things like emoticon :/), so I'm attaching it as a notepad.
Thank you.
I apologize for my poor English.
This is my MATLAB code.
Code:
--------------------------------------------------------------------------------------------------------
clear all; close all; clc;
%%
%Defining constants
e0 = 8.854*10^-12; % Permittivity of vacuum [F/m]
u0 = 4*pi*10^-7; % Permeability of vacuum [H/m]
c0 = 1/(e0*u0)^0.5; % Speed of light [m/s]
f0 = 10e6; % Frequency of Source [Hz]
w0 = 2*pi*f0; % Angular frequency [rad/s]
l0 = c0/f0; % lambda [m]
X = 100; % X size [m]
Y = 100; % Y size [m]
Z = 100; % Z size [m]
Mx = 100; % X cell
My = 100; % Y cell
Mz = 100; % Z cell
N = 2000; % Time step
dx = X/Mx; % Increment of X
dy = Y/My; % Increment of Y
dz = Z/Mz; % Increment of Z
dt = (dx^-2+dy^-2+dz^-2)^(-0.5)/c0; % Courant number (Magic time step)
x_axis = (1:Mx)*dx; % Physical length of X
y_axis = (1:My)*dy; % Physical length of X
sigx = zeros(Mx+1,My+1,Mz+1); % Conductivity x-axis
sigy = zeros(Mx+1,My+1,Mz+1); % Conductivity y-axis
sigz = zeros(Mx+1,My+1,Mz+1); % Conductivity z-axis (In 2D, sigz = 0)
e = e0*ones(Mx+1,My+1,Mz+1); % Permittivity matrix of PML region
u = u0*ones(Mx+1,My+1,Mz+1); % Permeability matrix of PML region
eta = sqrt(u0/e0); % Intrinsic impedance of PML region
%PML condition(polynomial grading)
p = 25; % PML grid
kx = 1; % Kappa x
ky = 1; % Kappa y
kz = 1; % Kappa z (In 2D, Kappa z = 1)
m = 4; % Grading order
R0 = exp(-16); % Reflection errors
sig_max = -((m+1)*log(R0))/(2*eta*p); % Sigma max
for n = 1:p
sigx(p+1-n,:,:) = (n/p)^m*sig_max;
sigy(:,p+1-n,:) = (n/p)^m*sig_max;
sigz(:,:,p+1-n) = (n/p)^m*sig_max;
sigx(end-p+n,:,:) = (n/p)^m*sig_max;
sigy(:,end-p+n,:) = (n/p)^m*sig_max;
sigz(:,:,end-p+n) = (n/p)^m*sig_max;
end
%Initializing dimensions
Dx = zeros(Mx,My+1,Mz+1);
Ex = zeros(Mx,My+1,Mz+1);
Dy = zeros(Mx+1,My,Mz+1);
Ey = zeros(Mx+1,My,Mz+1);
Dz = zeros(Mx+1,My+1,Mz);
Ez = zeros(Mx+1,My+1,Mz);
Bx = zeros(Mx+1,My,Mz);
Hx = zeros(Mx+1,My,Mz);
By = zeros(Mx,My+1,Mz);
Hy = zeros(Mx,My+1,Mz);
Bz = zeros(Mx,My,Mz+1);
Hz = zeros(Mx,My,Mz+1);
% monitor = zeros(Mx,My,N);
Wx = 2*e0*kx-sigx*dt; %Coefficient of update equations
Wy = 2*e0*ky-sigy*dt; %Coefficient of update equations
Wz = 2*e0*kz-sigz*dt; %Coefficient of update equations
Qx = 2*e0*kx+sigx*dt; %Coefficient of update equations
Qy = 2*e0*ky+sigy*dt; %Coefficient of update equations
Qz = 2*e0*kz+sigz*dt; %Coefficient of update equations
clims = [-70 -30];
% v = VideoWriter('2D FDTD PML1');
% open(v);
%%
%UPML Update equation by auxiliary differential equation
% Source
t = (1:N)*dt;
E_in = exp(-2 * pi^2 * (f0)^2 * (t - 1/(f0)).^2 );
figure(2);
plot(t, E_in);
i1 = 1:Mx;
j1 = 1:My;
k1 = 1:Mz;
i2 = 2:Mx;
j2 = 2:My;
k2 = 2:Mz;
% Define coefficients of update equations
T1 = (2*e0*dt);
Dx_C1 = Wy(i1,j2,k2)./Qy(i1,j2,k2);
Dx_C2 = T1./Qy(i1,j2,k2)/dy;
Dx_C3 = T1./Qy(i1,j2,k2)/dz;
Dy_C1 = Wz(i2,j1,k2)./Qz(i2,j1,k2);
Dy_C2 = T1./Qz(i2,j1,k2)/dz;
Dy_C3 = T1./Qz(i2,j1,k2)/dx;
Dz_C1 = Wx(i2,j2,k1)./Qx(i2,j2,k1);
Dz_C2 = T1./Qx(i2,j2,k1)/dx;
Dz_C3 = T1./Qx(i2,j2,k1)/dy;
Ex_C1 = Wz(i1,j2,k2)./Qz(i1,j2,k2);
Ex_C2 = 1./(Qz(i1,j2,k2).*e(i1,j2,k2)).*Qx(i1,j2,k2);
Ex_C3 = 1./(Qz(i1,j2,k2).*e(i1,j2,k2)).*Wx(i1,j2,k2);
Ey_C1 = Wx(i2,j1,k2)./Qx(i2,j1,k2);
Ey_C2 = 1./(Qx(i2,j1,k2).*e(i2,j1,k2)).*Qy(i2,j1,k2);
Ey_C3 = 1./(Qx(i2,j1,k2).*e(i2,j1,k2)).*Wy(i2,j1,k2);
Ez_C1 = Wy(i2,j2,k1)./Qy(i2,j2,k1);
Ez_C2 = 1./(Qy(i2,j2,k1).*e(i2,j2,k1)).*Qz(i2,j2,k1);
Ez_C3 = 1./(Qy(i2,j2,k1).*e(i2,j2,k1)).*Wz(i2,j2,k1);
Bx_C1 = Wy(:,j1,k1)./Qy(:,j1,k1);
Bx_C2 = T1./Qy(:,j1,k1)/dy;
Bx_C3 = T1./Qy(:,j1,k1)/dz;
By_C1 = Wz(i1,:,k1)./Qz(i1,:,k1);
By_C2 = T1./Qz(i1,:,k1)/dz;
By_C3 = T1./Qz(i1,:,k1)/dx;
Bz_C1 = Wx(i1,j1,:)./Qx(i1,j1,:);
Bz_C2 = T1./Qx(i1,j1,:)/dx;
Bz_C3 = T1./Qx(i1,j1,:)/dy;
Hx_C1 = Wz(:,j1,k1)./Qz(:,j1,k1);
Hx_C2 = 1./(Qz(:,j1,k1).*u(:,j1,k1)).*Qx(:,j1,k1);
Hx_C3 = 1./(Qz(:,j1,k1).*u(:,j1,k1)).*Wx(:,j1,k1);
Hy_C1 = Wx(i1,:,k1)./Qx(i1,:,k1);
Hy_C2 = 1./(Qx(i1,:,k1).*u(i1,:,k1)).*Qy(i1,:,k1);
Hy_C3 = 1./(Qx(i1,:,k1).*u(i1,:,k1)).*Wy(i1,:,k1);
Hz_C1 = Wy(i1,j1,:)./Qy(i1,j1,:);
Hz_C2 = 1./(Qy(i1,j1,:).*u(i1,j1,:)).*Qz(i1,j1,:);
Hz_C3 = 1./(Qy(i1,j1,:).*u(i1,j1,:)).*Wz(i1,j1,:);
figure(1);
for n = 1:N
t(n)
Dx_temp = Dx(:,j2,k2);
Dy_temp = Dy(i2,:,k2);
Dz_temp = Dz(i2,j2,:);
Ex_temp = Ex(:,j2,k2);
Ey_temp = Ey(i2,:,k2);
Ez_temp = Ez(i2,j2,:);
Bx_temp = Bx(:,j1,k1);
By_temp = By(i1,:,k1);
Bz_temp = Bz(i1,j1,:);
Hx_temp = Hx(:,j1,k1);
Hy_temp = Hy(i1,:,k1);
Hz_temp = Hz(i1,j1,:);
Dx(:,j2,k2) = Dx_C1.*Dx_temp + Dx_C2.*(Hz(:,j2,k2) - Hz(:,j2-1,k2)) - Dx_C3.*(Hy(:,j2,k2) - Hy(:,j2,k2-1));
Dy(i2,:,k2) = Dy_C1.*Dy_temp + Dy_C2.*(Hx(i2,:,k2) - Hx(i2,:,k2-1)) - Dy_C3.*(Hz(i2,:,k2) - Hz(i2-1,:,k2));
Dz(i2,j2,:) = Dz_C1.*Dz_temp + Dz_C2.*(Hy(i2,j2,:) - Hy(i2-1,j2,:)) - Dz_C3.*(Hx(i2,j2,:) - Hx(i2,j2-1,:));
Ex(:,j2,k2) = Ex_C1.*Ex_temp + Ex_C2.*Dx(:,j2,k2) - Ex_C3.*Dx_temp;
Ey(i2,:,k2) = Ey_C1.*Ey_temp + Ey_C2.*Dy(i2,:,k2) - Ey_C3.*Dy_temp;
Ez(i2,j2,:) = Ez_C1.*Ez_temp + Ez_C2.*Dz(i2,j2,:) - Ez_C3.*Dz_temp;
Bx(:,j1,k1) = Bx_C1.*Bx_temp - Bx_C2.*(Ez(:,j1+1,k1) - Ez(:,j1,k1)) + Bx_C3.*(Ey(:,j1,k1+1) - Ey(:,j1,k1));
By(i1,:,k1) = By_C1.*By_temp - By_C2.*(Ex(i1,:,k1+1) - Ex(i1,:,k1)) + By_C3.*(Ez(i1+1,:,k1) - Ez(i1,:,k1));
Bz(i1,j1,:) = Bz_C1.*Bz_temp - Bz_C2.*(Ey(i1+1,j1,:) - Ey(i1,j1,:)) + Bz_C3.*(Ex(i1,j1+1,:) - Ex(i1,j1,:));
Hx(:,j1,k1) = Hx_C1.*Hx_temp + Hx_C2.*Bx(:,j1,k1) - Hx_C3.*Bx_temp;
Hy(i1,:,k1) = Hy_C1.*Hy_temp + Hy_C2.*By(i1,:,k1) - Hy_C3.*By_temp;
Hz(i1,j1,:) = Hz_C1.*Hz_temp + Hz_C2.*Bz(i1,j1,:) - Hz_C3.*Bz_temp;
Ez(round(X/2),round(Y/2),round(Z/2)) = Ez(round(X/2),round(Y/2),round(Z/2)) - exp(-2 * pi^2 * (f0)^2 * (n*dt - 1/(f0)).^2 );
E_total = sqrt(abs(Ex(i1,j1,k1)).^2+abs(Ey(i1,j1,k1)).^2+abs(Ez(i1,j1,k1)).^2);
% subplot(1,2,1)
%pcolor(db(E_total(:,:,round(Mz/2))))
figure(1);
pcolor(y_axis, x_axis, (E_total(:,:,round(Z/2))))
shading flat
caxis([0 5]);
colorbar;
drawnow;
% frame = getframe(gcf);
% writeVideo(v,frame);
end
% close(v);
---------------------------------------------------------------------------------------------------------------------------
Last edited by a moderator: