%Metal heterostructure-based nanophotonic devices: finite-difference time-domain numerical simulations
%fig8
format long 
Ni=350;Nj=150; thick=25; w=5; d1=18;d2=18;pb=15;Lb=100;La=288;%pb是pml的層數,La是長，Lb是寬
c=299792458.0; pi=3.1415926; lambda0=780*10^(-9); f0=c/lambda0; omega=2*pi*f0;mu0=12.5663706144*10.^(-7); epsilon0=8.854187818*10.^(-12);
omegap2=1.37*10^(16);gamm2=0.15*10^(13); %Ag的參數
omegap1=2.27*10^(16);gamm1=0.15*10^(14);%al的參數
deltax=5*10^(-9);deltay=deltax; deltat=8*10^(-18); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%pml的設置
CA=1; CB=deltat./(epsilon0); CP=1; CQ=deltat./(mu0);
sigammax=5/(150*pi*deltax);sigammmax=sigammax*(mu0/epsilon0);  bs1=([1:pb]'/(pb)).^4;bs2=([0.5:pb-0.5]'/pb).^4;
bahalf=(1/(2*pb)^4)*sigammax;  bs3=(bs1*ones(1,Ni))';bs4=(ones(Nj,1)*bs1')';
%bs5=(bs2*ones(1,Ni))';bs6=(ones(Nj,1)*bs2')';
CAhalf=exp(-(bahalf*deltat/epsilon0));CBhalf=(1-CAhalf)./(deltay*bahalf); 
byu=bs3*sigammax; byd=rot90(byu,2); 
%ayu=bs5*sigammax;ayd=rot90(ayu,2);
bydm=(mu0/epsilon0)*byd; byum=(mu0/epsilon0)*byu;
byuu=bs4*sigammax; bydd=rot90(byuu,2);
%ayuu=bs6*sigammax; aydd=rot90(ayuu,2);
byddm=(mu0/epsilon0)*bydd;byuum=(mu0/epsilon0)*byuu;
CA3xu=exp(-(byu*deltat/epsilon0)); CB3xu=(1-CA3xu)./(deltay*byu);
CP3yu=exp(-(byum*deltat/mu0)); CQ3yu=(1-CP3yu)./(deltax*byum);%上邊界參數

CA3xd=exp(-(byd*deltat/epsilon0)); CB3xd=(1-CA3xd)./(deltay*byd);
CP3yd=exp(-(bydm*deltat/mu0)); CQ3yd=(1-CP3yd)./(deltax*bydm);%下邊界參數

CA4yu=exp(-(byuu*deltat/epsilon0)); CB4yu=(1-CA4yu)./(deltax*byuu);
CP4xu=exp(-(byuum*deltat/mu0)); CQ4xu=(1-CP4xu)./(deltay*byuum);%右邊界參數

CA4yd=exp(-(bydd*deltat/epsilon0)); CB4yd=(1-CA4yd)./(deltax*bydd);  
CP4xd=exp(-(byddm*deltat/mu0)); CQ4xd=(1-CP4xd)./(deltay*byddm);%左邊界參數
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Hz=zeros(Ni,Nj); Ex=zeros(Ni,Nj+1); Ey=zeros(Ni+1,Nj); Hzx=Hz; Hzy=Hz;Hz1=Hz;AHz=Hz; AEy=Ey;AEx=Ex;
Jx=zeros(La,Lb+1); Jy=zeros(La+1,Lb);
CA2x=Jx; CD2x=Jx; CA2y=Jy; CD2y=Jy;gammnewx=Jx; omegapnewx=Jx; gammnewy=Jy; omegapnewy=Jy;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for I=1:8
gammnewx(36*(I-1)+1:36*(I-1)+18,:)=gamm1+gammnewx(36*(I-1)+1:36*(I-1)+18,:);
gammnewx(36*(I-1)+19:36*I,:)=gamm2+gammnewx(36*(I-1)+19:36*I,:);
omegapnewx(36*(I-1)+1:36*(I-1)+18,:)=omegap1+omegapnewx(36*(I-1)+1:36*(I-1)+18,:);
omegapnewx(36*(I-1)+19:36*I,:)=omegap2+omegapnewx(36*(I-1)+19:36*I,:);

