41
LAMPIRAN
42
43
Lampiran 1. Diagram alir penelitian
Penelusuran Literatur
Sudah Siap
Penguasaan Software
Pembuatan dan Pengujian Program
Analisis Output
Penyusunan Laporan
Penentuan Parameter
44
45
Lampiran 2. Program coding untuk mensimulasikan proses perambatan gelombang elektromagnetik di dalam struktur sensor dengan menggunakan metode FDTD % % % % %
FDTD 2D dengan PML % optical biosensor kuasi kristal fotonik 1D dengan DUA ROD DEFEK Bahan penelitian Tesis Sumber berupa gelombang sinusoidal Dani
clear all; close all; clc; %----------------------------------------------------------------------% Kondisi Inisial. %----------------------------------------------------------------------Imax = 200; % Jumlah grid. Jmax = 400; Ez = zeros(Imax,Jmax); Dz = zeros(Imax,Jmax); Dz_hat = zeros(Imax,Jmax); Hy = zeros(Imax,Jmax); sumbu-y Hx = zeros(Imax,Jmax); sumbu-x fi1 fi2 fi3 gi2 gi3
= = = = =
zeros(Imax,1); ones(Imax,1); ones(Imax,1); ones(Imax,1); ones(Imax,1);
fj1 fj2 fj3 gj2 gj3
= = = = =
zeros(Jmax,1); ones(Jmax,1); ones(Jmax,1); ones(Jmax,1); ones(Jmax,1);
% % % %
vektor vektor vektor vektor
medan listrik flux listrik parameter medan magnet arah
% vektor medan magnet arah
iHx = zeros(Imax,Jmax); iHy = zeros(Imax,Jmax); c0 = 3e8; lambda=560e-9; freq=c0/lambda; excitation omega=2.0*pi*freq; k=2*pi/lambda; time=0; t0=50; spread=15;
%center wavelength of source excitation %center frequency of source
% inisial kondisi waktu
46
%----------------------------------------------------------------------% vektor konstanta dielektrik input. %----------------------------------------------------------------------e0 = 8.85e-12; er = ones(Imax,Jmax); %----------------------------------------------------------------------% konduktivitas input. konduktivity ohmik: J = sigma*E. %----------------------------------------------------------------------sigma = zeros(Imax,Jmax); %----------------------------------------------------------------------% parameter PML. %----------------------------------------------------------------------npml = 8; for i= 0:npml xnum = npml - i; % Dz xxn = xnum/npml; xn = .333*xxn^3; gi2(i+1) = 1/(1+xn); gi2(Imax-i) = 1/(1+xn); gi3(i+1) = (1-xn)/(1+xn); gi3(Imax-i) = (1-xn)/(1+xn); % for H_x dan H_y xxn = (xnum - .5)/npml; xn = .333*xxn^3; fi1(i+1) = xn; fi1(Imax-i) = xn; fi2(i+1) = 1/(1+xn); fi2(Imax-i) = 1/(1+xn); fi3(i+1) = (1-xn)/(1+xn); fi3(Imax-i) = (1-xn)/(1+xn); end for j= 0:npml xnum = npml-j; % Dz xxn = xnum/npml; xn = .333*xxn^3; gj2(j+1) = 1/(1+xn); gj2(Jmax-j) = 1/(1+xn); gj3(j+1) = (1-xn)/(1+xn); gj3(Jmax-j) = (1-xn)/(1+xn); % for H_x dan H_y xxn = (xnum - .5)/npml; xn = .333*xxn^3; fj1(j+1) = xn; fj1(Jmax-j) = xn; fj2(j+1) = 1/(1+xn); fj2(Jmax-j) = 1/(1+xn); fj3(j+1) = (1-xn)/(1+xn);
47
fj3(Jmax-j) = (1-xn)/(1+xn); end %----------------------------------------------------------------------% penentuan ukuran ddx dan dt. %----------------------------------------------------------------------ddx = 5.0e-8; % space increment dt = ddx/6e8; % time step %--------------------------------------------------------------------% Variasi indeks bias pada struktur sensor %--------------------------------------------------------------------er_A=(3.48)^2; % Si er_B=(1.44)^2; % SiO2 er_C= (1.7)^2; % Al2O3 for i=1:80 for j=30:380 er(i,j)=er_B; sigma(i,j)=0.0; end end for i=81:120 for j=30:380 sigma(i,j)=0.0; A1=zeros(i,j); A2=zeros(i,j); A3=zeros(i,j); A4=zeros(i,j); A5=zeros(i,j); A6=zeros(i,j); A7=zeros(i,j); A8=zeros(i,j); A9=zeros(i,j); A10=zeros(i,j); A11=zeros(i,j); radius = 10; radius2 = 12; xc=100; yc1=60; yc2=90; yc3=120; yc4=150; yc5=180; yc6=210; yc7=240; yc8=270; yc9=300; yc10=330; yc11=360; A1(i,j)=sqrt((i-xc)^2 + (j-yc1)^2); A2(i,j)=sqrt((i-xc)^2 + (j-yc2)^2); A3(i,j)=sqrt((i-xc)^2 + (j-yc3)^2); A4(i,j)=sqrt((i-xc)^2 + (j-yc4)^2); A5(i,j)=sqrt((i-xc)^2 + (j-yc5)^2); A6(i,j)=sqrt((i-xc)^2 + (j-yc6)^2); A7(i,j)=sqrt((i-xc)^2 + (j-yc7)^2); A8(i,j)=sqrt((i-xc)^2 + (j-yc8)^2); A9(i,j)=sqrt((i-xc)^2 + (j-yc9)^2); A10(i,j)=sqrt((i-xc)^2 + (j-yc10)^2); A11(i,j)=sqrt((i-xc)^2 + (j-yc11)^2); if A1(i,j)<=radius er(i,j)=er_B;
48
elseif A2(i,j)<=radius er(i,j)=er_B; elseif A3(i,j)<=radius er(i,j)=er_B; elseif A4(i,j)<=radius2 er(i,j)=er_C; elseif A5(i,j)<=radius er(i,j)=er_B; elseif A6(i,j)<=radius er(i,j)=er_B; elseif A7(i,j)<=radius er(i,j)=er_B; elseif A8(i,j)<=16 er(i,j)=1.45^2; elseif A9(i,j)<=radius er(i,j)=er_B; elseif A10(i,j)<=radius er(i,j)=er_B; elseif A11(i,j)<=radius er(i,j)=er_B; else er(i,j)=er_A; end
% Defek ke-1 (rod 4)
% Defek ke-2 (rod 8)
end end for i=120:200 for j=30:380 er(i,j)=er_B; sigma(i,j)=0.0; end end %----------------------------------------------------------------------% Penentuan vektor dan konstanta. %----------------------------------------------------------------------gb = (dt*sigma)./e0; gbc = zeros(Imax,Jmax); ga = 1./(er + gb + gbc); %----------------------------------------------------------------------% variabel awal. %----------------------------------------------------------------------T = 0; NSTEPS =2500; %----------------------------------------------------------------------% MAIN FDTD LOOP. %-----------------------------------------------------------------------
49
for n = 1:NSTEPS time= time+n*dt; t(n)=time; T = T + 1; %--------------------------------------------------------------------% Penentuan Dz dari Hy and Hx. %--------------------------------------------------------------------for i = 2:Imax for j = 2:Jmax Dz_hat_temp = gi3(i)*Dz_hat(i,j) + gi2(i)*.5*(Hy(i,j)-Hy(i1,j)-Hx(i,j)+Hx(i,j-1)); Dz(i,j) = gj3(j)*Dz(i,j) + gj2(j)*(Dz_hat_temp Dz_hat(i,j)); Dz_hat(i,j) = Dz_hat_temp; end end %--------------------------------------------------------------------% Pembentukan pulsa listrik. %--------------------------------------------------------------------Ic = 150; % domain komputasi Jc = 50; rtau=160.0e-12; tau=rtau/dt; delay=3*tau; for i = 1:Imax for j = 1:Jmax pulse = 50*sin(2*pi*(freq*dt*T-(freq/3e8)*ddx)); Dz(i,10) = Dz(i,10) + pulse; end end % %
source = 50*sin(2*pi*1500*1e6*dt*T); Dz(Ic,Jc) = Dz(Ic,Jc) + source;
% %
source = 50*sin(2*pi*(1500*1e6*dt*T-(1500*1e6/3e8)*ddx)); Dz(Ic,Jc) = Dz(Ic,Jc) + source;
% %
source = exp(-.5*((t0-T)/spread)^2 ); Dz(Ic,Jc) = Dz(Ic,Jc) + source; %--------------------------------------------------------------------% Penentuan Ez dari Dz. %--------------------------------------------------------------------Ez = ga.*Dz; Ez(1,:) = 0; Ez(Imax,:) = 0;
Ez(:,1) = 0; Ez(:,Jmax) = 0;
50
%--------------------------------------------------------------------% penentuan Hx dan Hy dari Ez. %--------------------------------------------------------------------for i = 1:Imax-1 for j = 1:Jmax-1 % Hy. dE = Ez(i+1,j) - Ez(i,j); iHy(i,j) = iHy(i,j) + fj1(j)*dE; Hy(i,j) = fi3(i)*Hy(i,j) + fi2(i)*.5*dE + fi2(i)*iHy(i,j); % Hx. dE = Ez(i,j) - Ez(i,j+1); iHx(i,j) = iHx(i,j) + fi1(i)*dE; Hx(i,j) = fj3(j)*Hx(i,j) + fj2(j)*.5*dE + fj2(j)*iHx(i,j); end end %--------------------------------------------------------------------%Perhitungan input dan output %--------------------------------------------------------------------for i=81:120 Q1=abs(Ez(i,28)); Q2=abs(Ez(i,382)); end for a=1:40 p(a)=a; Z1(a)=e0*Q1; Z2(a)=e0*Q2; end Az1(n)=trapz(p,Z1); Az2(n)=trapz(p,Z2);
timestep=int2str(n); figure(1) surf(real(Ez)); caxis([-3e4 3e4]); title(['Ez at time step = ',timestep]); xlabel('x-direction (*5E-2 µm)'); ylabel('y direction (*5E-2 µm)'); zlabel('Ez Field'); view(0,90) shading interp hold on hold off; pause(0.1)
51
end Energi_In=trapz(t,Az1) Energi_Out=trapz(t,Az2)
figure(2) surf(abs(er)); title('Struktur Berdasarkan Nilai Permitivitas') xlabel('x-direction (*5E-2 µm)'); ylabel('y direction (*5E-2 µm)'); zlabel('relative permitivity') view(0,90) shading interp figure(3) title('Power input'); plot(t,Az1); xlabel('Time(second)'); ylabel('Power Input(watt)'); figure(4) plot(t,Az2); title('Power Output er=11.00') xlabel('Time(second)'); ylabel('Power Output(watt)');
52
53
Lampiran 3. Proses pengukuran nilai rapat energi rata-rata W
Incident wave h
y
2 Q( t ) = ∫ε E( t) dy 0
defek st
1 defect
t
1 W = ∫ Q( t ) dt t0
nd
x
z
2 defect
Simulasi dilakukan beberapa kali dengan memvariasikan nilai indeks bias defek kedua
Grafik hubungan rapat energi rata-rata terhadap indeks bias
54
55
PUBLIKASI
56
57
PROCEEDING – The 4th Asian Physics Symposium 2010 (APS 2010), Bandung
58
59
