matlab在光纤中的应用

基于光纤模拟的matlab，程序代码，希望对大家有帮助呀

% this program is designed for solving the two dimensional photonic structure
% we are using the parameters from one of the paper
clear
warning off
epsa=1;
epsb=13;
%epsb=1.46^2; %%%Silica glass air system
a=1.0;
%f=0.39;
%R=sqrt（sqrt（3）.*f/（2*pi）*a^2）;
R=0.45;
i=sqrt（-1）;
f=2*pi/sqrt（3）*R^2/a^2;
%b1=2*pi/a*（1-sqrt（3）/3*i）;
%b2=2*pi/a*2*sqrt（3）/3*i;
b1=2*pi/a*[1 -1/sqrt（3） 0];
b2=2*pi/a*[0 2/sqrt（3） 0];
%n=input（‘please input n: ‘）;
n=5;
NumberofPW=（2*n+1）^2;

count=1;
for x=-n:n
  	for y=-n:n
     	G（count:）=x*b1+y*b2;
      count=count+1;
   end
end

for x=1:NumberofPW
  	for y=x+1:NumberofPW
   	eps2（xy）=（epsa-epsb）*2*f*besselj（1norm（G（x:）-G（y:））*R）./（norm（G（x:）-G（y:））*R）;
   	eps2（yx）=eps2（xy）;   
   end
  	eps2（xx）=f*epsa+（1-f）*epsb;
end
   
%k1=（0:0.1:1.0）/sqrt（3）.*i*2*pi/a;
%k2=（（0.1:0.1:1.0）./3+1/sqrt（3）*i）.*2.*pi./a;
%k3=（0.9:-0.1:0）.*（1.0/3.0+1/sqrt（3）*i）.*2.*pi./a;%-（1/3+1/sqrt（3）*i）*2*pi/a;
%k0=[k1 k2 k3];
%k0=k1;
%k0=（1/3+1/sqrt（3）*i）.*2*pi/a;
k0=zeros（length（0:0.1:1）+length（0.1:0.1:1）+length（0.9:-0.1:0）3）;
mm0=length（0:0.1:1）;
k0（1:mm02）=2*pi/a/sqrt（3）.*（0:0.1:1）‘;
mm0=length（0:0.1:1）;
mm1=length（0.1:0.1:1）;
k0（mm0+1:mm0+mm11）=2*pi/a/3.*（0.1:0.1:1）‘;
k0（mm0+1:mm0+mm12）=2*pi/a/sqrt（3）;
mm0=mm0+mm1+1;
mm1=length（0.9:-0.1:0）;
k0（mm0:mm0+mm1-11）=2*pi/a/3.*（0.9:-0.1:0）‘;
k0（mm0:mm0+mm1-12）=2*pi/a/sqrt（3）.*（0.9:-0.1:0）‘;
k0（:3）=0.0;

%k0=[1/3 1/sqrt（3） 9]*2*pi/a;
counter=1;
eps2=inv（eps2）;
for ii=1:size（k01）
   k=k0（ii:）;
   K（:1）=k（1）+G（:1）;
   K（:2）=k（2）+G（:2）;
   K（:3）=0;
   %%%find the unit cell perpendicular to it in xy plane
   %%%be sure to deal with the case of modulus（k）=0
   %%%NaN in this case
   e1=[K（:2）./modulus（K）-K（:1）./modulus（K）zeros（length（K）1）];
   e1（isnan（e1））=1/sqrt（2）; %%%when Kz is not zerokxky is zero choose arbitrary e1
   K（:3）=k（3）+G（:3）;
   %%%find the other perpendicular unit cell
   %%%be sure to deal with the case of modulus（k）=0
   %%%NaN in this case
   e2=cross（e1 K）;
   e2=[e2（:1）./modulus（e2）e2（:2）./modulus（e2）e2（:3）./modulus（e2）];
   e2（isnan（e2））=0;
   %%%form the equation matrix it should be 2N by 2N
   %%%in this case we will have no TE and TM decoupling
   M1=（[modulus（K）.*e2（:1）modulus（K）.*e2（:2）modulus（K）.*e2（:3）]*[modulus（K）.*e2（:1）modulus（K）.*e2（:2）modulus（K）.*e2（:3）]‘）.*eps2;
   M2=（[modulus（K）.*e1（:1）modulus（K）.*e1（:2）modulus（K）.*e1（:3）]*[modulus（K）.*e2（:1）modulus（K）.*e2（:2）modulus（K）.*e2（:3）]‘）.*eps2;
   M3=（[modulus（K）.*e2（:1）modulus（K）.*e2（:2）modulus（K）.*e2（:3）]*[modulus（K）.*e1（:1）modulus（K）.*e1（:2）modulus（K）.*e1（:3）]‘）.*eps2;
   M4=（[modulus（K）.*e1（:1）modulus（K）.*e1（:2）modulus（K）.*e1（:3）]*[modulus（K）.*e1（:1）modulus（K）.*e1（:2）modulus（K）.*e1（:3）]‘）.*eps2;
   M=[M1M2;M3M4];

	E=sort（abs（eig（M）））;
	freq（:counter）=sqrt（abs（E（1:20）））.*a./2./pi;
	display（sprintf（‘calculation of k=[%f%f%f] is finished‘k（1）k（2）k（3）））;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3279  2003-02-06 08:48  fibermode\bandkz.m

     文件        282  2003-02-06 08:47  fibermode\eigfun.m

     文件       1603  2004-10-13 16:56  fibermode\fibermode.m

     文件       5096  2003-02-06 08:48  fibermode\findmode.m

     文件        268  2003-02-06 08:48  fibermode\fixphase.m

     文件        517  2003-02-06 08:48  fibermode\fld.m

     文件       1024  2007-12-16 23:53  fibermode\Midamble.m

     文件        190  2003-02-06 08:48  fibermode\modulus.m

     文件       1881  2003-02-06 08:48  fibermode\proj.m

     文件       3682  2003-02-06 08:48  fibermode\ptdef5Rfld.m

     文件       3355  2003-02-06 08:49  fibermode\sizetest.m

     文件       2567  2003-02-06 08:49  fibermode\sqdef5.m

     文件       2708  2003-02-06 08:49  fibermode\sqideal.m

     目录          0  2008-03-26 22:22  fibermode

----------- ---------  ---------- -----  ----

                26452                    14