APENDIX A PERSAMAAN MAXWELL
Persamaan Maxwell dalam bentuk differensial. ⃗
Hukum Gauss
⃗ ⃗
Hukum Gauss untuk magnetisme
⃗ ⃗
Hukum Induksi Faraday
⃗
⃗
Hukum Ampere
⃗
⃗
(A.1) (A.2) ⃗
(A.3) ⃗
(A.4)
Persamaan Maxwell dalam bentuk integral Hukum Gauss
∮ ⃗
Hukum Gauss untuk magnetisme
∮ ⃗
Hukum Induksi Faraday
∮⃗
Hukum Ampere
∮⃗
⃗
∫
(A.5) (A.6)
∫⃗ ∫⃗
(A.7) ∫
(A.8)
Karena (teorema Gauss dan Stokes) ∮ ⃗
∮
∮ ⃗
∮
Untuk medan listrik pada persamaan (A.5) dan (A.7) dalam keadaan konservatif sehingga dapat ditulis :
⃗
⃗ ⃗
(A.9)
⃗
(A.10)
Selanjutnya hubungan vektornya ⃗
⃗
(A.11)
⃗
⃗ ⃗
(A.12)
dan
i
⃗
∫⃗
(A.13)
⃗ = kerapatan fluks listrik (C/m2) = rapat volume muatan (C/m3) ⃗ = Medan Listrik (V/m) ⃗ = tegangan listrik (V) = permitivitas dielektrik (F/m)
Hubungkan persamaan (A.9), (A.11) dan (A.12) akan didapat persamaan Poisson : ⃗
⃗⃗
(A.14a)
Jika konstan, ⃗ Jika
(A.14b)
maka menjadi persamaan Laplace : ⃗
⃗⃗
(A.15a)
⃗
(A.15b)
Karena (hukum Biot-Savart) ∮ ⃗
∮
∮⃗ ⃗ = fluks magnet (A/m) ⃗ = Medan magnet (T) Maka persamaan (A.6) dan(A.8) menjadi ⃗
⃗
(A.16)
Dan ⃗ ⃗
(A.17)
Dengan hubungan vector ⃗
= permeabilitas magnetic (H/m) =konduktivitas (mhos/m)
⃗
(A.18)
⃗
(A.19)
Medan magnet pada potensial A ⃗
⃗
(A.20)
Dengan menggunakan vektor identitas ⃗
(⃗
⃗ (⃗
)
)
Jika (A.17) dan (A.18) dihubungkan dalam kondisi tanpa arus (ruang hampa) maka ⃗
sehingga didapat persamaan Poisson : (A.21)
Ketika
maka (A.22)
i
APENDIX B PERSAMAAN GELOMBANG
Gelombang Elektromagnetik merambat secara transversal. Dengan medan listrik dan medan magnet yang merambat saling tegak lurus. Misalkan gelombang elektromagnet merambat ke sumbu-x dan arah gerak medan listrik ke sumbu-y sementara medan magnetik ke sumbu-z sehingga diperoleh ⃗ ⃗
⃗
⃗ ⃗
dan ⃗
⃗
⃗ . Maka persamaan Maxwell akan berubah menjadi : ⃗ ⃗ ⃗
(B.1) ⃗ ⃗
⃗
(B.2)
⃗
⃗
⃗ ⃗
(B.3)
⃗
⃗
⃗
(B.4) ⃗
⃗ ⃗
(B.5) ⃗⃗
⃗
(B.6)
Persamaan (B.4) dan persamaan (B.6) bisa disimpulkan medan listrik dan medan magnet hanya berpengaruh pada sumbu x dan gayut waktu. ⃗
⃗
⃗
⃗
Jika digabungkan maka menjadi : ⃗
⃗
Bandingkan dengan persamaan gelombang umum :
Atau ,
=0
Maka diperoleh nilai √ Dengan v = c = kecepatan cahaya = 3 x 108 m Dengan mengubah daerah waktu menjadi frekuensi maka
Dengan Maka persawaan Maxwell dapat dituliskan ⃗( )
( )⃗ ( )
⃗( )
( )⃗ ( )
Dengan vektor perambatan elektromagnetik sebagai berikut : ( )
⃗( )
⃗( )
( )
⃗( )
⃗( )
⃗( )
⃗( )
( )
( )⃗ ( )
(
( )
( )
) ( ) ⃗ ( )⃗ ( )
( )
(
) ( )⃗ ( )
( ) ( ) ⃗ ( )⃗ ( )
Maka dapat disimpulkan bahwa √ Dengan n adalah indeks bias.
√
i
APENDIX C METODE FDTD (FINITE DIFFERENCE TIME DOMAIN)
Metode ini menggunakan turunan fungsi yang dihubungkan dengan beda hingga seperti berikut (
( )
)
( )
(C.1)
Dengan dihubungkan ke dalam medan listrik dan medan magnet, akan dihasilkan persamaan FDTDnya misalkan untuk waktu penjalaran medan E dan medan H sebagai berikut untuk ⃗ (
untuk ⃗
)
dan letak pada sumbu z misalkan untuk ⃗ (
untuk ⃗
)
Telah diketahui persamaan Maxwell
Dengan notasi ⃗ ( ) )
⃗(
⃗
⃗
⃗
⃗
(C.2) (C.3)
) dan ⃗
(
)
⃗ ((
)
Maka didapat : ⃗
⃗ ( ⃗
⃗
⃗ ( )
) ( )
⃗ ( )
⃗
⃗
(
)
⃗
(
)
⃗
⃗
(
)
⃗
(
)
)
(
Dengan menghubungkan ke dalam persamaan (C.2) dan (C.3) persamaan diatas menjadi ⃗
( ⃗
⃗
)
(
⃗ ( )
( )
(⃗ (
) (⃗
( ⃗
Dengan
dan
⃗
)
⃗ ( )
) (
)
(C.4)
⃗
(C.5)
dalam model Drude diberikan dalam persamaan ( )
Substitusi model Drude dalam persamaan (C.5) akan didapat ⃗ ⃗ ⃗ ⃗
)⃗
( (
(
⃗
⃗
(
)
( )
(
(
)⃗
(
)⃗
)
)⃗
)
⃗
⃗
)⃗
(
(C.6)
nilai fungsi frekuensi dikonversikan menjadi fungsi waktu dengan hubungan dan
. Dengan perkalian ⃗
merupakan waktu rata-rata. Maka :
⃗
⃗
⃗
⃗
Dengan cara yang sama sebelumnya akan di dapat : ⃗
(⃗ ( ⃗
⃗ ⃗
⃗ )
(⃗
) (⃗
⃗
( ⃗
)
)
⃗
(C.7)
Nilai variabel tambahan dalam persamaan yang baru merupakan skalar yaitu (
)
(
)
(
)
i
( (
)
)
)
dan ⃗
(⃗
⃗
⃗
)
( ⃗
⃗
)
(⃗
(⃗
⃗
( ⃗
)
)
⃗
(
)
)
Maka persamaan FDTD gelombang yang merambat dalam arah z sebagai berikut ⃗
⃗ ( )
( )
(⃗ ( )
⃗ (
))
(
)
dan ⃗
⃗ ( )
( )
(⃗ ( )
⃗ (
))
(
)
Persamaan diatas memberikan algoritma FDTD untuk ⃗ dan ⃗ . Untuk indeks bias FDTD dihubungkan dengan transformasi Fourier sebagai pengubah ke FDTD : √ Untuk medium yang homogen nilai = 1 sehingga indeks bias hanya bergantung dan bergantung kepada medan listrik. ̃(
̃(
∫ ⃗(
)
)
∫ ⃗( )
)
(
(
)
)
(
)
Sehingga didapat bahwa hasil transformasi Fouriernya adalah ̃(
)
⃗( )
(
)
Indeks bias merupakan perbandingan kecepatan antara dua medium maka persamaan diatas dapat diubah menjadi : ̃(
)
̃(
)
|
( (
⃗( ) ⃗( ) ) | )
(
(
)
(
) )
(
|
)
⃗ ( ⃗ (
) )
|
Untuk simulasi dalam jubah silinder maka diperlukan persamaan dalam koordinat silinder sebagai berikut: ( ) ( ) ( )
(
)
(
)
(
)
(
)
Dari transformasi kartesius : (C.13) (
)
(C.14) (C.15)
Sehingga ⃗
( ⃗ )
(
⃗
)
(
)
(
)
(
)
⃗
( )
⃗ ⃗
( )
⃗
i
APENDIX D PROGRAM MATLAB clc clear all % Adimas Agung % 110801001 % Simulasi Parameter. SIZE = 4*1024; % Angka beda langkah SlabLeft = round(SIZE/3); % Letak sebelah kiri lempengan SlabRight = round(2*SIZE/3); % Letak sebelah kanan lempengan MaxTime = 4*SIZE; % Jumlah langkah waktu PulseWidth = round(SIZE/8); % Lebar Pulsa Gauss td = PulseWidth; % tempo pulsa. source = 10; % Letak Sumber SaveFields = 1; % 0. No, 1. Ya menyimpan. SnapshotInterval = 32; % jumlah waktu untuk memfoto (menangkap layar).
% Memilih Sumber Gelombang. % 1. Gaussian 2. Sine wave 3. Ricker wavelet SourceChoice = 2;
% Konstanta yang digunakan. c = 3e8; pi = 3.141592654; e0 = (1e-9)/(36*pi); u0 = (1e-7)*4*pi;
dt = 0.5e-11; dz = 3e-3; Sc = c * dt/dz
l = PulseWidth*dz; f = c/l fmax = 1/(2*dt) w = 2*pi*f; k0 = w/c; % bilangan gelombang.
% Parameter dalam Ricker wavelet. if SourceChoice == 3 fp = f; % puncak ke puncak frekuensi dr = PulseWidth*dt*2; % Waktu tunggu end
% Inisialisasi. Ex = zeros(SIZE, 3); % komponen x untuk medan E Dx = zeros(SIZE, 3); % komponen x untuk D Hy = zeros(SIZE, 3); % komponen y untuk medan H By = zeros(SIZE, 3); % komponen y medan B
% Medan datang dan transmisi. Exi = zeros(MaxTime, 1); Ext = zeros(MaxTime, 1); Extt = zeros(MaxTime, 1); x1 = SlabLeft+1; % letak observasi (posisi pengamatan).
% Penghitungan index Bias. Z1 = SlabLeft + 50; z1 = Z1*dz; Z2 = SlabLeft + 60; z2 = Z2*dz; Exz1 = zeros(MaxTime, 1); Exz2 = zeros(MaxTime, 1);
einf = ones(SIZE,1); einf(SlabLeft:SlabRight) = 1; % einf(Drude) atau er di lempengan. uinf = ones(SIZE,1); uinf(SlabLeft:SlabRight) = 1; % uinf(Drude) atau ur di lempengan. wpesq = zeros(SIZE,1); wpesq(SlabLeft:SlabRight) = 2*w^2; % nilai DNG(Drude) untuk wpe kuadrat pada lempengan. wpmsq = zeros(SIZE,1); wpmsq(SlabLeft:SlabRight) = 2*w^2; % nilai DNG(Drude) untuk wpm kuadrat pada lempengan. ge = zeros(SIZE,1); ge(SlabLeft:SlabRight) = w/32; % frekuensi listrik pada lempengan. gm = zeros(SIZE,1);
i
gm(SlabLeft:SlabRight) = w/32; % frekuensi magnet pada lempengan.
a0 = (4*dt^2)./(e0*(4*einf+dt^2*wpesq+2*dt*einf.*ge)); a = (1/dt^2)*a0; b = (1/(2*dt))*ge.*a0; c = (e0/dt^2)*einf.*a0; d = (-1*e0/4).*wpesq.*a0; e = (1/(2*dt))*e0*einf.*ge.*a0; am0 = (4*dt^2)./(u0*(4*uinf+dt^2*wpmsq+2*dt*uinf.*gm)); am = (1/dt^2)*am0; bm = (1/(2*dt))*gm.*am0; cm = (u0/dt^2)*uinf.*am0; dm = (-1*u0/4).*wpmsq.*am0; em = (1/(2*dt))*u0*uinf.*gm.*am0;
if SaveFields == 1 ExSnapshots = zeros(SIZE, MaxTime/SnapshotInterval); % Data yang akan diplot (bentuk kurva). frame = 1; end
n1 = 1; n2 = 2; linecount = 0; % waktu keluar dari pengulangan. tic % uji pengulangan untuk medan datang pada ruang bebas. for q = 0:MaxTime % Penghitungan Hy menggunakan turunan persamaan Hy. Beda waktu q. Hy(1:SIZE-1,n2) = Hy(1:SIZE-1,n1) + ( ( Ex(1:SIZE-1,n1) Ex(2:SIZE,n1) ) * dt/(u0*dz) );
% ABC untuk ukuran H. Hy(SIZE,n2) = Hy(SIZE-1,n1) + (Sc-1)/(Sc+1)*(Hy(SIZE-1,n2) Hy(SIZE,n1) );
% Penghitungan Ex menggunakan turunan persamaan Ex. Beda waktu q+1/2.
Ex(2:SIZE,n2) = Ex(2:SIZE, n1) + ( dt/(e0*dz)*(Hy(1:SIZE-1, n2) Hy(2:SIZE, n2)) );
% ABC untuk E di 1. Ex(1,n2) = Ex(2,n1) + (Sc-1)/(Sc+1)*(Ex(2,n2) - Ex(2,n1));
% Sumber. if SourceChoice == 1 Ex(source,n2) = Ex(source,n2) + exp( -1*((q-td)/(PulseWidth/4))^2 ) * Sc; elseif SourceChoice == 2 Ex(source,n2) = Ex(source,n2) + sin(2*pi*f*(q)*dt) * Sc; elseif SourceChoice == 3 Ex(source,n2) = Ex(source,n2) + (1-2*(pi*fp*(q*dt-dr))^2)*exp(1*(pi*fp*(q*dt-dr))^2) * Sc; end
Exi(q+1) = Ex(x1,n2); % Medan datang di kiri lempengan.
temp = n1; n1 = n2; n2 = temp; end % Pengulangan inisialisasi medan untuk simulasi. Ex = zeros(SIZE, 3); % komponen x untuk medan E Hy = zeros(SIZE, 3); % komponen y untuk medan H % Simulasi dimulai. fprintf ( 1, 'Simulation started... \n'); for q = 0:MaxTime
% Indikator. fprintf(1, repmat('\b',1,linecount)); linecount = fprintf(1, '%g %%', (q*100)/MaxTime );
Ex(:,3) = Ex(:,n2); Dx(:,3) = Dx(:,n2); Hy(:,3) = Hy(:,n2); By(:,3) = By(:,n2);
i
% Penghitungan Hy menggunakan turunan persamaan Hy. Beda waktu q. By(1:SIZE-1,n2) = By(1:SIZE-1,n1) + ( ( Ex(1:SIZE-1,n1) Ex(2:SIZE,n1) ) * dt/(dz) ); Hy(:,n2) = am.*(By(:,n2)-2*By(:,n1)+By(:,3))+bm.*(By(:,n2)By(:,3))+cm.*(2*Hy(:,n1)Hy(:,3))+dm.*(2*Hy(:,n1)+Hy(:,3))+em.*(Hy(:,3));
% ABC untuk ukuran. Hy(SIZE,n2) = Hy(SIZE-1,n1) + (Sc-1)/(Sc+1)*(Hy(SIZE-1,n2) Hy(SIZE,n1) ); By(SIZE,n2) = u0*Hy(SIZE,n2);
% Penghitungan Ex menggunakan persamaan diferensial Ex. Beda waktu q+1/2. Dx(2:SIZE,n2) = Dx(2:SIZE, n1) + ( dt/(dz)*(Hy(1:SIZE-1, n2) Hy(2:SIZE, n2)) ); Ex(:,n2) = a.*(Dx(:,n2)-2*Dx(:,n1)+Dx(:,3))+b.*(Dx(:,n2)Dx(:,3))+c.*(2*Ex(:,n1)Ex(:,3))+d.*(2*Ex(:,n1)+Ex(:,3))+e.*(Ex(:,3)); % ABC untuk E di 1. Ex(1,n2) = Ex(2,n1) + (Sc-1)/(Sc+1)*(Ex(2,n2) - Ex(2,n1)); Dx(1,n2) = e0*Ex(1,n2);
% Sumber. if SourceChoice == 1 Ex(source,n2) = Ex(source,n2) + exp( -1*((q-td)/(PulseWidth/4))^2 ) * Sc; elseif SourceChoice == 2 Ex(source,n2) = Ex(source,n2) + sin(2*pi*f*(q)*dt) * Sc; elseif SourceChoice == 3 Ex(source,n2) = Ex(source,n2) + (1-2*(pi*fp*(q*dt-dr))^2)*exp(1*(pi*fp*(q*dt-dr))^2) * Sc; end Dx(source,n2) = e0*Ex(source,n2);
if (SaveFields == 1 && mod(q,SnapshotInterval) == 0) ExSnapshots(:,frame) = Ex(:,n2);
frame=frame+1; end
Ext(q+1) = Ex(x1,n2); Extt(q+1) = Ex(SlabRight+10,n2);
% Medan untuk Penghitungan index bias. Exz1(q+1) = Ex(Z1, n2); Exz2(q+1) = Ex(Z2, n2);
temp = n1; n1 = n2; n2 = temp; end fprintf ( 1, '\nSimulation complete! \n'); toc % Pengolahan. Fs = 1/dt;
% pengambilan sampel frekuensi
T = dt;
% waktu sampel
L = length(Exi);
% Panjang sinyal
t = (0:L-1)*T;
% Vektor waktu
fspan = 100;
% Titik alur dalam daerah frekuensi
figure(1) subplot(211) plot(Fs*t, Exi, 'LineWidth', 2.0, 'Color', 'b') set(gca, 'FontSize', 10, 'FontWeight', 'b') title('Medan Listrik Datang', 'FontSize', 12, 'FontWeight', 'b') xlabel('waktu', 'FontSize', 11, 'FontWeight', 'b') grid on figure(2) subplot(211) plot(Fs*t, Ext, 'LineWidth', 2.0, 'Color', 'b') set(gca, 'FontSize', 10, 'FontWeight', 'b') title('Transmisi Medan Listrik', 'FontSize', 12, 'FontWeight', 'b') xlabel('waktu', 'FontSize', 11, 'FontWeight', 'b') grid on figure(3) subplot(211)
i
plot(Fs*t, Extt, 'LineWidth', 2.0, 'Color', 'b') set(gca, 'FontSize', 10, 'FontWeight', 'b') title('Transmisi Medan Listrik Melalui Lempeng', 'FontSize', 12, 'FontWeight', 'b') xlabel('waktu', 'FontSize', 11, 'FontWeight', 'b') grid on
NFFT = 2^nextpow2(L); % untuk tenaga kedua dari Exi % Medan datang dan transmisi. EXI = fft(Exi,NFFT)/L; EXT = fft(Ext,NFFT)/L; EXTT = fft(Extt,NFFT)/L; % Penghitungan Index Bias. EXZ1 = fft(Exz1,NFFT)/L; EXZ2 = fft(Exz2,NFFT)/L; f = Fs/2*linspace(0,1,NFFT/2+1); % Membuat alur spektrum amplitudo satu bagian. figure(1) subplot(212) EXIp = 2*abs(EXI(1:NFFT/2+1)); plot(f(1:fspan), EXIp(1:fspan), 'LineWidth', 2.0, 'Color', 'r') set(gca, 'FontSize', 10, 'FontWeight', 'b') title('Spektrum Amplitudo Exi(t)', 'FontSize', 12, 'FontWeight', 'b') xlabel('Frekuensi (Hz)', 'FontSize', 11, 'FontWeight', 'b') ylabel('|EXI(f)|', 'FontSize', 11, 'FontWeight', 'b') grid on figure(2) subplot(212) EXTp = 2*abs(EXT(1:NFFT/2+1)); plot(f(1:fspan), EXTp(1:fspan), 'LineWidth', 2.0, 'Color', 'r') set(gca, 'FontSize', 10, 'FontWeight', 'b') title('Spektrum Amplitudo Ext(t)', 'FontSize', 12, 'FontWeight', 'b') xlabel('Frekuensi (Hz)', 'FontSize', 11, 'FontWeight', 'b') ylabel('|EXT(f)|', 'FontSize', 11, 'FontWeight', 'b') grid on figure(3) subplot(212) EXTTp = 2*abs(EXTT(1:NFFT/2+1)); plot(f(1:fspan), EXTTp(1:fspan), 'LineWidth', 2.0, 'Color', 'r')
set(gca, 'FontSize', 10, 'FontWeight', 'b') title('Spektrum Amplitudo Extt(t)', 'FontSize', 12, 'FontWeight', 'b') xlabel('Frekuensi(Hz)', 'FontSize', 11, 'FontWeight', 'b') ylabel('|EXT(f)|', 'FontSize', 11, 'FontWeight', 'b') grid on % Koefisien Transmisi. figure(4) subplot(211) TAU = abs(EXT(1:NFFT/2+1)./EXI(1:NFFT/2+1)); plot(f(1:fspan), TAU(1:fspan), 'LineWidth', 2.0, 'Color', 'b') set(gca, 'FontSize', 10, 'FontWeight', 'b') title('Koefisien Transmisi', 'FontSize', 12, 'FontWeight', 'b') xlabel('Frekuensi (Hz)', 'FontSize', 11, 'FontWeight', 'b') ylabel('|EXT(f)/EXI(f)|', 'FontSize', 11, 'FontWeight', 'b') axis([-1 1 -2 2]) axis 'auto x' grid on subplot(212) plot(f(1:fspan), 1-TAU(1:fspan), 'LineWidth', 2.0, 'Color', 'b') set(gca, 'FontSize', 10, 'FontWeight', 'b') title('Koefisien Refleksi', 'FontSize', 12, 'FontWeight', 'b') xlabel('Frekuensi (Hz)', 'FontSize', 11, 'FontWeight', 'b') ylabel('1-|EXT(f)/EXI(f)|', 'FontSize', 11, 'FontWeight', 'b') axis([-1 1 -2 2]) axis 'auto x' grid on % Penghitungan Indeks Bias. nFDTD = (1/(1i*k0*(z1-z2))).*log(EXZ2(1:NFFT/2+1)./EXZ1(1:NFFT/2+1)); figure(5) subplot(211) plot(f(1:fspan), real(nFDTD(1:fspan)), 'LineWidth', 2.0, 'Color', 'b'); set(gca, 'FontSize', 10, 'FontWeight', 'b') title('Indeks Bias re(n)', 'FontSize', 12, 'FontWeight', 'b') xlabel('Frekuensi (Hz)', 'FontSize', 11) ylabel('re(n)', 'FontSize', 11) grid on subplot(212)
i
plot(f(1:fspan), imag(nFDTD(1:fspan)), 'LineWidth', 2.0, 'Color', 'r'); set(gca, 'FontSize', 10, 'FontWeight', 'b') title('Indeks Bias im(n)', 'FontSize', 12, 'FontWeight', 'b') xlabel('Frekuensi (Hz)', 'FontSize', 11, 'FontWeight', 'b') ylabel('im(n)', 'FontSize', 11, 'FontWeight', 'b') grid on
if SaveFields == 1 % Simulasi Animasi. for i=1:frame-1 figure (6) % Batas penyebaran. hold off plot([SlabLeft SlabLeft], [-1 1], 'Color', 'r'); hold on plot([SlabRight SlabRight], [-1 1], 'Color', 'r'); plot(ExSnapshots(:,i), 'LineWidth', 2.0, 'Color', 'b'); set(gca, 'FontSize', 10, 'FontWeight', 'b') axis([0 SIZE -1 1]) title('Simulasi Fungsi Waktu', 'FontSize', 12, 'FontWeight', 'b') xlabel('Beda Langkah(k)', 'FontSize', 11, 'FontWeight', 'b') ylabel('Medan Listrik (Ex)', 'FontSize', 11, 'FontWeight', 'b') grid on end end
APENDIX E PROGRAM INVISIBLE CLOAK clear all; close all; % Nilai Konstanta %******************************************************************* c_0 = 3.0e8; % Kecepatan cahaya mu_0 = 4.0*pi*1.0e-7; % Permeabilitas eps_0 = 8.8542e-12; % Permitivitas losstangent = 0.0; ie = 401; je = 400;
% # Besar sel dalam arah x % # Besar sel dalam arah y
ib = ie + 1; jb = je + 1; npmls = 20; ip = ie - npmls; jp = je - npmls; xc = round(ie/2); yc = round(je/2); is = 25; js = 25;
% Letak Sumber Gelombang
it = 30; jt = 30;
% Daerah Medan atas atau bawah % Daerah medan kiri atau kanan
freq = 6.0e14; omega = 2*pi*freq; k_0 = omega/c_0; dx = c_0/freq/400; dt = dx/c_0/sqrt(2);
% Beda waktu
R1 = c_0/freq/8/dx R2 = 1.43*R1;
% dimensi Jubah
nmax = 400000;
% Jumlah beda waktu
aimp = sqrt(mu_0/eps_0);
% impedansi gelombang
threshold = 0.1; step = ceil(nmax/10);
% FDTD error
% freq = 2.0e9; % omega = 2*pi*freq; tau = 1/freq; delay = 3*tau/dt;
i
N = round(tau/dt); M = round(c_0/freq/dx); source = zeros(1,nmax); clear j; ST = 300; st = 20e3; for n=1:nmax if n < ST*N x = 1.0 - (ST*N-n)/(ST*N); g = 10.0*x^3 - 15.0*x^4 + 6.0*x^5; % sumber(n) = g * sin(2*pi*freq*n*dt); source(n) = g * exp(j*2*pi*freq*n*dt); else % sumber(n) = sin(2*pi*freq*n*dt); source(n) = exp(j*2*pi*freq*n*dt); end % sumber (n) = exp(-(((n-delay)*dt)/tau)^2) .* exp(j*2*pi*freq*n*dt); end %******************************************************************* % Inisialisasi matriks untuk komponen medan %******************************************************************* Dx = zeros(ie,jb); Dx_h1 = zeros(ie,jb); Dx_h2 = zeros(ie,jb); Dxsynch = zeros(ie,jb); Dx_h1synch = zeros(ie,jb); Dx_h2synch = zeros(ie,jb); caDx = ones(ie,jb); cbDx = ones(ie,jb).*(dt/eps_0/dx); Ex = zeros(ie,jb); Ex_h1 = zeros(ie,jb); Ex_h2 = zeros(ie,jb); Ex_h1synch = zeros(ie,jb); Ex_h2synch = zeros(ie,jb); Dy = zeros(ib,je); Dy_h1 = zeros(ib,je); Dy_h2 = zeros(ib,je); Dysynch = zeros(ib,je); Dy_h1synch = zeros(ib,je); Dy_h2synch = zeros(ib,je); caDy = ones(ib,je); cbDy = ones(ib,je).*(dt/eps_0/dx); Ey = zeros(ib,je); Ey_h1 = zeros(ib,je); Ey_h2 = zeros(ib,je); Ey_h1synch = zeros(ib,je); Ey_h2synch = zeros(ib,je); Bz = zeros(ie,je); Bz_h1 = zeros(ie,je); Bz_h2 = zeros(ie,je); Bzx = zeros(ie,je); Bzy = zeros(ie,je); Hz_p = zeros(ie,je); daBzx = ones(ie,je); dbBzx = ones(ie,je).*(dt/mu_0/dx); daBzy = ones(ie,je); dbBzy = ones(ie,je).*(dt/mu_0/dx); daBz = ones(ie,je); dbBz = ones(ie,je).*(dt/mu_0/dx); Hz = zeros(ie,je); Hz_h1 = zeros(ie,je); Hz_h2 = zeros(ie,je); % % % %
HzsourceTF = zeros(length(it:ie-it),nmax); HztransTF = zeros(length(it:ie-it),nmax); HzsourceSF = zeros(length(it:ie-it),nmax); HztransSF = zeros(length(it:ie-it),nmax);
% Coefficients for cloak a0x = ones(ie,jb); b0x = ones(ie,jb); b0xy = zeros(ie,jb); a1x = zeros(ie,jb); b1x = zeros(ie,jb); a1xy = zeros(ie,jb); b1xy = zeros(ie,jb); a2x = zeros(ie,jb); b2x = zeros(ie,jb); a2xy = zeros(ie,jb); b2xy = zeros(ie,jb); a0y = ones(ib,je); b0y = ones(ib,je); b0yx = zeros(ib,je); a1y = zeros(ib,je); b1y = zeros(ib,je); a1yx = zeros(ib,je); b1yx = zeros(ib,je);
a2y = zeros(ib,je); b2y = zeros(ib,je); a2yx = zeros(ib,je); b2yx = zeros(ib,je); a0z = ones(ie,je); b0z = ones(ie,je); a1z = zeros(ie,je); b1z = zeros(ie,je); a2z = zeros(ie,je); b2z = zeros(ie,je); % Medan Datang Ex_inc = zeros(1,jb); Hzx_inc = zeros(1,jb); Ey_inc = zeros(1,ib); Hzy_inc = zeros(1,ib); caEx_inc = ones(1,jb); cbEx_inc = ones(1,jb).*(dt/eps_0/dx); caEy_inc = ones(1,ib); cbEy_inc = ones(1,ib).*(dt/eps_0/dx); daHzx_inc = ones(1,je); dbHzx_inc = ones(1,je).*(dt/mu_0/dx); daHzy_inc = ones(1,ie); dbHzy_inc = ones(1,ie).*(dt/mu_0/dx); %******************************************************************* % Jubah %******************************************************************* for i = 1 : ie for j = 1 : je r = sqrt((i - xc + 0.5).^2 + (j - yc).^2); if (r<=R2) && (r>=R1); erx = ((r-R1) / (r)); ethetax = ((r) / (r-R1)); sigma_px = ethetax*losstangent*((tan((omega*dt)/2))/(dt/2)); gamma_px = (2*losstangent*erx*sin((omega*dt)/2))/((1erx)*dt*cos((omega*dt)/2)); omega_px = sqrt((2*sin((omega*dt)/2)*(-2*(erx1)*sin((omega*dt)/2)+... losstangent*erx*gamma_px*dt*cos((omega*dt)/2)))/((dt^2)*((cos((omega* dt)/2))^2)));
sinx = (j - yc)/r; cosx = (i - xc + 0.5)/r; ax = ((cosx.^2) + (ethetax.*(sinx.^2))) ./ (dt^2); bx = ((gamma_px.*(cosx.^2)) + ((sigma_px + (ethetax*gamma_px)).*(sinx.^2))) ./ (2*dt); cx = (((cosx.^2).*((omega_px)^2)) + (sigma_px*gamma_px.*(sinx.^2))) ./ (4); wx = ((1-ethetax).*sinx.*cosx) ./ (dt^2); fx = ((gamma_px-sigma_px-(ethetax*gamma_px)).*sinx.*cosx) ./ (2*dt); vx = ((((omega_px)^2)-(sigma_px*gamma_px)).*sinx.*cosx) ./ (4); kx = 1 / (dt^2); lx = gamma_px / (2*dt);
i
tx = ((sinx.^2) + (ethetax.*(cosx.^2))) ./ (dt^2); qx = ((gamma_px.*(sinx.^2)) + ((sigma_px + (ethetax*gamma_px)).*(cosx.^2))) ./ (2*dt); px = (((sinx.^2).*((omega_px)^2)) + (sigma_px*gamma_px.*(cosx.^2))) ./ (4); Ax1 Bx1 Cx1 Dx1
= = = =
ax+bx+cx; Ax2 = (2.*ax)-(2.*cx); Ax3 = ax-bx+cx; wx+fx+vx; Bx2 = (2.*wx)-(2.*vx); Bx3 = wx-fx+vx; kx+lx; Cx2 = (2.*kx); Cx3 = kx-lx; tx+qx+px; Dx2 = (2.*tx)-(2.*px); Dx3 = tx-qx+px;
a0x(i, j) = Ax1 - (((Bx1).^2) ./ Dx1); a1x(i, j) = Ax2 - ((Bx1.*Bx2) ./ Dx1); a2x(i, j) = ((Bx1.*Bx3) ./ Dx1) - (Ax3);
a1xy(i, j) = Bx2 - ((Bx1.*Dx2) ./ Dx1); a2xy(i, j) = ((Bx1.*Dx3) ./ Dx1) - (Bx3);
b0x(i, j) = Cx1; b1x(i, j) = -Cx2; b2x(i, j) = Cx3;
b0xy(i, j) = -(Bx1.*Cx1) ./ Dx1; b1xy(i, j) = (Bx1.*Cx2) ./ Dx1; b2xy(i, j) = -(Bx1.*Cx3) ./ Dx1;
elseif (r
=R1) ery = ((r-R1) / (r)); ethetay = ((r) / (r-R1)); sigma_py = ethetay*losstangent*((tan((omega*dt)/2))/(dt/2)); gamma_py = (2*losstangent*ery*sin((omega*dt)/2))/((1ery)*dt*cos((omega*dt)/2)); omega_py = sqrt((2*sin((omega*dt)/2)*(-2*(ery1)*sin((omega*dt)/2)+...
losstangent*ery*gamma_py*dt*cos((omega*dt)/2)))/((dt^2)*((cos((omega* dt)/2))^2))); siny = (j - yc + 0.5)/r; cosy = (i - xc)/r; ay = ((cosy.^2) + (ethetay*(siny.^2))) ./ (dt^2); by = ((gamma_py*(cosy.^2)) + ((sigma_py + (ethetay*gamma_py)).*(siny.^2))) ./ (2*dt); cy = (((cosy.^2).*((omega_py)^2)) + (sigma_py*gamma_py.*(siny.^2))) ./ (4); wy = ((1-ethetay).*siny.*cosy) ./ (dt^2); fy = ((gamma_py-sigma_py-(ethetay*gamma_py)).*siny.*cosy) ./ (2*dt); vy = ((((omega_py)^2)-(sigma_py*gamma_py)).*siny.*cosy) ./ (4); ky = 1 / (dt^2); ly = gamma_py / (2*dt); ty = ((siny.^2) + (ethetay*(cosy.^2))) ./ (dt^2); qy = ((gamma_py*(siny.^2)) + ((sigma_py + (ethetay*gamma_py)).*(cosy.^2))) ./ (2*dt); py = (((siny.^2).*((omega_py)^2)) + (sigma_py*gamma_py.*(cosy.^2))) ./ (4);
Ay1 By1 Cy1 Dy1
= = = =
ay+by+cy; Ay2 = (2.*ay)-(2.*cy); Ay3 = ay-by+cy; wy+fy+vy; By2 = (2.*wy)-(2.*vy); By3 = wy-fy+vy; ky+ly; Cy2 = (2*ky); Cy3 = ky-ly; ty+qy+py; Dy2 = (2.*ty)-(2.*py); Dy3 = ty-qy+py;
a0y(i, j) = Dy1 - (((By1).^2) ./ Ay1); a1y(i, j) = Dy2 - ((By1.*By2) ./ Ay1); a2y(i, j) = ((By1.*By3) ./ Ay1) - (Dy3);
a1yx(i, j) = By2 - ((By1.*Ay2) ./ Ay1); a2yx(i, j) = ((By1.*Ay3) ./ Ay1) - (By3);
b0y(i, j) = Cy1; b1y(i, j) = -Cy2; b2y(i, j) = Cy3; b0yx(i, j) = -(By1.*Cy1) ./ Ay1; b1yx(i, j) = (By1.*Cy2) ./ Ay1; b2yx(i, j) = -(By1.*Cy3) ./ Ay1;
elseif (r
i
r = sqrt((i - xc).^2 + (j - yc).^2); if (r<=R2) && (r>=R1) muz =
(((r-R1) / (r)) * (((R2) / (R2-R1))^2));
if (muz < 1) gamma_pz = (2*losstangent*muz*sin((omega*dt)/2))/((1muz)*dt*cos((omega*dt)/2)); omega_pz = sqrt((2*sin((omega*dt)/2)*(-2*(muz1)*sin((omega*dt)/2)+... losstangent*muz*gamma_pz*dt*cos((omega*dt)/2)))/((dt^2)*((cos((omega* dt)/2))^2))); a0z(i, a1z(i, a2z(i, b0z(i, b1z(i, b2z(i,
j) j) j) j) j) j)
= = = = = =
1/dt^2 + gamma_pz/(2*dt) + ((omega_pz^2))/4; -2/dt^2 + ((omega_pz^2))/2; 1/dt^2 - gamma_pz/(2*dt) + ((omega_pz^2))/4; 1/dt^2 + gamma_pz/(2*dt); -2/dt^2; 1/dt^2 - gamma_pz/(2*dt);
elseif (muz >= 1) sigma_pz = mu_0*muz*losstangent*((tan((omega*dt)/2))/(dt/2)); daBzx(i, j) = (1-((sigma_pz*dt)/(2*mu_0*muz))) / (1+((sigma_pz*dt)/(2*mu_0*muz))); dbBzx(i, j) = (dt/(mu_0*muz)/dx) / (1+((sigma_pz*dt)/(2*mu_0*muz))); daBzy(i, j) = (1-((sigma_pz*dt)/(2*mu_0*muz))) / (1+((sigma_pz*dt)/(2*mu_0*muz))); dbBzy(i, j) = (dt/(mu_0*muz)/dx) / (1+((sigma_pz*dt)/(2*mu_0*muz))); daBz(i, j) = (1-((sigma_pz*dt)/(2*mu_0*muz))) / (1+((sigma_pz*dt)/(2*mu_0*muz))); dbBz(i, j) = (dt/(mu_0*muz)/dx) / (1+((sigma_pz*dt)/(2*mu_0*muz)));
end end end end
%******************************************************************* % Konstanta Berenger %******************************************************************** *** sigmax = -3.0*eps_0*c_0*log(1.0e-5)/(2.0*dx*npmls); rhomax = sigmax*(aimp^2); for m=1:npmls
sig(m) = sigmax*((m-0.5)/(npmls+0.5))^2; rho(m) = rhomax*(m/(npmls+0.5))^2; end %******************************************************************* % Koefisien Berenger %******************************************************************* for m=1:npmls re = sig(m)*dt/eps_0; rm = rho(m)*dt/mu_0; ca(m) = exp(-re); cb(m) = -(exp(-re)-1.0)/sig(m)/dx; da(m) = exp(-rm); db(m) = -(exp(-rm)-1.0)/rho(m)/dx; end %******************************************************************* % Inisialisasi matriks Berenger %******************************************************************* %<<<<<<<<<<<<<<<<<<<<<<<<<<<<< Medan Hz >>>>>>>>>>>>>>>>>>>>>>>>>>> for j=1:je % Depan & belakang for i=1:npmls m = npmls+1-i; daBzx(i,j) = da(m); dbBzx(i,j) = db(m); end for i=ip+1:ie m = i-ip; daBzx(i,j) = da(m); dbBzx(i,j) = db(m); end end for i=1:ie % kiri dan kanan for j=1:npmls m = npmls+1-j; daBzy(i,j) = da(m); dbBzy(i,j) = db(m); end for j=jp+1:je m = j-jp; daBzy(i,j) = da(m); dbBzy(i,j) = db(m); end end %<<<<<<<<<<<<<<<<<<<<<<<<<<<<< Medan Ex>>>>>>>>>>>>>>>>>>>>>>>>>>>>> for i=1:ie % Kiri dan kanan for j=2:npmls+1 m = npmls+2-j; caDx(i,j) = ca(m); cbDx(i,j) = cb(m); end for j=jp+1:je m = j-jp; caDx(i,j) = ca(m); cbDx(i,j) = cb(m); end end %<<<<<<<<<<<<<<<<<<<<<<<<<<<<< Medan Ey>>>>>>>>>>>>>>>>>>>>>>>>>>>>
i
for j=1:je for i=2:npmls+1 m = npmls+2-i; caDy(i,j) = ca(m); cbDy(i,j) = cb(m); end for i=ip+1:ie m = i-ip; caDy(i,j) = ca(m); cbDy(i,j) = cb(m); end end
% Depan dan belakang
% Algoritma 1D FDTD %******************************************************************* for j=1:npmls m = npmls+1-j; daHzx_inc(j) dbHzx_inc(j) end for j=jp+1:je m = j-jp; daHzx_inc(j) dbHzx_inc(j) end
= da(m); = db(m);
= da(m); = db(m);
for j=2:npmls+1 m = npmls+2-j; caEx_inc(j) = ca(m); cbEx_inc(j) = cb(m); end for j=jp+1:je m = j-jp; caEx_inc(j) = ca(m); cbEx_inc(j) = cb(m); end %----------------------------------------------for i=1:npmls m = npmls+1-i; daHzy_inc(i) = da(m); dbHzy_inc(i) = db(m); end for i=ip+1:ie m = i-ip; daHzy_inc(i) = da(m); dbHzy_inc(i) = db(m); end for i=2:npmls+1 m = npmls+2-i; caEy_inc(i) = ca(m); cbEy_inc(i) = cb(m); end for i=ip+1:ie m = i-ip; caEy_inc(i) = ca(m);
cbEy_inc(i) = cb(m); end
%******************************************************************* % inisialisasi Movie %******************************************************************* figure('position',[200 200 700 520]); set(gcf,'color','white'); rect = get(gcf,'position'); rect(1:2) = [0 0]; Dx_h1 = Dx; Dy_h1 = Dy; Bz_h1 = Bz; Ex_h1 = Ex; Ey_h1 = Ey; Hz_h1 = Hz; clear j; %******************************************************************* % Pengulangan FDTD %******************************************************************* % sy = 500; % Perambatan persamaan Gauss arah y % sx = 500; % Perambatan persamaan Gauss arah x % n = 1; ne = 1; err = 0; n = 1; ne = 1; err = 100; % Dengan (nthreshold) while (n<ST*N) || (nthreshold)
%******************************************************************* % Distribusi Medan %******************************************************************* Dx_h2 = Dx_h1; Dx_h1 = Dx; Dy_h2 = Dy_h1; Dy_h1 = Dy; Bz_h2 = Bz_h1; Bz_h1 = Bz; Ex_h2 = Ex_h1; Ex_h1 = Ex; Ey_h2 = Ey_h1; Ey_h1 = Ey; Hz_h2 = Hz_h1; Hz_h1 = Hz;
%sumber(n) = exp(-(((n-delay)*dt)/(0.8*tau))^2) .*sin(2*pi*freq*n*dt);
%******************************************************************* % Ex_inc datang (1-D FDTD) %******************************************************************* Ex_inc(2:jb) = caEx_inc(2:jb).*Ex_inc(2:jb) + cbEx_inc(2:jb).*(Hzx_inc(1:je)-Hzx_inc(2:jb));
%******************************************************************* % Medan Dx dan Dy field %******************************************************************* Dx( : ,2:je) = caDx( : ,2:je).*Dx( : ,2:je) + cbDx( : ,2:je).*(Hz( : ,2:je) - Hz( : ,1:je-1)); Dy(2:ie, : ) = caDy(2:ie, : ).*Dy(2:ie, : ) + cbDy(2:ie, : ).*(Hz(1:ie-1,: ) - Hz(2:ie, : ));
i
%******************************************************************* % Ruang penghubung D dan E %******************************************************************* Dysynch(1:ie, 2:je) = 0.25*(Dy(1:ie, 2:je)+Dy(2:ib, 2:je)+Dy(1:ie, 1:je-1)+Dy(2:ib, 1:je-1)); Dy_h1synch(1:ie, 2:je) = 0.25*(Dy_h1(1:ie, 2:je)+Dy_h1(2:ib, 2:je)+Dy_h1(1:ie, 1:je-1)+Dy_h1(2:ib, 1:je-1)); Dy_h2synch(1:ie, 2:je) = 0.25*(Dy_h2(1:ie, 2:je)+Dy_h2(2:ib, 2:je)+Dy_h2(1:ie, 1:je-1)+Dy_h2(2:ib, 1:je-1)); Ey_h1synch(1:ie, 2:je) = 0.25*(Ey_h1(1:ie, 2:je)+Ey_h1(2:ib, 2:je)+Ey_h1(1:ie, 1:je-1)+Ey_h1(2:ib, 1:je-1)); Ey_h2synch(1:ie, 2:je) = 0.25*(Ey_h2(1:ie, 2:je)+Ey_h2(2:ib, 2:je)+Ey_h2(1:ie, 1:je-1)+Ey_h2(2:ib, 1:je-1)); Dxsynch(2:ie, 2:je) = 0.25*(Dx(2:ie, 2:je)+Dx(2:ie, 3:jb)+Dx(1:ie-1, 2:je)+Dx(1:ie-1, 3:jb)); Dx_h1synch(2:ie, 2:je) = 0.25*(Dx_h1(2:ie, 2:je)+Dx_h1(2:ie, 3:jb)+Dx_h1(1:ie-1, 2:je)+Dx_h1(1:ie-1, 3:jb)); Dx_h2synch(2:ie, 2:je) = 0.25*(Dx_h2(2:ie, 2:je)+Dx_h2(2:ie, 3:jb)+Dx_h2(1:ie-1, 2:je)+Dx_h2(1:ie-1, 3:jb)); Ex_h1synch(2:ie, 2:je) = 0.25*(Ex_h1(2:ie, 2:je)+Ex_h1(2:ie, 3:jb)+Ex_h1(1:ie-1, 2:je)+Ex_h1(1:ie-1, 3:jb)); Ex_h2synch(2:ie, 2:je) = 0.25*(Ex_h2(2:ie, 2:je)+Ex_h2(2:ie, 3:jb)+Ex_h2(1:ie-1, 2:je)+Ex_h2(1:ie-1, 3:jb));
%******************************************************************* % Medan Ex dan Ey %******************************************************************* Ex(1:ie, 2:je) = ((b0x(1:ie,2:je).*Dx(1:ie,2:je) + b0xy(1:ie,2:je).*Dysynch(1:ie,2:je) + ... b1x(1:ie,2:je).*Dx_h1(1:ie,2:je) + b1xy(1:ie,2:je).*Dy_h1synch(1:ie,2:je) + ... a1x(1:ie,2:je).*Ex_h1(1:ie,2:je) + a1xy(1:ie,2:je).*Ey_h1synch(1:ie,2:je) + ... b2x(1:ie,2:je).*Dx_h2(1:ie,2:je) + b2xy(1:ie,2:je).*Dy_h2synch(1:ie,2:je) + ... a2x(1:ie,2:je).*Ex_h2(1:ie,2:je) + a2xy(1:ie,2:je).*Ey_h2synch(1:ie,2:je)) ... ./ (a0x(1:ie,2:je))); Ey(2:ie, 1:je) = ((b0y(2:ie, 1:je).*Dy(2:ie, 1:je) + 1:je).*Dxsynch(2:ie, 1:je) + ... b1y(2:ie, 1:je).*Dy_h1(2:ie, 1:je) 1:je).*Dx_h1synch(2:ie, 1:je) + ... a1y(2:ie, 1:je).*Ey_h1(2:ie, 1:je) 1:je).*Ex_h1synch(2:ie, 1:je) + ... b2y(2:ie, 1:je).*Dy_h2(2:ie, 1:je) 1:je).*Dx_h2synch(2:ie, 1:je) + ... a2y(2:ie, 1:je).*Ey_h2(2:ie, 1:je) 1:je).*Ex_h2synch(2:ie, 1:je)) ... ./ (a0y(2:ie, 1:je)));
b0yx(2:ie, + b1yx(2:ie, + a1yx(2:ie, + b2yx(2:ie, + a2yx(2:ie,
%******************************************************************* % Nilai Ex datang (1-D & 2-D)
%******************************************************************* Ex(it+1:ie-it, jt+1) = Ex(it+1:ie-it, jt+1) + (dt/eps_0/dx).*Hzx_inc( jt+1); %.*exp(-((it+1:ie-it)-ic).^2./(2*sy^2))'; % Left Ex(it+1:ie-it,je-jt+1) = Ex(it+1:ie-it,je-jt+1) (dt/eps_0/dx).*Hzx_inc(je-jt+1); %.*exp(-((it+1:ie-it)-ic).^2./(2*sy^2))'; % Right
%******************************************************************* % Nilai Ey datang(1-D & 2-D) %******************************************************************* Ey( it+1,jt+1:je-jt) = Ey( it+1,jt+1:je-jt) (dt/eps_0/dx).*Hzx_inc(jt+1:je-jt); %.*exp(-(( it+1)-ic).^2./(2*sy^2))'; % Bottom Ey(ie-it+1,jt+1:je-jt) = Ey(ie-it+1,jt+1:je-jt) + (dt/eps_0/dx).*Hzx_inc(jt+1:je-jt); %.*exp(-((ie-it+1)-ic).^2./(2*sy^2))'; % Top
%******************************************************************* % Penghitungan Hzx_inc (1-D FDTD) %******************************************************************* Hzx_inc(1:je) = daHzx_inc(1:je).*Hzx_inc(1:je) + dbHzx_inc(1:je).*(Ex_inc(1:je)-Ex_inc(2:jb));
%******************************************************************* % ..... Eksitasi Sumber Datang(1-D FDTD)..... %******************************************************************* Hzx_inc(js) = Hzx_inc(js) + 1*source(n);
%******************************************************************* % medan Bzx dan Bzy (Bz) %******************************************************************* Bzx(1:ie,: ) = daBzx(1:ie,: ).*Bzx(1:ie,: ) + dbBzx(1:ie,: ).*(Ey(1:ie,: ) - Ey(2:ib, :)); Bzy( :, 1:je) = daBzy( :, 1:je).*Bzy( :, 1:je) + dbBzy( :, 1:je).*(Ex( :, 2:jb) - Ex( :, 1:je)); Bz = Bzx + Bzy;
%******************************************************************* % Medan Hz %******************************************************************* Hz = (b0z.*Bz + b1z.*Bz_h1 + b2z.*Bz_h2 - (a1z.*Hz_h1 + a2z.*Hz_h2)) ./ a0z;
%******************************************************************* % Nilai Hz datang (1-D & 2-D) %******************************************************************* Hz(it+1:ie-it, jt+1) = Hz(it+1:ie-it, jt+1) (dt/mu_0/dx).*Ex_inc( jt+1); % Left Hz(it+1:ie-it,je-jt+1) = Hz(it+1:ie-it,je-jt+1) + (dt/mu_0/dx).*Ex_inc(je-jt+1); % Right
i
%HzsourceTF(1:length(it:ie-it),n)=Hz(it:ie-it,jt+1); %HztransTF(1:length(it:ie-it),n)=Hz(it:ie-it,je-jt-1); %HzsourceSF(1:length(it:ie-it),n)=Hz(it:ie-it,js); %HztransSF(1:length(it:ie-it),n)=Hz(it:ie-it,je-js);
%******************************************************************* % Visualisasi Medan %******************************************************************* if mod(n,10)==0 timestep = num2str(n); pcolor((1:je)./M,(1:ie)./M,real(Hz)); shading interp; axis image; colorbar; axis([1/M je/M 1/M ie/M]); %pcolor((1:je)./M,(1:ie)./M,real(Hz)/max(max(real(Hz)))); shading interp; axis image; colorbar; axis([1/M je/M 1/M ie/M]); line((R1*sin((0:360)*pi/180)+yc)./M,(R1*cos((0:360)*pi/180)+xc)./M,'c olor','k'); line((R2*sin((0:360)*pi/180)+yc)./M,(R2*cos((0:360)*pi/180)+xc)./M,'c olor','k'); %xlabel('y/\lambda'); ylabel('x/\lambda'); title(['amplitude of H_z']); xlabel('y/\lambda'); ylabel('x/\lambda'); title(['H_z at timestep ', timestep]); getframe(gcf,rect); end %******************************************************************* % Penghitungan fungsi eror untuk konvergensi pulsa %****************************************************************** if mod(n,N)==0 ne = n/N; if (20*log10(sum(sum(abs(Hz(it:ie-it,jt:je-jt)))))) > (20*log10(sum(sum(abs(Hz_p(it:ie-it,jt:je-jt)))))) err(ne) = 0; else err(ne) = (20*log10(sum(sum(abs(Hz(it:ie-it,jt:je-jt)))))) (20*log10(sum(sum(abs(Hz_p(it:ie-it,jt:je-jt)))))); end Hz_p = Hz; end
%******************************************************************* %Penghitungan fungsi eror untuk konvergensi pulsa %******************************************************************* if mod(n,N)==0 ne = n/N; if max(max(abs(Hz_p(npmls+10:ip-10,npmls+10:jp-10))))==0 err(ne) = 100; else
err(ne) = max(max(abs(abs(Hz(npmls+10:ip-10,npmls+10:jp10))-abs(Hz_p(npmls+10:ip-10,npmls+10:jp-10))))) / max(max(abs(Hz_p(npmls+10:ip-10,npmls+10:jp-10)))) * 100; end Hz_p = Hz; end %******************************************************************* % Menampilkan %******************************************************************* if mod(n,step)==0 %save results_2D.mat dx dt freq nmax ie je is js err N M Ex Ey Hz; c = round(clock); disp(['Current time step: ',int2str(n-mod(n,step)),', ',int2str(n/nmax*100),'% completed at ',num2str(c(4)),':',num2str(c(5)),':',num2str(c(6))]); end n = n + 1; end %******************************************************************* % Menyimpan %****************************************************************** save results_2DTEidealcloakdielcore3.mat dx dt freq n nmax ie je xc yc R1 R2 is js err N M Ex Ey Hz;
i
APENDIX F MENGUBAH PERSAMAAN-PERSAMAAN KE PROGRAM
Persamaan (4.40) : ( (
) )
(
( )
)
(
)
Diubah dalam program menjadi :
(
)
(
)
a = 4/4*e0*einf + e0*wpesq*dt^2 + e0*einf*ge*2*dt* a = 4/e0*4*einf+e0*dt^2*wpesq+e0*2*dt*einf*ge a = (1/dt^2)*(4*dt^2)/(e0*(4*einf+dt^2*wpesq+2*dt*einf*ge) a = (1/dt^2)*a0 a0 = (4*dt^2)/(e0*(4*einf+dt^2*wpesq+2*dt*einf*ge))
(
)
( (
) )
(
)
b = (ge*2*dt)/(4*e0*einf+e0*wpesq*dt^2+e0*einf*ge*2*dt) b = (ge*2*dt)/(e0*4*einf+e0*dt^2*wpesq+e0*2*dt*einf*ge) b = (1/(2*dt))*ge*(4*dt^2)/(e0*(4*einf+dt^2*wpesq+2*dt*einf*ge)) b = (1/(2*dt))*ge.*a0
(
)
(
)
c = 4*e0*einf/(4*e0*einf+e0*wpesq*dt^2+e0*einf*ge*2*dt) c = (e0/dt^2)*einf*(4*dt^2)/e0*(4*einf+dt^2*wpesq+2*dt*einf*ge) c = (e0/dt^2)*einf.*a0;
( ) ( )
(
)
d = (-e0*wpesq*dt^2)/(4*e0*einf+e0*wpesq*dt^2+e0*einf*ge*2*dt) d = (-1*e0/4)*wpesq*(4*dt^2)/e0*(4*einf+dt^2*wpesq+2*dt*einf*ge) d = (-1*e0/4)*wpesq*a0
( )
(
) (
)
e = e0*einf*ge2*dt/(4*e0*einf+e0*wpesq*dt^2+e0*einf*ge*2*dt) e = (1/(2*dt))*e0*einf*ge *(4*dt^2)/e0*(4*einf+dt^2*wpesq+2*dt*einf*ge) e = (1/(2*dt))*e0*einf*ge*a0
Selanjutnya untuk medan magnet persamaanya ( (
)
) (
( )
Diubah dalam program menjadi : i
(
) )
(
)
(
)
am = 4/4*u0*uinf + u0*wpesq*dt^2 + u0*uinf*gm*2*dt* am = 4/u0*4*uinf+u0*dt^2*wpesq+u0*2*dt*uinf*gm am = (1/dt^2)*(4*dt^2)/(u0*(4*uinf+dt^2*wpesq+2*dt*uinf*gm) am = (1/dt^2)*am0 am0 = (4*dt^2)/(u0*(4*uinf+dt^2*wpesq+2*dt*uinf*gm))
(
)
( (
) )
(
)
bm = (gm*2*dt)/(4*u0*uinf+u0*wpesq*dt^2+u0*uinf*gm*2*dt) bm = (gm*2*dt)/(u0*4*uinf+u0*dt^2*wpesq+u0*2*dt*uinf*gm) bm = (1/(2*dt))*gm*(4*dt^2)/(u0*(4*uinf+dt^2*wpesq+2*dt*uinf*gm)) bm = (1/(2*dt))*gm.*am0
(
)
(
)
cm = 4*u0*uinf/(4*u0*uinf+u0*wpesq*dt^2+u0*uinf*gm*2*dt) cm = (u0/dt^2)*uinf*(4*dt^2)/u0*(4*uinf+dt^2*wpesq+2*dt*uinf*gm) cm = (u0/dt^2)*uinf.*am0
( )
(
) (
)
dm = (-u0*wpesq*dt^2)/(4*u0*uinf+u0*wpesq*dt^2+u0*uinf*gm*2*dt) dm = (-1*u0/4)*wpesq*(4*dt^2)/u0*(4*uinf+dt^2*wpesq+2*dt*uinf*gm) dm = (-1*u0/4)*wpesq*am0
( (
)
)
(
)
em = u0*uinf*gm2*dt/(4*u0*uinf+u0*wpesq*dt^2+u0*uinf*gm*2*dt) em = (1/(2*dt))*u0*uinf*gm *(4*dt^2)/u0*(4*uinf+dt^2*wpesq+2*dt*uinf*gm) em = (1/(2*dt))*u0*uinf*gm*am0
persamaan (4.35) : ⃗
[
⃗
]
[
(⃗
]
⃗
)
Hy(1:SIZE-1,n2) = Hy(1:SIZE-1,n1) + ( ( Ex(1:SIZE-1,n1) - Ex(2:SIZE,n1) ) * dt/(u0*dz) )
persamaan (4.36) : ⃗
⃗
(⃗
[
]
⃗
[
])
Ex(2:SIZE,n2) = Ex(2:SIZE, n1) + ( dt/(e0*dz)*(Hy(1:SIZE-1, n2) - Hy(2:SIZE, n2)) )
persamaan (4.39), (4.41),(4.42) (4.43) dan (4.44) : ⃗
(⃗
⃗ ( ⃗
⃗ ⃗
(⃗
) (⃗
)
⃗
( ⃗
)
⃗
)
)
Ex(:,n2) = a.*(Dx(:,n2)-2*Dx(:,n1)+Dx(:,3))+b.*(Dx(:,n2)Dx(:,3))+c.*(2*Ex(:,n1)-Ex(:,3))+d.*(2*Ex(:,n1) +Ex(:,3)) +e.*(Ex(:,3))
⃗
(⃗
⃗
⃗
)
(⃗
( ⃗
⃗
)
(⃗
⃗
( ⃗
)
)
Hy(:,n2) = am.*(By(:,n2)-2*By(:,n1)+By(:,3))+bm.*(By(:,n2)By(:,3))+cm.*(2*Hy(:,n1)-Hy(:,3))+dm.*(2*Hy(:,n1) +Hy(:,3))+em.*(Hy(:,3))
⃗
( )
⃗ ( )
(⃗ ( )
⃗ (
))
By(1:SIZE-1,n2) = By(1:SIZE-1,n1) + ( ( Ex(1:SIZE-1,n1) – Ex(2:SIZE,n1) ) * dt/(dz) )
⃗
( )
⃗ ( )
(⃗ ( )
⃗ (
))
Dx(2:SIZE,n2) = Dx(2:SIZE, n1) + ( dt/(dz)*(Hy(1:SIZE-1, n2) –
i
⃗
)
Hy(2:SIZE, n2)) )
(
|
)
̃ ( ̃ (
) | )
nFDTD = (1/(1i*k0*(z1-z2))) .*log(EXZ2(1:NFFT/2+1)./EXZ1(1:NFFT/2+1))
( )
( )
erx = ((r-R1) / (r));
( ) ethetax =
( )
((r) / (r-R1));
⃗
[
⃗
⃗
⃗
⃗
⃗
⃗
⃗
⃗
]
⃗ ⃗
[
⃗
⃗
⃗
[
⃗
⃗
⃗
⃗
⃗
⃗ ⃗
⃗
⃗
⃗
⃗
⃗
(
)
]
⃗
⃗
⃗
⃗
⃗
⃗
]
⃗ ⃗
[
⃗
⃗ (
Ax1 = ax+bx+cx
⃗ )
⃗ ⃗
⃗ ]
(
)
Ax2 = (2.*ax)-(2.*cx) (
)
(
)
Ax3 = ax-bx+cx ax = ((cosx.^2) + (ethetax.*(sinx.^2))) ./ (dt^2); bx = ((gamma_px.*(cosx.^2)) + ((sigma_px +(ethetax*gamma_px)) .*(sinx.^2))) ./ (2*dt); cx = (((cosx.^2).*((omega_px)^2)) + (sigma_px*gamma_px .*(sinx.^2))) ./ (4); (
)
(
)
(
)
(
)
Bx1 = wx+fx+vx (
)
(
)
Bx2 = (2.*wx)-(2.*vx) (
)
(
)
Bx3 = wx-fx+vx; wx = ((1-ethetax).*sinx.*cosx) ./ (dt^2) fx = ((gamma_px-sigma_px-(ethetax*gamma_px)).*sinx.*cosx)./(2*dt) vx = (((omega_px)^2)-(sigma_px*gamma_px)).*sinx.*cosx) ./ (4)
Cx1 = kx+lx
Cx2 = (2.*kx)
Cx3 = kx-lx
i
kx = 1 / etheta0*(dt^2); lx = gamma_px / etheta0* (2*dt); (
)
(
)
(
)
Dx1 = tx+qx+px (
)
Dx2 = (2.*tx)-(2.*px) (
)
Dx3 = tx-qx+px tx = ((sinx.^2) + (ethetax.*(cosx.^2))) ./ (dt^2) qx = ((gamma_px.*(sinx.^2)) + ((sigma_px +(ethetax*gamma_px)) .*(cosx.^2)))./(2*dt) px = (((sinx.^2).*((omega_px)^2)) + (sigma_px*gamma_px .*(cosx.^2))) ./(4)
a0x(i, j) = Ax1 - (((Bx1).^2) ./ Dx1); a1x(i, j) = Ax2 - ((Bx1.*Bx2) ./ Dx1); a2x(i, j) = ((Bx1.*Bx3) ./ Dx1) - (Ax3);
a1xy(i, j) = Bx2 - ((Bx1.*Dx2) ./ Dx1); a2xy(i, j) = ((Bx1.*Dx3) ./ Dx1) - (Bx3);
b0x(i, j) = Cx1; b1x(i, j) = -Cx2; b2x(i, j) = Cx3;
b0xy(i, j) = -(Bx1.*Cx1) ./ Dx1; b1xy(i, j) = (Bx1.*Cx2) ./ Dx1; b2xy(i, j) = -(Bx1.*Cx3) ./ Dx1;
untuk bagian y : (
)
(
)
Ay1 = ay+by+cy (
)
Ay2 = (2.*ay)-(2.*cy) (
)
(
)
Ay3 = ay-by+cy ay = ((cosy.^2) + (ethetay.*(siny.^2))) ./ (dt^2); by = ((gamma_py.*(cosy.^2)) + ((sigma_py +(ethetay*gamma_py)) .*(siny.^2))) ./ (2*dt); cy = (((cosy.^2).*((omega_py)^2)) + (sigma_py*gamma_py .*(siny.^2))) ./ (4); (
)
(
)
(
)
(
)
By1 = wy+fy+vy (
)
(
)
By2 = (2.*wy)-(2.*vy) (
)
(
)
i
By3 = wy-fy+vy; wy = ((1-ethetay).*siny.*cosy) ./ (dt^2) fy = ((gamma_py-sigma_py-(ethetay*gamma_py)).*siny.*cosy)./(2*dt) vy = (((omega_py)^2)-(sigma_py*gamma_py)).*siny.*cosy) ./ (4)
Cy1 = ky+ly
Cy2 = (2.*ky)
Cy3 = ky-ly ky = 1 / etheta0*(dt^2); ly = gamma_py / etheta0* (2*dt); (
)
(
)
(
)
Dy1 = ty+qy+py (
)
Dy2 = (2.*ty)-(2.*py) (
)
Dy3 = ty-qy+py ty = ((siny.^2) + (ethetay.*(cosy.^2))) ./ (dt^2) qy = ((gamma_py.*(siny.^2)) + ((sigma_py +(ethetay*gamma_py)) .*(cosy.^2)))./(2*dt) py = (((siny.^2).*((omega_py)^2)) + (sigma_py*gamma_py .*(cosy.^2))) ./(4)
a0y(i, j) = Ay1 - (((By1).^2) ./ Dy1); a1y(i, j) = Ay2 - ((By1.*By2) ./ Dy1); a2y(i, j) = ((By1.*By3) ./ Dy1) - (Ay3);
a1yx(i, j) = By2 - ((By1.*Dy2) ./ Dy1); a2yx(i, j) = ((By1.*Dy3) ./ Dy1) - (By3);
b0y(i, j) = Cy1; b1y(i, j) = -Cy2; b2y(i, j) = Cy3;
b0yx(i, j) = -(By1.*Cy1) ./ Dy1; b1yx(i, j) = (By1.*Cy2) ./ Dy1; b2yx(i, j) = -(By1.*Cy3) ./ Dy1;
i
