P.R.E 67(046607) (2003) & J.A.P.94(2) 811 (2003)
程序名：TMM（Plane Wave Based Transfer Matrix Method）
常數(shù)：
pi=3.141592653589793 圓周率
error=0.001 相對誤差 用于判斷計(jì)算結(jié)果的準(zhǔn)確性
m=15:每一個(gè)層單元（Layer）的離散層數(shù)（Slice）
輸入?yún)?shù)：
a:REAL(8)，晶格常數(shù)(micrometer)
f：REAL(8)，填充率 圓形的介質(zhì)柱／空氣孔半徑與晶格常數(shù)的比值r/a 范圍：0.0～0.5
die1：REAL(8)，周期介支柱／空氣孔介電常數(shù)(相對值)
die2:REAL(8)，背景的介電常數(shù)(相對值)
theta_k，phi_k:REAL(8)，入射波矢的極化角和方位角
theta_E，phi_E:REAL(8)，入射波的電場矢量的極化角和方位角

內(nèi)部參數(shù)：
a1,a2:REAL(8)，分別為TM,TK方向離散化得到的光柵(slice)中的周期長度
b1,b2:REAL(8)，用a1,a2表示的周期光柵的的Bragg波矢(倒格矢)
c1,c2:REAL(8)，c1是每層光柵中對應(yīng)于介質(zhì)柱或空氣孔部分的寬度；c2是當(dāng)Layer與Layer發(fā)生重疊時(shí)屬于相鄰層單元
              的介質(zhì)柱或空氣孔部分的寬度
h:REAL(8)，每一層光柵的厚度
sh_i_1:REAL(8)，臨時(shí)存放層單元重疊部分的厚度，用于判斷重疊部分是否離散完畢
y01,y02:REAL(8)，兩種情況下相鄰光柵的平移距離
theta_i:REAL(8)，圓柱／圓孔離散化過程中采用角度平均得到的第i個(gè)離散化角
ko,kox,koy,koz：REAL(8)，入射波真空波矢，z方向?yàn)椴▊鞑シ较颍瑈在周期結(jié)構(gòu)平面內(nèi)，
                          x平行于介支柱／空氣孔
Eox,Eoy,Eoz:REAL(8)，入射電磁波的電場分量
RT_result:REAL(8)，2000*3矩陣，用于存放波數(shù)(第1列)反射率(第2列)透射率(第3列)
temp:REAL(8)，反射率／透射率的相對誤差
sum:REAL(8)，所有頻率處的反射率／透射率的總和
temp_result:REAL(8)，2000*2矩陣，用于臨時(shí)存放反射率（第1列）／透射率（第2列）
k0:INTEGER，頻率域的離散點(diǎn)
n:INTEGER，EF0,EFz等的維數(shù)，包括E的X,Y分量的Bragg波矢數(shù)目
nn:INTEGER，E的X或Y分量的Gragg波矢數(shù)目
i:INTEGER，循環(huán)變量
m_unit:INTEGER，在TM和TK方向的重復(fù)單元數(shù)
RT:INTEGER，RT=1計(jì)算反射譜，RT=2計(jì)算透射譜
t1:REAL(8),n*n矩陣，物理意義同文獻(xiàn),是E的一階微分方程與磁場的系數(shù)矩陣
p:REAL(8),n*n矩陣，物理意義同文獻(xiàn)，是T1和T2矩陣的乘積
sa:COMPLEX(8),n*n矩陣，光柵層中的電場E與向前和向后傳播的電場的系數(shù)矩陣
ta:COMPLEX(8),n*n矩陣，光柵層中的磁場H與向前和向后傳播的電場的系數(shù)矩陣
beta:COMPLEX(8),n*n矩陣，傳輸系數(shù)(本征值)矩陣，對角元是傳輸系數(shù)
S0:COMPLEX(8),n*n矩陣，空氣層中的電場E與向前和向后傳播的電場的系數(shù)矩陣,物理意義同sa
T0:COMPLEX(8),n*n矩陣，空氣層中的電場H與向前和向后傳播的電場的系數(shù)矩陣,物理意義同ta
s:COMPLEX(8),2n*2n矩陣，單層光柵層(Slice)的傳遞矩陣(散射矩陣S)
layer1:COMPLEX(8),2n*2n矩陣，第1個(gè)層單元(Layer)的傳遞矩陣(散射矩陣S)
layer2:COMPLEX(8),2n*2n矩陣，第2個(gè)層單元(Layer)的傳遞矩陣(散射矩陣S)
slab:COMPLEX(8),2n*2n矩陣，重復(fù)周期(unit,有2個(gè)Layer組成)的傳遞矩陣(散射矩陣S)
EF0:COMPLEX(8),n*1向量，入射電磁波的列向量
EFn:COMPLEX(8),n*1向量，反射／透射電磁波的列向量(由RT決定)
Ez:COMPLEX(8),nn*1向量，反射／透射電磁波的Z分量的列向量(由RT決定)