gammnewy(36*(I-1)+1:36*(I-1)+18,:)=gamm1+gammnewy(36*(I-1)+1:36*(I-1)+18,:);
gammnewy(36*(I-1)+19:36*I,:)=gamm2+gammnewy(36*(I-1)+19:36*I,:);
gammnewy(36*(I-1)+19,:)=1/2*(gamm1+gamm2);
gammnewy(36*I+1,:)=1/2*(gamm1+gamm2);

omegapnewy(36*(I-1)+1:36*(I-1)+18,:)=omegap1+omegapnewy(36*(I-1)+1:36*(I-1)+18,:);
omegapnewy(36*(I-1)+19:36*I,:)=omegap2+omegapnewy(36*(I-1)+19:36*I,:);
omegapnewy(36*(I-1)+19,:)=1/2*(omegap1+omegap2);
omegapnewy(36*I+1,:)=1/2*(omegap1+omegap2);

end
gammnewy(La+1,:)=gamm2;
omegapnewy(La+1,:)=omegap2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CA2x=(1-0.5*gammnewx*deltat)./(1+0.5*gammnewx*deltat);  CD2x=epsilon0*omegapnewx.^2*deltat./(1+0.5*gammnewx*deltat);
CA2y=(1-0.5*gammnewy*deltat)./(1+0.5*gammnewy*deltat);  CD2y=epsilon0*omegapnewy.^2*deltat./(1+0.5*gammnewy*deltat);
%%%%%%%%%%%%%%CA2=(1-0.5*gamm*deltat)/(1+0.5*gamm*deltat);  CD2=epsilon0*omegap^2*deltat/(1+0.5*gamm*deltat);
judgx=ones(La,Lb+1); judgy=ones(La+1,Lb);
judgx(:,49:52)=0*judgx(:,49:52);judgx(:,1)=0.5*judgx(:,1); judgx(:,Lb+1)=0.5*judgx(:,Lb+1);
judgx(:,48)=0.5*judgx(:,48);judgx(:,53)=0.5*judgx(:,53);
judgy(:,48:52)=0*judgy(:,48:52); judgy(1,:)=0.5*judgy(1,:); judgy(La+1,:)=0.5*judgy(La+1,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%金屬范圍的參數設置
aT=ceil(1/f0/deltat); m=20;  T=10000; t0=8*aT;tuot=2*aT;

for t=1:T
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%E update   
if t<m*aT; Xon=1.0-(m*aT-t)/(m*aT); Gon=10*Xon.^3-15*Xon.^4+6*Xon.^5; ft=Gon*sin(omega*t*deltat);
   elseif t>=m*aT; 
  ft=sin(omega*t*deltat); 
  end      %源點引入
  ft2=exp(-(4*pi)*((t-t0)/tuot)^2)*sin(omega*t*deltat);
   
   Ey(25,26:125)=Ey(25,26:125)-0.5*ft; 
    Ex(:,pb+2:Nj-pb)=Ex(:,pb+2:Nj-pb)+(CB/deltay).*(Hz(:,pb+2:Nj-pb)-Hz(:,pb+1:Nj-pb-1)); %update Ex
   Ex(:,2:pb)=CA3xd(:,2:pb).*Ex(:,2:pb)+CB3xd(:,2:pb).*(Hz(:,2:pb)-Hz(:,1:pb-1)); %下邊界 
   Ex(:,pb+1)=CAhalf*Ex(:,pb+1)+(CBhalf).*(Hz(:,pb+1)-Hz(:,pb)); %下邊界連接處
   Ex(:,Nj-pb+2:Nj)=CA3xu(:,2:pb).*Ex(:,Nj-pb+2:Nj)+CB3xu(:,2:pb).*(Hz(:,Nj-pb+2:Nj)-Hz(:,Nj-pb+1:Nj-1));%上邊界
   Ex(:,Nj-pb+1)=CAhalf*Ex(:,Nj-pb+1)+CBhalf*(Hz(:,Nj-pb+1)-Hz(:,Nj-pb)); %上邊界連接處
Ex(:,1)=Ex(:,2);Ex(:,Nj+1)=Ex(:,Nj);
   Ey(pb+2:Ni-pb,:)=Ey(pb+2:Ni-pb,:)-(CB/deltax)*(Hz(pb+2:Ni-pb,:)-Hz(pb+1:Ni-pb-1,:));%update Ey
   Ey(2:pb,:)=CA4yd(2:pb,:).*Ey(2:pb,:)-CB4yd(2:pb,:).*(Hz(2:pb,:)-Hz(1:pb-1,:)); %左邊界
   Ey(pb+1,:)=CAhalf*Ey(pb+1,:)-(CBhalf).*(Hz(pb+1,:)-Hz(pb,:)); %左邊界連接處
   Ey(Ni-pb+2:Ni,:)=CA4yu(2:pb,:).*Ey(Ni-pb+2:Ni,:)-CB4yu(2:pb,:).*(Hz(Ni-pb+2:Ni,:)-Hz(Ni-pb+1:Ni-1,:));%右邊界
   Ey(Ni-pb+1,:)=CAhalf*Ey(Ni-pb+1,:)-CBhalf*(Hz(Ni-pb+1,:)-Hz(Nj-pb,:)); %右邊界連接處
Ey(1,:)=Ey(2,:);Ey(Ni+1,:)=Ey(Ni,:);

 
   
   
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%E update  21068521
Ex(57:344,26:126)=Ex(57:344,26:126)-(deltat/epsilon0)*Jx;
Ey(57:345,26:125)=Ey(57:345,26:125)-(deltat/epsilon0)*Jy;
Jx=CA2x.*Jx+CD2x.*Ex(57:344,26:126).*judgx; 
Jy=CA2y.*Jy+CD2y.*Ey(57:345,26:125).*judgy;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%metal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%H update
Hz1(pb+1:Ni-pb,pb+1:Nj-pb)=Hz(pb+1:Ni-pb,pb+1:Nj-pb)-(CQ/deltay)*((Ey(pb+2:Ni-pb+1,pb+1:Nj-pb)-Ey(pb+1:Ni-pb,pb+1:Nj-pb))-(Ex(pb+1:Ni-pb,pb+2:Nj-pb+1)-Ex(pb+1:Ni-pb,pb+1:Nj-pb)));%update Hz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%hzx update
Hzx(pb+1:Ni-pb,1:pb)=CP*Hzx(pb+1:Ni-pb,1:pb)-(CQ/deltay)*(Ey(pb+2:Ni-pb+1,1:pb)-Ey(pb+1:Ni-pb,1:pb)); %下邊界Hzx
Hzx(pb+1:Ni-pb,Nj-pb+1:Nj)=CP*Hzx(pb+1:Ni-pb,Nj-pb+1:Nj)-(CQ/deltay)*(Ey(pb+2:Ni-pb+1,Nj-pb+1:Nj)-Ey(pb+1:Ni-pb,Nj-pb+1:Nj)); %上邊界Hzx
Hzx(1:pb,:)=CP4xd(1:pb,:).*Hzx(1:pb,:)-CQ4xd(1:pb,:).*((Ey(2:pb+1,:)-Ey(1:pb,:))); %左邊界Hzx
Hzx(Ni-pb+1:Ni,:)=CP4xu(1:pb,:).*Hzx(Ni-pb+1:Ni,:)-CQ4xu(1:pb,:).*((Ey(Ni-pb+2:Ni+1,:)-Ey(Ni-pb+1:Ni,:))); %右邊界Hzx

Hzy(:,1:pb)=CP3yd(:,1:pb).*Hzy(:,1:pb)+CQ3yd(:,1:pb).*((Ex(:,2:pb+1)-Ex(:,1:pb))); %下邊界Hzy
Hzy(:,Nj-pb+1:Nj)=CP3yu(:,1:pb).*Hzy(:,Nj-pb+1:Nj)+CQ3yu(:,1:pb).*((Ex(:,Nj-pb+2:Nj+1)-Ex(:,Nj-pb+1:Nj))); %上邊界Hzy
Hzy(1:pb,pb+1:Nj-pb)=CP*Hzy(1:pb,pb+1:Nj-pb)+(CQ/deltay)*(Ex(1:pb,pb+2:Nj-pb+1)-Ex(1:pb,pb+1:Nj-pb)); %左邊界
Hzy(Ni-pb+1:Ni,pb+1:Nj-pb)=CP*Hzy(Ni-pb+1:Ni,pb+1:Nj-pb)+(CQ/deltay)*(Ex(Ni-pb+1:Ni,pb+2:Nj-pb+1)-Ex(Ni-pb+1:Ni,pb+1:Nj-pb));%右邊界
Hz=Hz1+Hzx+Hzy; 

%Hz(100,30:55)=Hz(100,30:55)+0.5*(sqrt(epsilon0/mu0))*ft;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cichangfenbu=(abs(Hz(32:319,26:125))').^2;
imagesc(cichangfenbu);colorbar;pause(0.001)
end