Implementation of a TFSF source in four sides of FDTD
6 views (last 30 days)
Show older comments
Hi all,
I am trying to correct my code by the implementation of the TFSF source in four sides of the FDTD problem, however, I have an issue in auxiliary griding of the code. I attached my code here, and would highly appreciate any help.
clear all; clc; close all;
% Grid Dimension in x (xdim) and y (ydim) directions
xN=400;
yN=xN;
%Total no of time steps
tN=600;
%Position of the source (center of the domain)
x_s=100;
y_s=100;
%Courant stability factor
S=1/(2^0.5);
% S=1/sqrt(3);
% Parameters of free space (permittivity and permeability and speed of
% light) are all not 1 and are given real values
epsilon0=(1/(36*pi))*1e-9;
mu0=4*pi*1e-7;
c=3e+8;
frequency=5e+13;
lambda=c/frequency;
% Spatial grid step length (spatial grid step= 1 micron and can be changed)
delx=lambda/10;
% Temporal grid step obtained using Courant condition
delt=S*delx/c;
% Initialization of field matrices
Ez=zeros(xN,yN);
Hy=zeros(xN,yN);
Hx=zeros(xN,yN);
Ezf=zeros(xN,yN);
Hyf=zeros(xN,yN);
Hxf=zeros(xN,yN);
Ezinc=zeros(xN,yN);
Hyinc=zeros(xN,yN);
Hxinc=zeros(xN,yN);
% Initialization of permittivity and permeability matrices
epsilon=epsilon0*ones(xN,yN);
mu=mu0*ones(xN,yN);
% Initializing electric and magnetic conductivity matrices
sigma=4e-4*ones(xN,yN);
sigma_star=4e-4*ones(xN,yN);
%Choice of nature of source
gaussian=0;
sine=0;
% The user can give a frequency of his choice for sinusoidal (if sine=1 above) waves in Hz
impulse=0;
%Choose any one as 1 and rest as 0. Default (when all are 0): Unit time step
%Multiplication factor matrices for H matrix update to avoid being calculated many times
%in the time update loop so as to increase computation speed
A=((mu-0.5*delt*sigma_star)./(mu+0.5*delt*sigma_star));
B=(delt/delx)./(mu+0.5*delt*sigma_star);
%Multiplication factor matrices for E matrix update to avoid being calculated many times
%in the time update loop so as to increase computation speed
C=((epsilon-0.5*delt*sigma)./(epsilon+0.5*delt*sigma));
D=(delt/delx)./(epsilon+0.5*delt*sigma);
F=(delt./epsilon)./(1+0.5*sigma*delt./epsilon)+ones(xN,yN); % Jsource coeffecient
Cb=(c*delt-delx)/(c*delt+delx);
one_way=1; %1 for One way BC, 0 for PEC
jsource=0; %1=Jsource ON, 0=Jsource OFF
%!-----------------Source parameters---------------------------------%
t=(1:tN)*delt;
Zs=50; % position index of source
delta=20*delt;
% Update loop begins
for n=1:1:tN
if jsource==1
% J source input
if n<=30
f(x_s,y_s)=sin(2*pi*frequency*(n)*delt)*exp(-(n-30)^2/30^2);
else
f(x_s,y_s)=sin(2*pi*frequency*(n)*delt);
end
end
%if source is impulse or unit-time step
if impulse==1
if n==1
Ez(x_s,y_s)=1;
else
Ez(x_s,y_s)=0;
end
end
%% Update fields
%Vector update instead of for-loop for Hy and Hx fields
Ez(Zs+1,Zs+2:yN-Zs-2)=(1-exp(-((delt*n).^2)/(delta^2)))*cos(frequency*delt*(n-3));
Hy(1:xN-1,1:yN-1)=A(1:xN-1,1:yN-1).*Hy(1:xN-1,1:yN-1)+B(1:xN-1,1:yN-1).*(Ez(1+1:xN-1+1,1:yN-1)-Ez(1:xN-1,1:yN-1));
Hx(1:xN-1,1:yN-1)=A(1:xN-1,1:yN-1).*Hx(1:xN-1,1:yN-1)-B(1:xN-1,1:yN-1).*(Ez(1:xN-1,1+1:yN-1+1)-Ez(1:xN-1,1:yN-1));
Ez(1+1:xN-1+1,1+1:yN-1+1)=C(1+1:xN-1+1,1+1:yN-1+1).*Ez(1+1:xN-1+1,1+1:yN-1+1)+D(1+1:xN-1+1,1+1:yN-1+1).*(Hy(1+1:xN-1+1,1+1:yN-1+1)-Hy(1:xN-1,1+1:yN-1+1)-Hx(1+1:xN-1+1,1+1:yN-1+1)+Hx(1+1:xN-1+1,1:yN-1));
% Ezinc(Zs,xN/2-50:xN/2+50)=cos(frequency*delt*(n-3));
Ezinc(Zs+1,Zs+2:yN-Zs-2)=(1-exp(-((delt*n).^2)/(delta^2)))*cos(frequency*delt*(n-3));
% Ezinc(Zs+1,Zs+2:yN-Zs-2)=(1-exp(-((delt*n).^2)/(delta^2)))*cos(frequency*delt*(n-3));
% Hyinc(Zs,:)=Am*cos(frequency.*(delt*n+s));
Hyinc(1:xN-1,1:yN-1)=A(1:xN-1,1:yN-1).*Hyinc(1:xN-1,1:yN-1)+B(1:xN-1,1:yN-1).*(Ezinc(1+1:xN-1+1,1:yN-1)-Ezinc(1:xN-1,1:yN-1));
Hxinc(1:xN-1,1:yN-1)=A(1:xN-1,1:yN-1).*Hxinc(1:xN-1,1:yN-1)-B(1:xN-1,1:yN-1).*(Ezinc(1:xN-1,1+1:yN-1+1)-Ezinc(1:xN-1,1:yN-1));
Ezinc(1+1:xN-1+1,1+1:yN-1+1)=C(1+1:xN-1+1,1+1:yN-1+1).*Ezinc(1+1:xN-1+1,1+1:yN-1+1)+D(1+1:xN-1+1,1+1:yN-1+1).*(Hyinc(1+1:xN-1+1,1+1:yN-1+1)-Hyinc(1:xN-1,1+1:yN-1+1)-Hxinc(1+1:xN-1+1,1+1:yN-1+1)+Hxinc(1+1:xN-1+1,1:yN-1));
Hy(Zs-1,Zs:xN-Zs)=Hy(Zs-1,Zs:xN-Zs)-B(Zs,Zs:xN-Zs).*Ezinc(Zs,Zs:xN-Zs);
Hy(xN-Zs+1,Zs:xN-Zs)=Hy(xN-Zs+1,Zs:xN-Zs)+B(xN-Zs,Zs:xN-Zs).*Ezinc(xN-Zs,Zs:xN-Zs);
Hx(Zs:xN-Zs,Zs-1)=Hx(Zs:xN-Zs,Zs-1)+B(Zs:xN-Zs,Zs).*Ezinc(Zs:xN-Zs,Zs);
Hx(Zs:xN-Zs,xN-Zs+1)=Hx(Zs:xN-Zs,xN-Zs+1)-B(Zs:xN-Zs,xN-Zs).*Ezinc(Zs:xN-Zs,xN-Zs);
Ez(Zs,Zs:xN-Zs)=Ez(Zs,Zs:xN-Zs)-D(Zs,Zs:xN-Zs).*Hyinc(Zs-1,Zs:xN-Zs);
Ez(xN-Zs,Zs:xN-Zs)=Ez(xN-Zs,Zs:xN-Zs)+D(xN-Zs-1,Zs:xN-Zs).*Hyinc(xN-Zs-1,Zs:xN-Zs);
Ez(Zs:xN-Zs,Zs)=Ez(Zs:xN-Zs,Zs)+D(Zs:xN-Zs,Zs).*Hxinc(Zs:xN-Zs,Zs-1);
Ez(Zs:xN-Zs,xN-Zs)=Ez(Zs:xN-Zs,xN-Zs)-D(Zs:xN-Zs,xN-Zs).*Hxinc(Zs:xN-Zs,xN-Zs+1);
if jsource ==1
Ez(x_s,y_s)=Ez(x_s,y_s)-F(x_s,y_s)*f(x_s,y_s); %J_source update
end
%% temporary storage of fields for getting Ez(n+1) (required for One way BC)
Hyf=Hy;
Hxf=Hx;
Ezf=Ez;
Hyf(1:xN-1,1:yN-1)=A(1:xN-1,1:yN-1).*Hyf(1:xN-1,1:yN-1)+B(1:xN-1,1:yN-1).*(Ezf(1+1:xN-1+1,1:yN-1)-Ezf(1:xN-1,1:yN-1));
Hxf(1:xN-1,1:yN-1)=A(1:xN-1,1:yN-1).*Hxf(1:xN-1,1:yN-1)-B(1:xN-1,1:yN-1).*(Ezf(1:xN-1,1+1:yN-1+1)-Ezf(1:xN-1,1:yN-1));
Ezf(1+1:xN-1+1,1+1:yN-1+1)=C(1+1:xN-1+1,1+1:yN-1+1).*Ezf(1+1:xN-1+1,1+1:yN-1+1)+D(1+1:xN-1+1,1+1:yN-1+1).*(Hyf(1+1:xN-1+1,1+1:yN-1+1)-Hyf(1:xN-1,1+1:yN-1+1)-Hxf(1+1:xN-1+1,1+1:yN-1+1)+Hxf(1+1:xN-1+1,1:yN-1));
if jsource==1
Ezf(x_s,y_s)=Ezf(x_s,y_s)-F(x_s,y_s)*f(x_s,y_s);
end
%% Boundary condition
if one_way==1
Ez(1,:)=Ez(2,:)+Cb.*(Ezf(2,:)-Ez(1,:)); %at x=0
Ez(xN,:)=Ez(xN-1,:)+Cb.*(Ezf(xN-1,:)-Ez(xN,:)); %at x=l
Ez(:,1)=Ez(:,2)+Cb.*(Ezf(:,2)-Ez(:,1)); %at y=0
Ez(:,yN)=Ez(:,yN-1)+Cb.*(Ezf(:,yN-1)-Ez(:,yN)); %at y=l
else
% Perfect Electric Conductor boundary condition
Ez(1,:)=0;
Ez(xN,:)=0;
Ez(:,1)=0;
Ez(:,yN)=0;
end
%% Hard Source type
if impulse==0 && jsource==0
% If unit-time step
% if gaussian==0 && sine==0
% Ez(x_s,y_s)=0;
% end
%if sine
if sine==1
tstart=1;
N_lambda=c/(frequency*delx);
Ez(x_s,y_s)=sin(((2*pi*(c/(delx*N_lambda))*(n-tstart)*delt)));
% Ez(:,y_s)=sin(((2*pi*(c/(delx*N_lambda))*(n-tstart)*delt)));
end
%if gaussian
if gaussian==1
if n<=42
Ez(x_s,y_s)=(10-15*cos(n*pi/20)+6*cos(2*n*pi/20)-cos(3*n*pi/20))/32;
else
Ez(x_s,y_s)=0;
end
end
end
%% Movie plot of Ez
if jsource==1
imagesc(delx*(1:1:xN)*1e+6,(1e+6*delx*(1:1:yN))',Ez',[-0.05 0.05]); colormap jet;
else
imagesc(delx*(1:1:xN)*1e+6,(1e+6*delx*(1:1:yN))',Ez',[-1 1]); colormap jet;colorbar
end
title(['\fontsize{20}E_{z} at t = ',num2str(round(n*delt*1e+15)),' fs']);
xlabel('x (in um)','FontSize',20);
ylabel('y (in um)','FontSize',20);
set(gca,'FontSize',20);
getframe;
end
0 Comments
Answers (0)
See Also
Categories
Find more on Blue in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!