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

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





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







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






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






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




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






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






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





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





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






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






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