前情提要
本文Ince高斯光束的模拟是在基模高斯光束的基础上进行的,其中大部分基本参量与基模高斯光束的一致,如有需要可以参考我的另一篇文章《基模高斯光束的MATLAB仿真模拟》
基模高斯光束的MATLAB仿真模拟http://t.csdnimg.cn/ON1S3
公式部分
Ince高斯光束(IG)是高斯光束在椭圆-双曲坐标系(简称椭圆坐标系)下的一种模式,其复振幅的分离变量带入傍轴波动方程(PWE)后满足Ince方程,且可以通过调节椭圆率实现向厄密高斯光束和拉盖尔高斯光束模的转换。
首先相比于基模高斯光束,IG的基本参量多了以下三项:
- 阶数------------------u
- 阶数------------------v
- 椭圆率---------------e(0~∞)
由于Ince方程的本征解分为偶Ince多项式和奇Ince多项式,Ince高斯光束自然有偶模式和奇模式,电场表达式如下所示:
其中,和分别为偶模式和奇模式的归一化系数,和分别为(u,v)阶的偶Ince多项式和奇Ince多项式,并且对于偶恩斯多项式要求0 <=v<= u,对于奇恩斯多项式要求1<=v<= u,u和v具有相同的奇偶性,即(-1)^(v-u)=1 。
Ince高斯光束的电场表达式(E_Ince)与基模高斯光束的电场表达式(E_Gauss)的关系如下:
E_Ince=C_IG*C_vu(etha,e)*C_vu(i*xhi,e)*E_Gauss*exp((-u+1)*Phi))/w0 偶模式
E_Ince=S_IG*S_vu(etha,e)*S_vu(i*xhi, e)*E_Gauss*exp((-u+1)*Phi)/w0 奇模式
依据以上公式便可以编写模拟程序了。
代码部分
首先确定基本参量并进行简单计算:
%% 基本参量
parity=input('Ince高斯光束模式(奇1或偶0):');%选择偶/奇模式
lambda=input('入射光波波长(mm)='); %光束的输入波长mm
z=input('传播距离='); %光束传播距离
v=input('阶数v=');
u=input('阶数u='); %Ince多项式的阶数
e=input('椭圆率e='); %椭圆坐标系的椭圆率
k=2*pi/lambda; %传播常数,又称波数
w0=1.5*lambda; %高斯光束的中心束腰半径
ZR=k*w0^2/2; %瑞利长度或共焦参数
wz=w0*sqrt(1+(z/ZR)^2); %光斑范围半径,传播到z处的束宽
Rz=z+ZR*ZR/z; %高斯球面波的波前曲率半径
Phi=atan(z/ZR); %古伊相位
由于Ince多项式的条件:偶恩斯多项式要求0 <=v<= u,对于奇恩斯多项式要求1<=v<= u ,(-1)^(v-u)=1。需要对输入的参数进行校验:
%% 输入校验
if parity==0
if (v<0)||(v>u); error('ERROR:偶恩斯多项式: m,p 应满足 0<=m<=p'); end
else
if (v<1)||(v>u); error('ERROR:奇恩斯多项式: m,p 应满足1<=m<=p'); end
end
if (-1)^(v-u)~=1; error('ERROR: m,p 的奇偶性应相同,即(-1)^(m-p)=1');end
然后建立坐标系(椭圆坐标系):
%% 坐标参数
L=5*lambda; %坐标尺寸,[-L,L]*[-L,L]
N=1001; %取样点数
f0 = sqrt(e/2)*w0; %z为0处的椭圆焦距,e为椭圆率
f = f0*wz/w0; %z处的椭圆焦距
[xhi,etha,x1,y1]=mesh_elliptic(f0,L,N);%椭圆坐标系
r = sqrt(x1.^2 + y1.^2); %极坐标系
其中椭圆坐标系函数mesh_elliptic(f,L,N)如下:
%% 笛卡尔坐标系转椭圆坐标系
function [xi,eta,x1,y1]=mesh_elliptic(f,L,N)
[x1,y1]=meshgrid(linspace(-L,L,N));
xi=zeros(N); eta=zeros(N); % 初始化
% 计算第一象限
en = acosh((x1(1:(N+1)/2,(N+1)/2:N)+1i*y1(1:(N+1)/2,(N+1)/2:N))/f);
ee = real(en);
nn = imag(en);
nn = nn + (nn<0)*2*pi;
xi(1:(N+1)/2,(N+1)/2:N)=ee;
eta(1:(N+1)/2,(N+1)/2:N)=nn;
% 计算其他象限
xi(1:(N+1)/2 , 1:(N-1)/2)=fliplr(xi(1:(N+1)/2,(N+3)/2:N));
xi((N+3)/2:N,1:N)=flipud(xi(1:(N-1)/2,1:N));
eta(1:(N+1)/2,1:(N-1)/2)=pi-fliplr(eta(1:(N+1)/2,(N+3)/2:N));
eta((N+3)/2:N,1:N)=pi+rot90((eta(1:(N-1)/2,1:N)),2);
return
接下来根据公式编写程序实现IG的表达:
%% Ince高斯光束的表达式
E_Gauss=(w0./wz).*exp(-(r.^2./wz^2)).*exp(1i.*(k.*z+k.*r.^2./(2*Rz)-Phi));%基模高斯光束的完整的电场表达式
if z==0
if parity==0 %奇偶校验
E_Ince=CInce(u,v,e,etha).*CInce(u,v,e,1i*xhi).*exp(-(r/w0).^2); % 偶
else
E_Ince=SInce(u,v,e,etha).*SInce(u,v,e,1i*xhi).*exp(-(r/w0).^2); % 奇
end
else
if parity==0
E_Ince=(CInce(u,v,e,etha).*CInce(u,v,e,1i*xhi)).*E_Gauss.*exp((-u+1).*Phi); % 偶
else
E_Ince=(SInce(u,v,e,etha).*SInce(u,v,e,1i*xhi)).*E_Gauss.*exp((-u+1).*Phi); % 奇
end
end
if parity == 0
if mod(u,2)==0
[C0,~,coef]=CInce(u,v,e,0);
[Cp,~,~]=CInce(u,v,e,pi/2);
Norm = (-1)^(v/2)*sqrt(2)*gamma(u/2+1)*coef(1) *sqrt(2/pi)/w0/C0/Cp; %计算归一化系数
else
[C0,~,coef]=CInce(u,v,e,0);
[~,~,~,DCp]=CInce(u,v,e,pi/2);
Norm = (-1)^((v+1)/2) * gamma((u+1)/2+1) * sqrt(4*e/pi) * coef(1) / w0 / C0 / DCp; %计算归一化系数
end
else
if mod(u,2)==0
[~,~,coef,dS0]=SInce(u,v,e,0);
[~,~,~,dSp]=SInce(u,v,e,pi/2);
Norm = (-1)^(v/2)*sqrt(2)*e*gamma((u+2)/2+1)*coef(1) *sqrt(2/pi)/w0/ dS0 / dSp; %计算归一化系数
else
[Sp,~,coef,~]=SInce(u,v,e,pi/2);
[~,~,~,dS0]=SInce(u,v,e,0);
Norm = (-1)^((v-1)/2) * gamma((u+1)/2+1) * sqrt(4*e/pi) * coef(1) / w0 / Sp / dS0;%计算归一化系数
end
end
E_Ince = E_Ince*Norm; %乘以归一化系数
I_Ince=E_Ince.*conj(E_Ince); %计算Ince高斯光束的光强
normI_Ince=I_Ince/max(max(I_Ince)); %归一化光强
其中,(u,v)阶的偶Ince多项式和奇Ince多项式函数如下:
%% 偶Ince多项式
function [IP,eta,coef,dIP]=CInce(u,v,q,z)
%% 输入校验
if (v<0)||(v>u); error('ERROR:"v"范围错误, 0<=v<=u'); end
if (-1)^(v-u)~=1; error('ERROR: (u,v)必须有相同的奇偶性, 即(-1)^(v-u)=1'); end
[largo,ancho]=size(z); % 转向量形式
z=transpose(z(:));
normalization=1;
%% 计算系数
if mod(u,2)==0
%===== u为偶数=======
j=u/2; N=j+1; n=v/2+1;
% 矩阵
M=diag(q*(j+(1:N-1)),1) + diag([2*q*j,q*(j-(1:N-2))],-1) + diag([0,4*((0:N-2)+1).^2]);
if u==0; M=0; end
% 特征值和特征向量
[A,ets]=eig(M);
ets=diag(ets);
[ets,index]=sort(ets);
A=A(:,index);
% 归一化
if normalization==0
N2=2*A(1,n).^2+sum(A(2:N,n).^2);
NS=sign(sum(A(:,n)));
A=A/sqrt(N2)*NS;
else
mv=(2:2:u).';
N2=sqrt(A(1,n)^2*2*gamma(u/2+1)^2+sum((sqrt(gamma((u+mv)/2+1).*gamma((u-mv)/2+1) ).*A(2:u/2+1,n)).^2 ));
NS=sign(sum(A(:,n)));
A=A/N2*NS;
end
% ince多项式
r=0:N-1;
[R,X]=meshgrid(r,z);
IP=cos(2*X.*R)*A(:,n);
dIP=-2*R.*sin(2*X.*R)*A(:,n);
eta=ets(n);
else
%===== u为偶数=======
j=(u-1)/2; N=j+1; n=(v+1)/2;
% 矩阵
M=diag(q/2*(u+(2*(0:N-2)+3)),1)+diag(q/2*(u-(2*(1:N-1)-1)),-1) + diag([q/2+u*q/2+1,(2*(1:N-1)+1).^2]);
% 特征值和特征向量
[A,ets]=eig(M);
ets=diag(ets);
[ets,index]=sort(ets);
A=A(:,index);
% 归一化
if normalization==0
N2=sum(A(:,n).^2);
NS=sign(sum(A(:,n)));
A=A/sqrt(N2)*NS;
else
mv=(1:2:u).';
N2=sqrt(sum( ( sqrt(gamma((u+mv)/2+1).*gamma((u-mv)/2+1) ).*A(:,n)).^2 ));
NS=sign(sum(A(:,n)));
A=A/N2*NS;
end
% ince多项式
r=2*(0:N-1)+1;
[R,X]=meshgrid(r,z);
IP=cos(X.*R)*A(:,n);
dIP=-R.*sin(X.*R)*A(:,n);
eta=ets(n);
end
coef=A(:,n);
IP=reshape(IP,[largo,ancho]); % 重构输出结果数组,满足原始结构
dIP=reshape(dIP,[largo,ancho]);
return
%% 奇恩斯多项式
function [IP,eta,coef,dIP]=SInce(u,v,q,z)
%% 输入校验
if (v<1)||(v>u); error('ERROR:"v"范围错误, 1<=v<=u'); end
if (-1)^(v-u)~=1; error('ERROR: (u,v)必须有相同的奇偶性,即(-1)^(v-u)=1'); end
[largo,ancho]=size(z); % 转向量形式
z=transpose(z(:));
normalization=1;
%% 计算系数
if mod(u,2)==0
%======u为偶数=======
j=u/2; N=j+1; n=v/2;
% 矩阵
M=diag(q*(j+(2:N-1)),1)+diag(q*(j-(1:N-2)),-1) + diag(4*((0:N-2)+1).^2);
% 特征值和特征向量
[A,ets]=eig(M);
ets=diag(ets);
[ets,index]=sort(ets);
A=A(:,index);
% 归一化
r=1:N-1;
if normalization==0
N2=sum(A(:,n).^2);
NS=sign(sum(r.*transpose(A(:,n))));
A=A/sqrt(N2)*NS;
else
mv=(2:2:u).';
N2=sqrt(sum((sqrt(gamma((u+mv)/2+1).*gamma((u-mv)/2+1) ).*A(:,n)).^2 ));
NS=sign(sum(r.*A(:,n)'));
A=A/N2*NS;
end
% ince多项式
[R,X]=meshgrid(r,z);
IP=sin(2*X.*R)*A(:,n);
dIP=2*R.*cos(2*X.*R)*A(:,n);
eta=ets(n);
else
%========u为奇数=======
j=(u-1)/2; N=j+1; n=(v+1)/2;
% 矩阵
M=diag(q/2*(u+(2*(0:N-2)+3)),1)+diag(q/2*(u-(2*(1:N-1)-1)),-1) + diag([-q/2-u*q/2+1,(2*(1:N-1)+1).^2]);
% 特征值和特征向量
[A,ets]=eig(M);
ets=diag(ets);
[ets,index]=sort(ets);
A=A(:,index);
% 归一化
r=2*(0:N-1)+1;
if normalization==0
N2=sum(A(:,n).^2);
NS=sign(sum(r.*transpose(A(:,n))));
A=A/sqrt(N2)*NS;
else
mv=(1:2:u).';
N2=sqrt(sum( ( sqrt(gamma((u+mv)/2+1).*gamma((u-mv)/2+1) ).*A(:,n)).^2 ));
NS=sign(sum(r.*A(:,n)'));
A=A/N2*NS;
end
% ince多项式
[R,X]=meshgrid(r,z);
IP=sin(X.*R)*A(:,n);
dIP=R.*cos(X.*R)*A(:,n);
eta=ets(n);
end
coef=A(:,n);
IP=reshape(IP,[largo,ancho]); % 重构输出结果
dIP=reshape(dIP,[largo,ancho]);
return
最后,绘制Ince高斯光束的图像:
figure()
pcolor(x1,y1,normI_Ince); %二维图像
set(gca,'fontname','times new roman','fontsize',16);
xlabel('\itx(\it\lambda)','fontname','times new roman','fontsize',24);
ylabel('\ity(\it\lambda)','fontname','times new roman','fontsize',24);
shading interp; colormap hot; colorbar;
figure()
pcolor(x1,y1,angle(E_Ince)); %位相分布
set(gca,'fontname','times new roman','fontsize',16);
xlabel('\itx(\it\lambda)','fontname','times new roman','fontsize',24);
ylabel('\ity(\it\lambda)','fontname','times new roman','fontsize',24);
clim([-pi,pi]); %颜色映射范围
shading interp; colormap hot; colorbar;
figure()
surf(x1,y1,normI_Ince); %三维图像
set(gca,'fontname','times new roman','fontsize',16);
xlabel('\itx(\it\lambda)','fontname','times new roman','fontsize',24);
ylabel('\ity(\it\lambda)','fontname','times new roman','fontsize',24);
clim([0,1]);
shading interp; colormap hot; colorbar;
theAxes=axis; %保存坐标变量
值得注意的是,绘制位相分布图的时候不能将颜色映射范围(clim)归一化,应与相位范围一致,如设为[-pi,pi]。如果输出的位相图只有最大值和最小值对应的颜色,很可能就是这里范围设置过小。
结果部分
Ince高斯光束的相位分布是阶跃式的:
当椭圆率趋近于无穷时,Ince高斯光束将变换成厄密高斯光束的模式,如下图所示:
当椭圆率趋近于零时,Ince高斯光束将倾向于拉盖尔高斯光束的模式,不同之处IG是相位是阶跃的,LG是渐变的,如下图所示:
以上实现了在matlab中对Ince高斯光束的模拟。
参考文献
Miguel A. Bandres (2023). Ince Gaussian Beam (https://www.mathworks.com/matlabcentral/fileexchange/46222-ince-gaussian-beam), MATLAB Central File Exchange.