c dotaz

[fatma1000]

toto je kód pro matlab upml
Potřebuji převést jej do c
d = 8;% hloubky PML regionu (d v počtu buněk)
NXA% = 100 x = 50 NYA původní výpočetní domény
NXA = 100; NYA = 50;
NX = NXA 2 * d,% Y Dimension
NY = NYA 2 * d;% x rozměr
EY (1: NX, 1: NY) = 0;% Ey složka
DY (1: NX, 1: NY) = 0;% Dy složka
DY1 (1: NX, 1: NY) = 0;% předchozího časového kroku
EX (1: NX, 1: NY) = 0;% Ex složka
DX (1: NX, 1: NY) = 0;% Dx složka
DX1 (1: NX, 1: NY) = 0;% předchozího časového kroku
HZ (NX, NY) = 0;% Hz složka
BZ (NX, NY) = 0;% BZ složka
BZ1 (NX, NY) = 0;% předchozího časového krokuje = NX / 2; js = NY / 2;cc = 2.99792458e8;
m0 = 4.0 * pi * 1.0e-7;% permability volného prostoru
E0 = 1,0 / (cc * cc * m0);% permitivita volného prostoru
dx = 1.5E-2; dy = 1.5E-2;% velikosti buněk
dt = dx / (2.0 * cc);% času velikost kroku

% c = 1/sqrt (m0 * E0)% rychlosti světla
N = 700;% času počet kroků
H (N) = 0;% excitační funkce
W (N) = 0;% pozorovala bod
% Koeficientu matice
cadx (1: NX, 1: NY) = 1;
cbdx (1: NX, 1: NY) = dt;
Cady (1: NX, 1: NY) = 1;
cbdy (1: NX, 1: NY) = dt;
caex (1: NX, 1: NY) = 1;
cbex (1: NX, 1: NY) = 1/e0;
ccex (1: NX, 1: NY) = 1/e0;
caey (1: NX, 1: NY) = 1;
cbey (1: NX, 1: NY) = 1/e0;
ccey (1: NX, 1: NY) = 1/e0;
cabz (1: NX, 1: NY) = 1;
cbbz (1: NX, 1: NY) = dt;
cahz (1: NX, 1: NY) = 1;
cbhz (1: NX, 1: NY) = 1/m0;
cchz (1: NX, 1: NY) = 1/m0;
sigma_x (1: NX, 1: NY) = 0;
sigma_y (1: NX, 1: NY) = 0;

% zdroj
pro n = 1: N
Pokud n * dt <1E-9
H = (1 / 320) * (10-15 * cos (2 * pi * 1e9 * dt * n) 6 * cos (4 * pi * 1e9 * dt * n) ...
-cos (6 * pi * 1e9 * dt * n));
konec
konec
% caculation pro coefficints v PML
m = 2;% variace (pořadí)
R = 1E-5;% reflecion koeficient

sm =- (sqrt (e0/m0 ))*(( m 1) * log (R) / (2 * d * dx));

pro i = 1: d
sigma (i) = sm * ((i-1) / d) ^ m;
sigmas (i) = sm * ((i-0.5) / d) ^ m;
konec
%
PEC% po vrstvě PML
EX (1: NX, NY) = 0;% PEC
EX (1: NX, 1) = 0;% PEC
EY (1,1: NY) = 0;% PEC
EY (NX, 1: NY) = 0;% PECpro ii = 1: NX
pro jj = 1: NY

