環(huán)形光波導(dǎo)耦合器件程序

10度：%Thisprogramcalculatethewaveguidebyaphotoniccrystal.%Fortwodimensioncase.%OnlyforTMcase（Ez,Hx,Hy）,inthisversion試驗版本.clear;%tic%Initialparametersandotherthings.W=0.36;%Normalizedfrequency就是能帶圖中的縱坐標(biāo)，單位為（a/lamda）%Thefollowingparametersarecontrolparameters.WaveGuide=1;%IfthisprogramisforWaveGuide?Ifso,pleasespecify1.IsMovie=0;%Ifyouwanttoplaymovie,pleaseuse1.IsFigure=1;%Ifitwillplotthefigures?Ifso,pleasespecify1.WantToSeeEp=1;%DoyouwanttoseethedistrubutionofEp（TE極化，S極化）？Ifso,pleasespecify1.%EndofdefiningcontrolparametersMLatx=11;%HowmanyLatticecellinxdirection#格MLaty=11;%HowmanyLatticecellinydirection.NMlat=21;%ThegirdnumberineachLatticeCell.每個柵格中包含的fdtd網(wǎng)格個數(shù)%SHOULDBEODDINTEGER!!!!奇數(shù)整數(shù)個ifmod（NMlat,2）==0NMlat=NMlat+1;end%Forceittobeaoddinteger!NTx=MLatx*NMlat+1;%ItisthenumberoftheGridalongxaxis總的fdtd網(wǎng)格坐標(biāo)點（所以要加1）NTy=MLaty*NMlat+1;%ItisthenumberoftheGridalongyaxis.ifWaveGuide==1Nrow=5;%TherownumberofcoloumnsbetweenthePMLboundaryandthewaveguide.（此處定義，現(xiàn)在還不知所措）endNPML=12;%HowmanyPMLlayerswillbeusedinourcomputation.PML的層數(shù)NTimeSteps=5000;%TotalnumberofTimeStepsMeach=20;%DefinetheintervalforplotfiguresifIsFigure==1.%Thisalsoworksforsavingintervals.每隔20個時間步顯示一次。Zmax=0.8;%Themaximumvalueforzaxiswhenplottingfigures.Colormax=2.8;%Themaximumvalueforcolormapwhenplottingfigures.%Someconstantsmu0=4*pi*1.0e-7;%EpsilonZero,ifusingGaussUnit,itequalsto1.e0=8.85*1e-12;%MuZero,ifusingGaussUnit,itequalsto1.c=1/sqrt(mu0*e0);%Thelightspeed.factor=mu0/e0;%Thefactorbetweenconductivityandpermeability.%Permeability=Conductivity*factor,inPML.PML中的匹配條件a=1;%e-6; %Thelatticeconstant.柵格常數(shù)。W=W*(2*pi*c/a);%frequencyDx=a/NMlat;%Deltax.dx=1/(一個柵格中的網(wǎng)格數(shù))Dy=Dx;%Deltay.Dt=1/sqrt(1/(Dx*Dx)+1/(Dy*Dy))/c;%Timeinterval%tic%Inthefollowingpartmwedefinethedielectric%constants:,定義FDTD中平面波導(dǎo)與空氣的介電常數(shù)Ep=ones(NTx-1,NTy-1)*e0;%FDTD總的網(wǎng)格中的介質(zhì)常數(shù)，真空。Ep_cell=ones(NMlat,NMlat)*e0;%每個柵格中的fdtd介電常數(shù)，真空。x=-(NMlat-1)/2*Dx:Dx:(NMlat-1)/2*Dx;%先對柵格操作，然后進(jìn)行repmat,復(fù)制并平鋪陣列。[X,Y]=meshgrid(x);%把x方向上的向量拓展為X,Y兩個矩陣。X=X';Y=Y';flag=find(sqrt(X.A2+Y.A2)<R);%Ep_cell(flag)=e0*ea;Ep=repmat(Ep_cell,MLatx,MLaty);%拓展到真?zhèn)€空間了。angle=pi/18;%夾角，可修改；ifWaveGuide==1eb=1;%ThedielectricconstantinthewaveguideEp((MLatx-Nrow-6)*NMlat+1:(MLatx-Nrow+5)*NMlat,(Nrow-5)*NMlat+1:(Nrow+6)*NMlat)=e0*eb;Ep(1:(MLatx-Nrow-1)*NMlat,Nrow*NMlat+1:(Nrow+1)*NMlat)=e0*ea;length=round((NMlat+1)/2*cot(angle));j=(MLaty*NMlat+1)/2;fori=1:cot(angle):(MLatx*NMlat+1)/2+length+(NMlat+1)/2k=round(i);left=(MLaty*NMlat+1)/2-(NMlat+1)/2-length+i-1;%常數(shù)可修改right=left+length+55;%常數(shù)可修改if(right>MLatx*NMlat)right=MLatx*NMlat;endEp(left:right,j)=e0*ea;j=j-1;endj=(MLaty*NMlat+1)/2;fori=1:cot(angle):(MLatx*NMlat+1)/2+length+(NMlat+1)/2k=round(i);left=(MLaty*NMlat+1)/2-(NMlat+1)/2-length+i-1;%常數(shù)可修改right=left+length+55;%常數(shù)可修改if(right>MLatx*NMlat)right=MLatx*NMlat;endEp(left:right,j)=e0*ea;j=j+1;endend%tocifWantToSeeEp==1x=0:Dx:(NMlat*MLatx-1)*Dx;x=x-(NMlat*MLatx-1)*Dx/2;y=0:Dy:(NMlat*MLaty-1)*Dy;y=y-(NMlat*MLaty-1)*Dy/2;[X,Y]=meshgrid(x,y);X=X';Y=Y';surf(X,Y,Ep/e0);shadinginterp;view(0,90);axis([min(x),max(x),min(y),max(y)])axisoff;disp('Pressanykeytocontinue...');pauseend%EndofdefiningtheEp.15度：%Thisprogramcalculatethewaveguidebyaphotoniccrystal.%Fortwodimensioncase.%OnlyforTMcase（Ez,Hx,Hy）,inthisversion試驗版本.clear;%tic%Initialparametersandotherthings.W=0.36;%Normalizedfrequency就是能帶圖中的縱坐標(biāo)，單位為（a/lamda）%Thefollowingparametersarecontrolparameters.WaveGuide=1;%IfthisprogramisforWaveGuide?Ifso,pleasespecify1.IsMovie=0;%Ifyouwanttoplaymovie,pleaseuse1.IsFigure=1;%Ifitwillplotthefigures?Ifso,pleasespecify1.WantToSeeEp=1;%DoyouwanttoseethedistrubutionofEp（TE極化，S極化）？Ifso,pleasespecify1.%EndofdefiningcontrolparametersMLatx=11;%HowmanyLatticecellinxdirection#格MLaty=11;%HowmanyLatticecellinydirection.NMlat=21;%ThegirdnumberineachLatticeCell.每個柵格中包含的fdtd網(wǎng)格個數(shù)%SHOULDBEODDINTEGER!!!!奇數(shù)整數(shù)個ifmod（NMlat,2）==0NMlat=NMlat+1;end%Forceittobeaoddinteger!NTx=MLatx*NMlat+1;%ItisthenumberoftheGridalongxaxis總的fdtd網(wǎng)格坐標(biāo)點（所以要加1）NTy=MLaty*NMlat+1;%ItisthenumberoftheGridalongyaxis.ifWaveGuide==1Nrow=5;%TherownumberofcoloumnsbetweenthePMLboundaryandthewaveguide.（此處定義，現(xiàn)在還不知所措）endNPML=12;%HowmanyPMLlayerswillbeusedinourcomputation.PML的層數(shù)NTimeSteps=5000;%TotalnumberofTimeStepsMeach=20;%DefinetheintervalforplotfiguresifIsFigure==1.%Thisalsoworksforsavingintervals.每隔20個時間步顯示一次。Zmax=0.8;%Themaximumvalueforzaxiswhenplottingfigures.Colormax=2.8;%Themaximumvalueforcolormapwhenplottingfigures.%Someconstantsmu0=4*pi*1.0e-7;%EpsilonZero,ifusingGaussUnit,itequalsto1.e0=8.85*1e-12;%MuZero,ifusingGaussUnit,itequalsto1.c=1/sqrt(mu0*e0);%Thelightspeed.factor=mu0/e0;%Thefactorbetweenconductivityandpermeability.%Permeability=Conductivity*factor,inPML.PML中的匹配條件a=1;%e-6; %Thelatticeconstant.柵格常數(shù)。W=W*(2*pi*c/a);%frequencyDx=a/NMlat;%Deltax.dx=1/(一個柵格中的網(wǎng)格數(shù))Dy=Dx;%Deltay.Dt=1/sqrt(1/(Dx*Dx)+1/(Dy*Dy))/c;%Timeinterval%tic%Inthefollowingpartmwedefinethedielectric%constants:,定義FDTD中平面波導(dǎo)與空氣的介電常數(shù)Ep=ones(NTx-1,NTy-1)*e0;%FDTD總的網(wǎng)格中的介質(zhì)常數(shù)，真空。Ep_cell=ones(NMlat,NMlat)*e0;%每個柵格中的fdtd介電常數(shù)，真空。x=-(NMlat-1)/2*Dx:Dx:(NMlat-1)/2*Dx;%先對柵格操作，然后進(jìn)行repmat,復(fù)制并平鋪陣列。[X,Y]=meshgrid(x);%把x方向上的向量拓展為X,Y兩個矩陣。X=X';Y=Y';flag=find(sqrt(X.A2+Y.A2)<R);%Ep_cell(flag)=e0*ea;Ep=repmat(Ep_cell,MLatx,MLaty);%拓展到真?zhèn)€空間了。angle=pi/12;%夾角，可修改；ifWaveGuide==1eb=1;%ThedielectricconstantinthewaveguideEp((MLatx-Nrow-6)*NMlat+1:(MLatx-Nrow+5)*NMlat,(Nrow-5)*NMlat+1:(Nrow+6)*NMlat)=e0*eb;Ep(1:(MLatx-Nrow-1)*NMlat,Nrow*NMlat+1:(Nrow+1)*NMlat)=e0*ea;length=round((NMlat+1)/2*cot(angle));j=(MLaty*NMlat+1)/2;fori=1:cot(angle):(MLatx*NMlat+1)/2+length+(NMlat+1)/2k=round(i);left=(MLaty*NMlat+1)/2-(NMlat+1)/2-length+i-1;%常數(shù)可修改right=left+length+40;%常數(shù)可修改if(right>MLatx*NMlat)right=MLatx*NMlat;endEp(left:right,j)=e0*ea;j=j-1;endj=(MLaty*NMlat+1)/2;fori=1:cot(angle):(MLatx*NMlat+1)/2+length+(NMlat+1)/2k=round(i);left=(MLaty*NMlat+1)/2-(NMlat+1)/2-length+i-1;%常數(shù)可修改right=left+length+40;%常數(shù)可修改if(right>MLatx*NMlat)right=MLatx*NMlat;endEp(left:right,j)=e0*ea;j=j+1;endend%tocifWantToSeeEp==1x=0:Dx:(NMlat*MLatx-1)*Dx;x=x-(NMlat*MLatx-1)*Dx/2;y=0:Dy:(NMlat*MLaty-1)*Dy;y=y-(NMlat*MLaty-1)*Dy/2;[X,Y]=meshgrid(x,y);X=X';Y=Y';surf(X,Y,Ep/e0);shadinginterp;view(0,90);axis([min(x),max(x),min(y),max(y)])axisoff;disp('Pressanykeytocontinue...');pauseend%EndofdefiningtheEp.20度：%Thisprogramcalculatethewaveguidebyaphotoniccrystal.%Fortwodimensioncase.%OnlyforTMcase（Ez,Hx,Hy）,inthisversion試驗版本.clear;%tic%Initialparametersandotherthings.W=0.36;%Normalizedfrequency就是能帶圖中的縱坐標(biāo)，單位為（a/lamda）%Thefollowingparametersarecontrolparameters.WaveGuide=1;%IfthisprogramisforWaveGuide?Ifso,pleasespecify1.IsMovie=0;%Ifyouwanttoplaymovie,pleaseuse1.IsFigure=1;%Ifitwillplotthefigures?Ifso,pleasespecify1.WantToSeeEp=1;%DoyouwanttoseethedistrubutionofEp（TE極化，S極化）？Ifso,pleasespecify1.%EndofdefiningcontrolparametersMLatx=11;%HowmanyLatticecellinxdirection#格MLaty=11;%HowmanyLatticecellinydirection.NMlat=21;%ThegirdnumberineachLatticeCell.每個柵格中包含的fdtd網(wǎng)格個數(shù)%SHOULDBEODDINTEGER!!!!奇數(shù)整數(shù)個ifmod（NMlat,2）==0NMlat=NMlat+1;end%Forceittobeaoddinteger!NTx=MLatx*NMlat+1;%ItisthenumberoftheGridalongxaxis總的fdtd網(wǎng)格坐標(biāo)點（所以要加1）NTy=MLaty*NMlat+1;%ItisthenumberoftheGridalongyaxis.ifWaveGuide==1Nrow=5;%TherownumberofcoloumnsbetweenthePMLboundaryandthewaveguide.（此處定義，現(xiàn)在還不知所措）endNPML=12;%HowmanyPMLlayerswillbeusedinourcomputation.PML的層數(shù)NTimeSteps=5000;%TotalnumberofTimeStepsMeach=20;%DefinetheintervalforplotfiguresifIsFigure==1.%Thisalsoworksforsavingintervals.每隔20個時間步顯示一次。Zmax=0.8;%Themaximumvalueforzaxiswhenplottingfigures.Colormax=2.8;%Themaximumvalueforcolormapwhenplottingfigures.%Someconstantsmu0=4*pi*1.0e-7;%EpsilonZero,ifusingGaussUnit,itequalsto1.e0=8.85*1e-12;%MuZero,ifusingGaussUnit,itequalsto1.c=1/sqrt(mu0*e0);%Thelightspeed.factor=mu0/e0;%Thefactorbetweenconductivityandpermeability.%Permeability=Conductivity*factor,inPML.PML中的匹配條件a=1;%e-6; %Thelatticeconstant.柵格常數(shù)。W=W*(2*pi*c/a);%frequencyDx=a/NMlat;%Deltax.dx=1/(一個柵格中的網(wǎng)格數(shù))Dy=Dx;%Deltay.Dt=1/sqrt(1/(Dx*Dx)+1/(Dy*Dy))/c;%Timeinterval%tic%Inthefollowingpartmwedefinethedielectric%constants:,定義FDTD中平面波導(dǎo)與空氣的介電常數(shù)Ep=ones(NTx-1,NTy-1)*e0;%FDTD總的網(wǎng)格中的介質(zhì)常數(shù)，真空。Ep_cell=ones(NMlat,NMlat)*e0;%每個柵格中的fdtd介電常數(shù)，真空。x=-(NMlat-1)/2*Dx:Dx:(NMlat-1)/2*Dx;%先對柵格操作，然后進(jìn)行repmat,復(fù)制并平鋪陣列。[X,Y]=meshgrid(x);%把x方向上的向量拓展為X,Y兩個矩陣。X=X';Y=Y';flag=find(sqrt(X.A2+Y.A2)<R);%Ep_cell(flag)=e0*ea;Ep=repmat(Ep_cell,MLatx,MLaty);%拓展到真?zhèn)€空間了。angle=pi/9;%夾角，可修改；ifWaveGuide==1eb=1;%ThedielectricconstantinthewaveguideEp((MLatx-Nrow-6)*NMlat+1:(MLatx-Nrow+5)*NMlat,(Nrow-5)*NMlat+1:(Nrow+6)*NMlat)=e0*eb;Ep(1:(MLatx-Nrow-1)*NMlat,Nrow*NMlat+1:(Nrow+1)*NMlat)=e0*ea;length=round((NMlat+1)/2*cot(angle));j=(MLaty*NMlat+1)/2;fori=1:cot(angle):(MLatx*NMlat+1)/2+length+(NMlat+1)/2k=round(i);left=(MLaty*NMlat+1)/2-(NMlat+1)/2-length+i-1;%常數(shù)可修改right=left+length+30;%常數(shù)可修改if(right>MLatx*NMlat)right=MLatx*NMlat;endEp(left:right,j)=e0*ea;j=j-1;endj=(MLaty*NMlat+1)/2;fori=1:cot(angle):(MLatx*NMlat+1)/2+length+(NMlat+1)/2k=round(i);left=(MLaty*NMlat+1)/2-(NMlat+1)/2-length+i-1;%常數(shù)可修改right=left+length+30;%常數(shù)可修改if(right>MLatx*NMlat)right=MLatx*NMlat;endEp(left:right,j)=e0*ea;j=j+1;endend%tocifWantToSeeEp==1x=0:Dx:(NMlat*MLatx-1)*Dx;x=x-(NMlat*MLatx-1)*Dx/2;y=0:Dy:(NMlat*MLaty-1)*Dy;y=y-(NMlat*MLaty-1)*Dy/2;[X,Y]=meshgrid(x,y);X=X';Y=Y';surf(X,Y,Ep/e0);shadinginterp;view(0,90);axis([min(x),max(x),min(y),max(y)])axisoff;disp('Pressanykeytocontinue...');pauseend%EndofdefiningtheEp.30度：%Thisprogramcalculatethewaveguidebyaphotoniccrystal.%Fortwodimensioncase.%OnlyforTMcase（Ez,Hx,Hy）,inthisversion試驗版本.clear;%tic%Initialparametersandotherthings.W=0.36;%Normalizedfrequency就是能帶圖中的縱坐標(biāo)，單位為（a/lamda）%Thefollowingparametersarecontrolparameters.WaveGuide=1;%IfthisprogramisforWaveGuide?Ifso,pleasespecify1.IsMovie=0;%Ifyouwanttoplaymovie,pleaseuse1.IsFigure=1;%Ifitwillplotthefigures?Ifso,pleasespecify1.WantToSeeEp=1;%DoyouwanttoseethedistrubutionofEp（TE極化，S極化）？Ifso,pleasespecify1.%EndofdefiningcontrolparametersMLatx=11;%HowmanyLatticecellinxdirection#格MLaty=11;%HowmanyLatticecellinydirection.NMlat=21;%ThegirdnumberineachLatticeCell.每個柵格中包含的fdtd網(wǎng)格個數(shù)%SHOULDBEODDINTEGER!!!!奇數(shù)整數(shù)個ifmod（NMlat,2）==0NMlat=NMlat+1;end%Forceittobeaoddinteger!NTx=MLatx*NMlat+1;%ItisthenumberoftheGridalongxaxis總的fdtd網(wǎng)格坐標(biāo)點（所以要加1）NTy=MLaty*NMlat+1;%ItisthenumberoftheGridalongyaxis.ifWaveGuide==1Nrow=5;%TherownumberofcoloumnsbetweenthePMLboundaryandthewaveguide.（此處定義，現(xiàn)在還不知所措）endNPML=12;%HowmanyPMLlayerswillbeusedinourcomputation.PML的層數(shù)NTimeSteps=5000;%TotalnumberofTimeStepsMeach=20;%DefinetheintervalforplotfiguresifIsFigure==1.%Thisalsoworksforsavingintervals.每隔20個時間步顯示一次。Zmax=0.8;%Themaximumvalueforzaxiswhenplottingfigures.Colormax=2.8;%Themaximumvalueforcolormapwhenplottingfigures.%Someconstantsmu0=4*pi*1.0e-7;%EpsilonZero,ifusingGaussUnit,itequalsto1.e0=8.85*1e-12;%MuZero,ifusingGaussUnit,itequalsto1.c=1/sqrt(mu0*e0);%Thelightspeed.factor=mu0/e0;%Thefactorbetweenconductivityandpermeability.%Permeability=Conductivity*factor,inPML.PML中的匹配條件a=1;%e-6; %Thelatticeconstant.柵格常數(shù)。W=W*(2*pi*c/a);%frequencyDx=a/NMlat;%Deltax.dx=1/(一個柵格中的網(wǎng)格數(shù))Dy=Dx;%Deltay.Dt=1/sqrt(1/(Dx*Dx)+1/(Dy*Dy))/c;%Timeinterval%tic%Inthefollowingpartmwedefinethedielectric%constants:,定義FDTD中平面波導(dǎo)與空氣的介電常數(shù)Ep=ones(NTx-1,NTy-1)*e0;%FDTD總的網(wǎng)格中的介質(zhì)常數(shù)，真空。Ep_cell=ones(NMlat,NMlat)*e0;%每個柵格中的fdtd介電常數(shù)，真空。x=-(NMlat-1)/2*Dx:Dx:(NMlat-1)/2*Dx;%先對柵格操作，然后進(jìn)行repmat,復(fù)制并平鋪陣列。[X,Y]=meshgrid(x);%把x方向上的向量拓展為X,Y兩個矩陣。X=X';Y=Y';flag=find(sqrt(X.A2+Y.A2)<R);%Ep_cell(flag)=e0*ea;Ep=repmat(Ep_cell,MLatx,MLaty);%拓展到真?zhèn)€空間了。angle=pi/6;%夾角，可修改；ifWaveGuide==1eb=1;%ThedielectricconstantinthewaveguideEp((MLatx-Nrow-6)*NMlat+1:(MLatx-Nrow+5)*NMlat,(Nrow-5)*NMlat+1:(Nrow+6)*NMlat)=e0*eb;Ep(1:(MLatx-Nrow-1)*NMlat,Nrow*NMlat+1:(Nrow+1)*NMlat)=e0*ea;length=round((NMlat+1)/2*cot(angle));j=(MLaty*NMlat+1)/2;fori=1:cot(angle):(MLatx*NMlat+1)/2+length+(NMlat+1)/2k=round(i);left=(MLaty*NMlat+1)/2-(NMlat+1)/2-length+i-1;%常數(shù)可修改right=left+length+20;%常數(shù)可修改if(right>MLatx*NMlat)right=MLatx*NMlat;endEp(left:right,j)=e0*ea;j=j-1;endj=(MLaty*NMlat+1)/2;fori=1:cot(angle):(MLatx*NMlat+1)/2+length+(NMlat+1)/2k=round(i);left=(MLaty*NMlat+1)/2-(NMlat+1)/2-length+i-1;%常數(shù)可修改right=left+length+20;%常數(shù)可修改if(right>MLatx*NMlat)right=MLatx*NMlat;endEp(left:right,j)=e0*ea;j=j+1;endend%tocifWantToSeeEp==1x=0:Dx:(NMlat*MLatx-1)*Dx;x=x-(NMlat*MLatx-1)*Dx/2;y=0:Dy:(NMlat*MLaty-1)*Dy;y=y-(NMlat*MLaty-1)*Dy/2;[X,Y]=meshgrid(x,y);X=X';Y=Y';surf(X,Y,Ep/e0);shadinginterp;view(0,90);axis([min(x),max(x),min(y),max(y)])axisoff;disp('Pressanykeytocontinue...');pauseend%EndofdefiningtheEp.80度：%Thisprogramcalculatethewaveguidebyaphotoniccrystal.%Fortwodimensioncase.%OnlyforTMcase（Ez,Hx,Hy）,inthisversion試驗版本.clear;%tic%Initialparametersandotherthings.W=0.36;%Normalizedfrequency就是能帶圖中的縱坐標(biāo)，單位為（a/lamda）%Thefollowingparametersarecontrolparameters.WaveGuide=1;%IfthisprogramisforWaveGuide?Ifso,pleasespecify1.IsMovie=0;%Ifyouwanttoplaymovie,pleaseuse1.IsFigure=1;%Ifitwillplotthefigures?Ifso,pleasespecify1.WantToSeeEp=1;%DoyouwanttoseethedistrubutionofEp（TE極化，S極化）？Ifso,pleasespecify1.%EndofdefiningcontrolparametersMLatx=11;%HowmanyLatticecellinxdirection光子晶體柵格MLaty=11;%HowmanyLatticecellinydirection.NMlat=21;%ThegirdnumberineachLatticeCell.每個柵格中包含的fdtd網(wǎng)格個數(shù)%SHOULDBEODDINTEGER!!!!奇數(shù)整數(shù)個,目的就是為了將光子晶體的原點定位在網(wǎng)格中心ifmod（NMlat,2）==0NMlat=NMlat+1;end%Forceittobeaoddinteger!NTx=MLatx*NMlat+1;%ItisthenumberoftheGridalongxaxis總的fdtd網(wǎng)格坐標(biāo)點NTy=MLaty*NMlat+1;%ItisthenumberoftheGridalongyaxis.ifWaveGuide==1Nrow=5;%TherownumberofcoloumnsbetweenthePMLboundaryandthewaveguide.（此處定義，現(xiàn)在還不知所措）endNPML=12;%HowmanyPMLlayerswillbeusedinourcomputation.PML的層數(shù)NTimeSteps=4000;%TotalnumberofTimeStepsMeach=20;%DefinetheintervalforplotfiguresifIsFigure==1.%Thisalsoworksforsavingintervals.每隔20個時間步顯示一次。R=0.2;%Theradiusofdielectriccolumns，光子晶體的半徑，可調(diào)Zmax=0.6;%Themaximumvalueforzaxiswhenplottingfigures.Colormax=1.6;%Themaximumvalueforcolormapwhenplottingfigures.%Someconstantsmu0=4*pi*1.0e-7;%EpsilonZero,ifusingGauss