FDTD C zdrojový kód z knihy

[tgfet2002]

1D FDTD
2D TDTD
3D FDTD

 

[nazanin]

Můžete odkazovat na odkaz:
http://www.edaboard.com/viewtopic.php?t=74057
to ebook je elektromagnetické simulace FDTD (podle Sullivan)

a odkaz:
http://www.edaboard.com/viewtopic.php?t=149266
to ebook je elektromagnetické podle FDTD (podle Kunz)

a kód MATLAB z odkazu:
http://www.mathworks.com/matlabcentral/fileexchange/loadCategory.do
Musíte hledat asi FDTD potom začít stahovat kódu MATLABu.
To je příklad 1D o metodu FDTD:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%
% Scott Hudson, WSU Tri-Cities
% 1D elektromagnetické konečný-rozdíl čas-domain (FDTD) program.
% Předpokládá Ey a Hz oblasti komponentů rozmnožovací v x směr.
% Pole, permitivita, permeabilita, vodivost a
% Jsou funkce x.Zkuste změnit hodnotu "profil".
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%

Zavřete všechny;
Vymazat vše;

L = 5;% domény délka v metrech
N = 505;% # prostorových vzorků v doméně
Ledek = 800;% # iterací vykonávat
fs = 300e6;% Zdroj frekvence v Hz
ds = L / N,%, prostorové krok v metrech
dt = ds/300e6;% "magic časovým krokem"
eps0 = 8.854e-12,% permitivita volného prostoru
mu0 = pi * 4e-7;% propustnost volného prostoru
x = linspace (0, L, N);% x souřadnice prostorových vzorků
showWKB = 0;% if = 1 pak show WKB appromination na konci

stupnice faktory% pro E a H
ae = ones (N, 1) * dt / (DS * eps0);
AM = ones (N, 1) * dt / (DS * mu0);
as = ones (N, 1);
epsr = ones (N, 1);
Mury = ones (N, 1);
sigma = nuly (N, 1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%
% Zde specifikovat epsilon, sigma, a MU profily.Jsem
% Předdefinovaných několik zajímavých příkladů.Zkuste profile = 1,2,3,4,5,6 v pořadí.
% Lze definovat epsr (i), Mur (i) (relativní permitivita a permeabilita)
% A sigma (i) (vodivost), které mají být cokoli chcete.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%
profile = 1;
pro i = 1: N
epsr (i) = 1;
Mur (i) = 1;
w1 = 0,5;
W2 = 1,5;
if (profil == 1)% dielektrické okno
if (abs (x (i)-L / 2) <0,5) epsr (i) = 4; konec
konec
if (profil == 2)% dielektrickou okno s hladký přechod
if (abs (x (i)-L / 2) <1,5) epsr (i) = 1 3 * (1 cos (pi * (abs (x (i)-L / 2)-W1) / (W2 -W1))) / 2; konec
if (abs (x (i)-L / 2) <0,5) epsr (i) = 4; konec
konec
if (profil == 3)% dielektrické diskontinuitu
if (x (i)> L / 2) epsr (i) = 9; konec
konec
if (profil == 4)% dielektrickou disontinuity s 1/4-wave odpovídající vrstvou
if (x (i)> (L/2-0.1443)) epsr (i) = 3; konec
if (x (i)> L / 2) epsr (i) = 9; konec
konec
if (profil == 5)% dirigování poloviny prostoru
if (x (i)> L / 2) Sigma (i) = 0,005; konec
konec
if (profil == 6)% sinusový dielektrické
epsr (i) = 1 sin (2 * pi * x (i) / L) ^ 2;
konec
konec
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%AE = ae. / epsr;
AM = pozměňovací návrh / Mur;
ae = ae. / (1 dt * (sigma. / epsr) / (2 * eps0));
as = (1-dt * (sigma. / epsr) / (2 * eps0 ))./( 1 dt * (sigma. / epsr) / (2 * eps0));

% Spiknutí permitivita, permeabilita, vodivost a profily
Obrázek (1)
subplot (3,1,1);
plot (x, epsr);
sítě o;
osy ([3 * DS L min (epsr) * 0.9 max (epsr) * 1.1]);
Název ('relativní permitivity');
subplot (3,1,2);
plot (x, Mur);
sítě o;
osy ([3 * DS L min (Mur) * 0.9 max (Mur) * 1.1]);
Název ('relativní permeabiliity');
subplot (3,1,3);
plot (x, sigma);
sítě o;
osy ([3 * DS L min (sigma) * 0,9 až 0,001 max (sigma) * 1.1 0.001]);
Název ('vodivost');

% Inicializovat pole na nulu
Hz = nuly (N, 1);
Ey = nuly (N, 1);
Obrázek (2);
set (gcf, 'Doublebuffer', 'o');% stanovené dvojitou vyrovnávací paměť o pro hladší grafiku
plot (EY);
sítě o;

ITER = 1: ledek
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%
% Během příštích 10 nebo tak řádků kódu tam, kde jsme skutečně integrovat Maxwellovy
% Rovnice.Všechny ostatní program je v podstatě účetnictví a kreslení.
% "Hladký zapnout" sinusový zdroj
Ey (3) = Ey (3) 2 * (1-exp (- ((ITER-1) / 50) ^ 2)) * sin (2 * pi * fs * dt * iter);
Hz (1) = Hz (2);% absorpční okrajové podmínky pro odešel-rozmnožovací vlny
pro i = 2: N-1% update pole H
Hz (i) = Hz (i)-já (i) * (Ey (i 1)-Ey (i));
konec
Ey (N) = Ey (N-1);% absorpční okrajové podmínky pro nárok rozmnožovací vlny
pro i = 2: N-1% update E pole
Ey (i) = v (i) * Ey (i)-ae (i) * (Hz (i)-Hz (i-1));
konec
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%
Obrázek (2)
odrazit
plot (x, Ey, 'b');
osy ([3 * ds L -2 2]);
sítě o;
Název ('E (modrá) a 377 * H (červená)');
vydrž
plot (x, 377 Hz *, 'r');
xlabel ('x (m)');
pauza (0);
ITER
konec

% WKB predikce pole
if (showWKB == 1)
fáze = cumsum ((epsr). ^ 0,5) * ds;
beta0 = 2 * pi FS / (300e6);
Teorie = sin (2 * pi * fs * (ledek 4) * DT-beta0 * fáze). / (epsr. ^ 0,25);
vstup ('stisknutím klávesy ENTER show WKB teorií');
plot (x, teorie, 'b.');
Teorie = sin (2 * pi * fs * (ledek 4) * DT-beta0 * fáze) .* (epsr. ^ 0,25);
plot (x, teorie, 'r.');
Název ('E (modrá), 377 * H (červená), WKB teorie (bodů)');
konec