sigma_x (ii, jj) = 0;
sigma_y (ii, jj) = 0;
sigmas_x (ii, jj) = 0;
sigmas_y (ii, jj) = 0;
if (ii <= d)
sigma_x (ii, jj) = sigma (d-ii 1);
sigmas_x (ii, jj) = sigmas (d-ii 1);
konec
if (ii> NX-d)
sigma_x (ii, jj) = sigma (ii-NX d);
sigmas_x (ii, jj) = sigmas (ii-NX d);
konecif (jj <= d)
sigma_y (ii, jj) = sigma (d-jj 1);
sigmas_y (ii, jj) = sigmas (d-jj 1);
konec
if (jj> NY-d)
sigma_y (ii, jj) = sigma (jj-NY d);
sigmas_y (ii, jj) = sigmas (jj-NY d);
konec
cbex (ii, jj) = (2 * e0 sigma_x (ii, jj) * dt) / (2 * e0 * E0);
ccex (ii, jj) = (2 * E0-sigma_x (ii, jj) * dt) / (2 * e0 * E0);
caey (ii, jj) = (2 * E0-sigma_x (ii, jj) * dt) / (2 * e0 sigma_x (ii, jj) * dt);
cbey (ii, jj) = (2 * e0 sigma_y (ii, jj) * dt) / (2 * e0 * e0 sigma_x (ii, jj) * e0 * dt);
ccey (ii, jj) = (2 * E0-sigma_y (ii, jj) * dt) / (2 * e0 * e0 sigma_x (ii, jj) * e0 * dt);
cadx (ii, jj) = (2 * E0-sigma_y (ii, jj) * dt) / (2 * e0 sigma_y (ii, jj) * dt);
cbdx (ii, jj) = 2 * e0 * dt / (2 * e0 sigma_y (ii, jj) * dt);
cahz (ii, jj) = (2 * E0-sigmas_y (ii, jj) * dt) / (2 * e0 sigma_y (ii, jj) * dt);
cbhz (ii, jj) = (2 * E0) / (2 * e0 sigmas_y (ii, jj) * dt) / m0;
cchz (ii, jj) = (2 * E0) / (2 * e0 sigmas_y (ii, jj) * dt) / m0;
cabz (ii, jj) = (2 * E0-sigmas_x (ii, jj) * dt) / (2 * e0 sigmas_x (ii, jj) * dt);
cbbz (ii, jj) = (2 * e0 * dt) / (2 * e0 sigmas_x (ii, jj) * dt);

konec
konecpro n = 1: N

DX1 = DX;
DY1 = DY;
pro ii = 1: NX
pro jj = 1: NY
if (jj> 1)
DX (ii, jj) = cadx (ii, jj) * DX (ii, jj) cbdx (ii, jj) * (HZ (ii, jj)-HZ (ii, jj-1)) / dy;
jiné
DX (ii, jj) = 0;
konec
konec
konecpro ii = 1: NX
pro jj = 1: NY
if (ii> 1)
DY (ii, jj) = Cady (ii, jj) * DY (ii, jj)-cbdy (ii, jj) * (HZ (ii, jj)-HZ (II-1, jj)) / dx;
jiné
DY (ii, jj) = 0;
konec

konec
konec

pro ii = 1: NX
pro jj = 1: NY
EX (ii, jj) = caex (ii, jj) * EX (ii, jj) cbex (ii, jj) * DX (ii, jj)-ccex (ii, jj) * DX1 (ii, jj);

konec
konec
pro ii = 1: NX
pro jj = 1: NY
EY (ii, jj) = caey (ii, jj) * EY (ii, jj) cbey (ii, jj) * DY (ii, jj)-ccey (ii, jj) * DY1 (ii, jj);
konec
konecBZ1 = BZ;
pro ii = 1: NX-1
pro jj = 1: NY-1
BZ (ii, jj) = cabz (ii, jj) * BZ (ii, jj)-cbbz (ii, jj) * ((EY (ii 1, jj)-EY (ii, jj)) / dx-( EX (ii, jj 1)-EX (ii, jj)) / dy);
konec
konec

pro ii = 1: NX-1
pro jj = 1: NY-1
HZ (ii, jj) = cahz (ii, jj) * HZ (ii, jj) cbhz (ii, jj) * BZ (ii, jj)-cchz (ii, jj) * BZ1 (ii, jj);
konec
konec
HZ (je, js) = HZ (je, js) H ;
% W = HZ (tj. D 2);
W = HZ (tj. D 2);

konec
n = 1: N;
plot (n, W )