外部函數(shù)：
RT_fun(ko,kox,koy,koz,b1,n,nn,EFn,Ez)子函數(shù)。
用于計(jì)算反射率(RT=1)和透射率(RT=2)
輸入?yún)?shù)如上所述
輸出:REAL(8)，RT=1為反射率；RT=2為透射率

overlap(i,a,f,m,h,sh_i_1,c1,c2,KM)子程序
用于計(jì)算c1,c2,h和判別圓孔／圓柱離散是否結(jié)束
輸入?yún)?shù)
i:INTEGER,離散化過程中的循環(huán)次數(shù)
a:REAL(8)，光柵周期
f:REAL(8)，填充率
m:INTEGER，離散化過程中的光柵離散總數(shù)
KM:INTEGER,KM=1計(jì)算TM方向，KM=2計(jì)算TK方向
輸出參數(shù)
h:REAL(8)，光柵層厚度
sh_i_1:REAL(8)，臨時(shí)存放層單元重疊部分的厚度，用于判斷重疊部分是否離散完畢
c1,c2:REAL(8)，c1是每層光柵中對應(yīng)于介質(zhì)柱或空氣孔部分的寬度；c2是當(dāng)Layer與Layer發(fā)生重疊時(shí)屬于相鄰層單元
              的介質(zhì)柱或空氣孔部分的寬度

Ez_array(ko,kox,koy,b,n,nn,EF0,slab,T0,EFn,Ez,RT)子程序
用于計(jì)算反射／透射電磁波的電場的X,Y,Z分量
輸入?yún)?shù)    
ko,kox,koy:REAL(8),入射波真空波矢，z方向?yàn)椴▊鞑シ较颍瑈在周期結(jié)構(gòu)平面內(nèi)，
                    x平行于介支柱／空氣孔
b:REAL(8),光柵周期的倒格矢
n,nn:INTEGER,矩陣維數(shù)和Bragg階數(shù)
EF0:COMPLEX(8),n*1向量，入射電磁波的列向量
slab:COMPLEX(8),2n*2n矩陣，重復(fù)周期(unit,有2個(gè)Layer組成)的傳遞矩陣(散射矩陣S)
T0:COMPLEX(8),n*n矩陣，空氣層中的電場H與向前和向后傳播的電場的系數(shù)矩陣,物理意義同ta
輸出參數(shù)
EFn:COMPLEX(8),n*1向量，反射／透射電磁波的列向量(由RT決定)
Ez:nn*1向量，反射／透射電磁波的Z分量的列向量(由RT決定)

s_layer2(n,layer1,b,y0,layer2)子程序
用于計(jì)算第2層Layer的散射矩陣
輸入?yún)?shù)
n:INTEGER,矩陣維數(shù)
b:REAL(8),光柵周期的倒格矢
y0:REAL(8),Y方向的平移距離
layer1:COMPLEX(8),2n*2n矩陣，第1層層單元的散射矩陣
輸出參數(shù)
layer2:COMPLEX(8),2n*2n矩陣，第2層層單元的散射矩陣

s_layer1(n,s,layer)子矩陣
通過散射矩陣S的迭代關(guān)系計(jì)算若干層的傳遞矩陣
輸入?yún)?shù)
n:INTEGER,矩陣維數(shù)
s:COMPLEX(8),2n*2n矩陣，單層光柵／層單元的傳遞矩陣
輸入輸出參數(shù)
layer:COMPLEX(8),2n*2n矩陣，輸入時(shí)為上一次迭代的結(jié)果，表示前i層的傳遞矩陣，輸出時(shí)為
      迭代后的傳遞矩陣

s_matrix(h,n,S0,T0,sa,ta,beta,s)子程序
用于計(jì)算傳遞矩陣(散射矩陣S)
輸入?yún)?shù)
h:REAL(8),光柵層(Slice)厚度
n:INTEGER,矩陣維數(shù)
S0:COMPLEX(8),n*n矩陣，空氣層中的電場E與向前和向后傳播的電場的系數(shù)矩陣,物理意義同sa
T0:COMPLEX(8),n*n矩陣，空氣層中的電場H與向前和向后傳播的電場的系數(shù)矩陣,物理意義同ta
sa:COMPLEX(8),n*n矩陣，光柵層中的電場E與向前和向后傳播的電場的系數(shù)矩陣，也是傳遞矩陣的子矩陣
ta:COMPLEX(8),n*n矩陣，光柵層中的磁場H與向前和向后傳播的電場的系數(shù)矩陣，也是傳遞矩陣的子矩陣
輸入輸出參數(shù)
beta:COMPLEX(8),n*n矩陣，輸入時(shí)為本征值矩陣，輸出時(shí)為e指數(shù)，表示相應(yīng)矩陣元上的位相因子
輸出參數(shù)
s:COMPLEX(8),2n*2n矩陣，層單元／整個(gè)光子晶體的傳遞矩陣

