MATLAB：BBO晶体的相位匹配角计算

MATLAB：BBO晶体的相位匹配角计算

通过matlab来辅助计算BBO晶体的相位匹配角以及deff系数。
本代码只能计算SHG。
实际使用的角度与计算角度不一致的主要因素是晶体温度导致。

%BBO单轴晶体
%实际角度还需按照使用时晶体温度来优化。
%因为温度改变了晶体的折射率。


clcsell
clear all
close all

lambda1=1064/10^3;        %在此处输入基频光波长，单位微米
lambda2=532/10^3;         %在此处输入倍频光波长，单位微米
                          %sellmeier公式的单位是微米
Ano=2.7359;               %ABCD分别对应sellmeier公式的4个系数
Bno=0.01878;
Cno=0.01822;
Dno=0.01354 ;
Ane=2.3753;
Bne=0.01224;
Cne=0.01667;
Dne=0.01516;
%基频光波长的o光  e光折射率
no=sqrt(Ano+Bno/(lambda1^2-Cno)-(Dno*lambda1^2));
ne=sqrt(Ane+Bne/(lambda1^2-Cne)-(Dne*lambda1^2));
fprintf('基频光波长的折射率  no=%6.4f',no);
fprintf('       ne=%6.4f\n',ne);
% 倍频光波长的o光  e光折射率
n2o=sqrt(Ano+Bno/(lambda2^2-Cno)+(Dno*lambda2^2));
n2e=sqrt(Ane+Bne/(lambda2^2-Cne)+(Dne*lambda2^2));
fprintf('倍频光波长的折射率  no=%6.4f',n2o);
fprintf('       ne=%6.4f\n',n2e);

% 求Ⅰ类相位匹配theta角        (第一种方法：菲涅尔公式 通过K值求)
a=asin(sqrt((no^(-2)-n2o^(-2))/(n2e^(-2)-n2o^(-2))));%求得结果为弧度
b=a*180/pi;%将弧度转换为角度
%角度转弧度: π/180×角度 ;弧度变角度:180/π×弧度
fprintf('Ⅰ类相位匹配的θ角为 θ=%6.2f\n',b);

% % % Ⅰ类相位匹配theta角          (第二种方法：相位匹配条件化简)
% % c=asin(sqrt(((n2o^2)*(n2e^2)-(n2e^2)*(no^2))/((n2o^2)*(no^2)-(n2e^2)*(no^2))));
% % d=c*180/pi;
% % fprintf('Ⅰ类相位匹配的θ角为 θ=%6.2f\n',d)

%BBO晶体Ⅱ类相位匹配
%相位匹配条件：     1/2[n_e^ω (θ_m )+n_o^ω]= n_e^2ω (θ_m )
%设，sinθ^2为x    cosθ^2为1-x
syms x
f1=2/sqrt(((1-x)/n2o^2)+(x/n2e^2))-1/sqrt(((1-x)/no^2)+x/ne^2)-no;
f2=solve(f1);              %输出所有解，包括实数和复数
f3=double(f2);             %将符号型转换为数值型
f4=asin(sqrt(f3));         %将数值型进行asin求值
f5=f4*180/pi;              %转换为角度显示
A=f5==real(f5);             %选取实数，如果是实数为1，不是就为0
fprintf('Ⅱ类相位匹配的θ角为 θ=%6.2f\n ',f5(A)) %输出为1的实数


%BBO φ角计算，deff系数计算
theta=b;                   %引用Ⅰ类相位匹配theta角度
theta2=f5(A);              %引用Ⅱ类相位匹配theta角度
fai1 = 0:0.001:90; 
d36=0.44;
d11=5.8*d36;     %d36(KDP)
d31=0.05*d11;
d22=0.05*d11;    %此处应为d22＜0.05*d11，为求近似值选择使用＝
deff1=d31*sind(theta)+(d11*cosd(3*fai1)-d22*sind(3*fai1)*cosd(theta));
deff2=(d11*sind(3*fai1)+d22*cosd(3*fai1))*cosd(theta2)*cosd(theta2);
plot(fai1,deff1,'red');hold on
plot(fai1,deff2,'green'); grid on; 
xlabel( 'φ角' )
ylabel( '非线性系数' )
legend('Ⅰ类相位匹配deff','Ⅱ类相位匹配deff','Location', 'northeast' )
title( '非线性系数随φ角的变化图' )


[deff1max, fai1max] = max(deff1);
[deff2max, fai2max] = max(deff2);
fai1=fai1max*0.001;
fai2=fai2max*0.001;
fprintf('Ⅰ类相位匹配最大非线性系数  deff=%0.3f\n',deff1max')
fprintf(' Ⅰ类相位匹配               θ=%0.2f',theta);
fprintf('   φ=%0.2f\n',fai1);
fprintf(' Ⅱ类相位匹配最大非线性系数  deff=%0.3f\n',deff2max')
fprintf(' Ⅱ类相位匹配               θ=%0.2f',theta2);
fprintf('   φ=%0.2f\n',fai2);

如果遇到该情况：

表示该波长不满足Ⅱ类相位匹配条件。
可将代码中关于计算Ⅱ类相位匹配的屏蔽掉。%%屏蔽。

附：
物理光学特性：

应用参考：