Continue to Site

Welcome to EDAboard.com

Welcome to our site! EDAboard.com is an international Electronics Discussion Forum focused on EDA software, circuits, schematics, books, theory, papers, asic, pld, 8051, DSP, Network, RF, Analog Design, PCB, Service Manuals... and a whole lot more! To participate you need to register. Registration is free. Click here to register now.

FDTD C source code from the book

Status
Not open for further replies.

tgfet2002

Newbie level 2
Joined
Jun 6, 2006
Messages
2
Helped
1
Reputation
2
Reaction score
1
Trophy points
1,283
Activity points
1,288
fdtd source code

1D FDTD
2D TDTD
3D FDTD
 

fdtd c

you can refer to link:

this ebook is electromagnetic simulation by FDTD (by sullivan)

and link:

this ebook is electromagnetic by FDTD( by kunz)

and code matlab from link:
https://www.mathworks.com/matlabcentral/fileexchange/loadCategory.do
you must search about FDTD then start download code matlab.
this is example 1D about method FDTD:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Scott Hudson, WSU Tri-Cities
%1D electromagnetic finite-difference time-domain (FDTD) program.
%Assumes Ey and Hz field components propagating in the x direction.
%Fields, permittivity, permeability, and conductivity
%are functions of x. Try changing the value of "profile".
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;
clear all;

L = 5; %domain length in meters
N = 505; %# spatial samples in domain
Niter = 800; %# of iterations to perform
fs = 300e6; %source frequency in Hz
ds = L/N; %spatial step in meters
dt = ds/300e6; %"magic time step"
eps0 = 8.854e-12; %permittivity of free space
mu0 = pi*4e-7; %permeability of free space
x = linspace(0,L,N); %x coordinate of spatial samples
showWKB=0; %if =1 then show WKB appromination at end

%scale factors for E and H
ae = ones(N,1)*dt/(ds*eps0);
am = ones(N,1)*dt/(ds*mu0);
as = ones(N,1);
epsr = ones(N,1);
mur= ones(N,1);
sigma = zeros(N,1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Here we specify the epsilon, sigma, and mu profiles. I've
%predefined some interesting examples. Try profile = 1,2,3,4,5,6 in sequence.
%You can define epsr(i), mur(i) (relative permittivity and permeability)
%and sigma(i) (conductivity) to be anything you want.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
profile = 1;
for i=1:N
epsr(i) = 1;
mur(i) = 1;
w1 = 0.5;
w2 = 1.5;
if (profile==1) %dielectric window
if (abs(x(i)-L/2)<0.5) epsr(i)=4; end
end
if (profile==2)%dielectric window with smooth transition
if (abs(x(i)-L/2)<1.5) epsr(i)=1+3*(1+cos(pi*(abs(x(i)-L/2)-w1)/(w2-w1)))/2; end
if (abs(x(i)-L/2)<0.5) epsr(i)=4; end
end
if (profile==3) %dielectric discontinuity
if (x(i)>L/2) epsr(i) = 9; end
end
if (profile==4) %dielectric disontinuity with 1/4-wave matching layer
if (x(i)>(L/2-0.1443)) epsr(i) = 3; end
if (x(i)>L/2) epsr(i) = 9; end
end
if (profile==5) %conducting half space
if (x(i)>L/2) sigma(i) = 0.005; end
end
if (profile==6) %sinusoidal dielectric
epsr(i) = 1+sin(2*pi*x(i)/L)^2;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


ae = ae./epsr;
am = am./mur;
ae = ae./(1+dt*(sigma./epsr)/(2*eps0));
as = (1-dt*(sigma./epsr)/(2*eps0))./(1+dt*(sigma./epsr)/(2*eps0));

%plot the permittivity, permeability, and conductivity profiles
figure(1)
subplot(3,1,1);
plot(x,epsr);
grid on;
axis([3*ds L min(epsr)*0.9 max(epsr)*1.1]);
title('relative permittivity');
subplot(3,1,2);
plot(x,mur);
grid on;
axis([3*ds L min(mur)*0.9 max(mur)*1.1]);
title('relative permeabiliity');
subplot(3,1,3);
plot(x,sigma);
grid on;
axis([3*ds L min(sigma)*0.9-0.001 max(sigma)*1.1+0.001]);
title('conductivity');

%initialize fields to zero
Hz = zeros(N,1);
Ey = zeros(N,1);
figure(2);
set(gcf,'doublebuffer','on'); %set double buffering on for smoother graphics
plot(Ey);
grid on;

for iter=1:Niter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The next 10 or so lines of code are where we actually integrate Maxwell's
%equations. All the rest of the program is basically bookkeeping and plotting.
%"smooth turn on" sinusoidal source
Ey(3) = Ey(3)+2*(1-exp(-((iter-1)/50)^2))*sin(2*pi*fs*dt*iter);
Hz(1) = Hz(2); %absorbing boundary conditions for left-propagating waves
for i=2:N-1 %update H field
Hz(i) = Hz(i)-am(i)*(Ey(i+1)-Ey(i));
end
Ey(N) = Ey(N-1); %absorbing boundary conditions for right propagating waves
for i=2:N-1 %update E field
Ey(i) = as(i)*Ey(i)-ae(i)*(Hz(i)-Hz(i-1));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
hold off
plot(x,Ey,'b');
axis([3*ds L -2 2]);
grid on;
title('E (blue) and 377*H (red)');
hold on
plot(x,377*Hz,'r');
xlabel('x (m)');
pause(0);
iter
end

%WKB prediction for the fields
if (showWKB==1)
phase = cumsum((epsr).^0.5)*ds;
beta0 = 2*pi*fs/(300e6);
theory = sin(2*pi*fs*(Niter+4)*dt-beta0*phase)./(epsr.^0.25);
input('press enter to show WKB theory');
plot(x,theory,'b.');
theory = sin(2*pi*fs*(Niter+4)*dt-beta0*phase).*(epsr.^0.25);
plot(x,theory,'r.');
title('E (blue), 377*H (red), WKB theory (points)');
end
 

Status
Not open for further replies.

Part and Inventory Search

Welcome to EDABoard.com

Sponsor

Back
Top