ST0(b,kox,koy,ko,n,S0,T0)子程序
用于計(jì)算空氣層的傳遞矩陣的子矩陣S0,T0(物理意義如前述)
輸入?yún)?shù)
b:REAL(8),光柵周期的倒格矢
kox,koy,ko:REAL(8)，入射電磁波的波矢
n:INTEGER,矩陣維數(shù)
輸出參數(shù)
S0:COMPLEX(8),n*n矩陣，空氣層中的電場E與向前和向后傳播的電場的系數(shù)矩陣,物理意義同sa
T0:COMPLEX(8),n*n矩陣，空氣層中的電場H與向前和向后傳播的電場的系數(shù)矩陣,物理意義同ta

STa(ko,t1,p,beta,sa,ta,n)子程序
用于計(jì)算光柵層中的傳遞矩陣的子矩陣sa,ta(物理意義如前述)
輸入?yún)?shù)
ko:REAL(8)，入射電磁波的波矢
t1:REAL(8),n*n矩陣，物理意義同文獻(xiàn),是E的一階微分方程與磁場的系數(shù)矩陣
p:REAL(8),n*n矩陣，物理意義同文獻(xiàn)，是T1和T2矩陣的乘積
n:INTEGER,矩陣維數(shù)
輸出參數(shù)
beta:COMPLEX(8),n*n矩陣，傳輸系數(shù)(本征值)矩陣，對角元是傳輸系數(shù)
sa:COMPLEX(8),n*n矩陣，光柵層中的電場E與向前和向后傳播的電場的系數(shù)矩陣，也是傳遞矩陣的子矩陣
ta:COMPLEX(8),n*n矩陣，光柵層中的磁場H與向前和向后傳播的電場的系數(shù)矩陣，也是傳遞矩陣的子矩陣

construc(a,c1,c2,b,kox,koy,ko,die1,die2,n,t1,p)子程序
用于構(gòu)筑系數(shù)矩陣t1,本征矩陣p
輸入?yún)?shù)
a:REAL(8)，光柵周期
c1:REAL(8)，每層光柵中對應(yīng)于介質(zhì)柱或空氣孔部分的寬度；
c2:REAL(8)，當(dāng)Layer與Layer發(fā)生重疊時(shí)屬于相鄰層單元的介質(zhì)柱或空氣孔部分的寬度
b:REAL(8)，光柵周期的倒格矢
ko,kox,koy:REAL(8)，入射電磁波的波矢
die1:REAL(8)，周期圓柱／圓孔的介電常數(shù)
die2:REAL(8)，基體的介電常數(shù)
n:INTEGER,矩陣維數(shù)
輸出參數(shù)
t1:REAL(8),n*n矩陣，物理意義同文獻(xiàn),是E的一階微分方程與磁場的系數(shù)矩陣
p:REAL(8),n*n矩陣，物理意義同文獻(xiàn)，是T1和T2矩陣的乘積

delta(i,j)子函數(shù)
用于計(jì)算delta函數(shù)
輸入?yún)?shù)：
i:INTEGER，Bragg波矢階數(shù)
j:INTEGER，Bragg波矢階數(shù)
輸出參數(shù)：
delta:REAL(8)，i=j為1.0；i/=j為0.0

ftrans(a,b,c1,c2,die1,die2,k,l,n)子函數(shù)
用于計(jì)算光柵周期分布的介電函數(shù)的Fourier系數(shù)
輸入?yún)?shù)：
a:REAL(8)，光柵周期
b:REAL(8)，光柵的Bragg波矢
c1:REAL(8)，周期圓柱／圓孔在光柵(slice)中的寬度
C2:REAL(8)，相鄰光柵的周期圓柱／圓孔在本光柵中的交疊部分的寬度
die1:REAL(8)，周期圓柱／圓孔的介電常數(shù)
die2:REAL(8)，基體的介電常數(shù)
k:INTEGER,Bragg波矢的階數(shù)
l:INTEGER,Bragg波矢的階數(shù)
n:INTEGER，n=1時(shí)計(jì)算介電函數(shù)的Fourier系數(shù)；n=2時(shí)計(jì)算介電函數(shù)倒數(shù)的Fourier系數(shù)
輸出參數(shù)：
ftrans:REAL(8)，周期函數(shù)的Fourier系數(shù)