04-bpm-标量光束传播方法

光束传播方法(BeamPropagationMethod,BPM)丁卫强哈尔滨工业大学物理系•各向异性介质中电磁波的本征模求解2上节课内容回顾xxxyxzyxyyyzzxzyzz~exp~expEitixyztkr22=0kkEkEE2222222222222220xxxyxzxyxyyyzyzxzyzzzEEE波动方程41ˆexpzApixyztE41ˆexpAqixyztHˆˆˆ/qckpˆˆˆˆzkxyz•传输矩阵的建立3上节课内容回顾,xyE,xyH1nzz在交界面上连续，可得如下四个方程其中：z0z1z2zNz1nz4上节课内容回顾偏振相关传播矩阵，位相相关波导结构脊型波导，掩埋式波导，平面波导5沿传播方向，折射率均匀/或者变化很小；但是传播方向的距离很长！光纤、介质波导、集成光学结构……通常用光束传播方法比较合适在这类结构中，光波的传输用亥母霍兹方程精确描述，并且沿传播方向的电磁场都可以分解为慢变振幅项与位相振荡项。6光束传播方法横向：快速傅立叶变换（FFT-BPM），有限差分（FD-BPM），以及有限元（FEM-BPM）传播方向上：将空间上相邻的两层之间的场用一个“传播算符”联系起来，这样从第一层开始递推，可以得到整个空间的场。…z0z…zzM传递矩阵联系7二维FD-BPM算法00,=rjjEHHE时谐场的Maxwell方程出发假设，沿y方向结构是均匀的，即，则：0y6个电磁场分量自然分解为2组：一组只包含（Ey，Hx，Hz）一组只包含（Ex，Ez，Hy）TE模TM模tjTE模的Ey控制方程：8TE模01yxEHjz得到：2220220yyryEEkEzxTE模的Hx控制方程：将方程3对z求导：2202yxzrEHHjzzxz/0rz另外，由于0H2220220xxrxHHkHzx01yzEHjx22000k其中0yxzHHHxzyxzHHzx9TM模的控制方程Ex的控制方程：2202yxzHEEjzzxz2000rxrxjjEkE0rDE0rxrzEExz1zrxrEEzx220210xrxrxrEEkEzxx最终得到：电场高斯定理：/0rz10TM模的控制方程：Hy01yxrHEjz2200011yyyyrrHHHkHzzxx再次利用假设：结构沿z方向变化非常缓慢2211yyrrHHzzz220210yyrryrHHkHzxx最后得到：01yzrHEjx/0rz11TE模的有限差分离散2220220yyryEEkEzx假设，沿传播方向Ey表示为慢变振幅和快速振荡的位相两部分：,,,,jzyExyzxyze其中0effnk22222exp2expexpyEjzjjzjzzzz2222expyEjzxx22220222reffjknzzx带入Ey控制方程得：如果令，则得到近轴近似22/0z12有限差分离散222022reffjknzx,xpxzlz令：,lpxz,lrrxzp21/21/211222pppppxxxx右边第一项：右边第二项：222200reffreffpknkpn带入方程得到：11220222ppppreffpjkpnzx22101wpxreffpepkpn其中系数定义为：21wx22xx21ex连续空间中，只取有限个离散的点13左边对z差分离散：122llppjjzz注意：离散中心在l+1/22210111112212121110011222wpxreffpepllllwpwpllllllllllxreffpxreffpepepkpnkpnkpn右边对z的差分中心也应该在l+1/2，1112121111012210144lllllllwpxreffpeplllllllwpxreffpepjkpnzjkpnz将l+1和l项分别整理到方程左右两侧，得：11/22lllxxx利用平均值近似14TM模的差分离散220210yyrryrHHkHzxx得到近轴近似方程：22012rreffrjknzxx,,,,,llprrxpxzlzxzxzp利用记号：111/21/11211pppprrrrrpxxxpxpx右边第一项为：P-1/2P+1/2同样的方法，令且22/0yHz,,,,,jzyHxyzxyze15111/2,1/222rrrrrrpppppp1122222111122111rrpprrrrrrprrrrppppppxxppppppx222200reffreffpknkpn221012wpxreffpepjkpnz中值近似：1112211rpppprrrrpxppxppx111/21/11211pppprrrrrpxxxpxpx11wpxpep右边第二项22012rreffrjknzxxTM模的差分离散左边差分为：122llppjjzz差分中心为l+1/2同样的方法，利用平均值近似，将l+1/2的值用l和l+1代替，同时将l+1和l项分别整理到方程两侧，得：1112121111012210144lllllllwpxreffpeplllllllwpxreffpepjkpnzjkpnz222211,,11rrwerrrrppppppxxxew其中系数分别为：现在要解决的问题就是利用z处的值求得处的未知值17BPM算法的线性方程组lzz1l1112121111012210144lllllllwpxreffpeplllllllwpxreffpepjkpnzjkpnz例如：其中：也就是：取遍所有的p值，得到p个方程。所以BPM算法归结为求解一个三元一次方程组18边界条件边界的反射干扰计算区域的场1234…M-2M-1MM+100x1x2x3x1MxMx1Mx,expxxzAzjkx对于左边界，沿x方向传播的波可以写为：则0，1，2三点处的值分别为：计算区域边界格点计算时，信息不够1,,1ppp19透明边界条件（TBC）2110expxjkx2110xxxxx211expxjkx令：01expxjkx同时：则：11lnxkjx同样的道理，对于右边边界，可以用同样的方法处理：1expMMxjkx其中：1lnxMkjx1MMM可得：对于左边界(p=1)：由于左边界的值φ1受φ0的影响，即：20边界条件处理011112111/Lll1110121111lllABCD111111lllLwLA1112'111llBCD11111lllpppApBpCpDp中，令p=1在方程21右边边界(p=M)：111/MMRMllMM右边的方程简化为：111'llMMAMBMDM最终得到矩阵方程：算法程序流程图设置结构参数,rxz确定工作频率00,确定离散步长,xz确定入射波源0z对传播方向循环l=1:L传播方向循环结束对横向格点循环p=2:P-1构造出系数矩阵求解方程得到l+1处的值对p=1和p=P应用边界条件做图显示，及其他后处理程序实例closeall;clear;clc;n0=1;dz=0.4e-6;dx=0.2e-6;x=(-50:50)*dx;lam0=1e-6;neff=n0;k0=2*pi/lam0;beta=neff*k0;epsr=1;23%DefinethematrixesLx=numel(x);Lz=200;EE=zeros(Lx,Lz);M=zeros(Lx,Lx);D0=zeros(Lx,1);%Definethegridsize;alphaw=1/dx^2;alphax=-2/dx^2;alphae=1/dx^2;%Boundaryvalueatz=0phi=exp(-x.^2/(4*lam0).^2);forL=1:Lz%对传播方向的格点循环p=1;gammaL=phi(1)/phi(2);B1=-alphaw*gammaL+(-alphax+4j*beta/dz-k0^2*(epsr-neff^2));C1=-alphae;D1=alphaw*gammaL*phi(1)+...(alphax+4j*beta/dz+k0^2*(epsr-neff^2))*phi(1)+...alphae*phi(2);M(p,p:(p+1))=[B1,C1];D0(p)=D1;forp=2:Lx-1A=-alphaw;B=-alphax+4j*beta/dz-k0^2*(epsr-neff^2);C=-alphae;D=alphaw*phi(p-1)+...(alphax+4j*beta/dz+k0^2*(epsr-n