Implementation of a TFSF source in four sides of FDTD
6 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commenti
Risposte (0)
Vedere anche
Categorie
Scopri di più su Blue in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!