对于相位匹配来说，计算相位匹配角是一个繁琐的过程。
以下是通过matlab来辅助计算单轴晶体的相位匹配角。
可以简单计算单轴晶体的一类相位匹配以及二类相位匹配。

懒得检查代码了，各位大神欢迎指出错误。

%负单轴晶体一类相位匹配
%首先采用KDP晶体举例，相位匹配条件为
%   ne(532)(θ)=no(1064)
%θ的角度为波矢K与Z轴的夹角
%φ的角度为波矢K在xoy面上的投影与X轴的夹角
%单轴晶体主轴可以垂直于xy平面上任意选择，为方便起见，选择y轴方向使给
%定的波法线方向位于yoz平面，Z轴为光轴方向。
%使基频光为O光，倍频光为e光。
%要求基频光与倍频光折射率相等，即无色散。clc
clear alllamda1=1064/10^3;         %sellmeier公式的单位是微米
lamda2=532/10^3;
Ano=2.259276;
Bno=0.01008956;
Cno=0.012942625;
Dno=13.00522;
Ane=2.132668;
Bne=0.008637494;
Cne=0.012281043;
Dne=3.227992;%1064nm波长的o光  e光折射率
no=sqrt(Ano+Bno/(lamda1^2-Cno)+(Dno*lamda1^2)/(lamda1^2-400));
ne=sqrt(Ane+Bne/(lamda1^2-Cne)+(Dne*lamda1^2)/(lamda1^2-400));
fprintf('基频光1064nm波长的折射率 no')
disp(no)
fprintf('基频光1064nm波长的折射率 ne')
disp(ne)% 532nm波长的o光  e光折射率
n2o=sqrt(Ano+Bno/(lamda2^2-Cno)+(Dno*lamda2^2)/(lamda2^2-400));
n2e=sqrt(Ane+Bne/(lamda2^2-Cne)+(Dne*lamda2^2)/(lamda2^2-400));
fprintf('倍频光532nm波长的折射率  no')
disp(n2o)
fprintf('倍频光532nm波长的折射率  ne')
disp(n2e)% 求theta的角度         (第一种方法：菲涅尔公式 通过K值求)
b=asin(sqrt((no^(-2)-n2o^(-2))/(n2e^(-2)-n2o^(-2))));
%求得结果为弧度
theta=b*180/pi;
%将弧度转换为角度
%角度转弧度: π/180×角度 ;弧度变角度:180/π×弧度
fprintf('方法一：Ⅰ类相位匹配的θ角为')
disp(theta)
%输出的为一类相位匹配的θ角（角度），此时φ为45度。% 求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('方法二：Ⅰ类相位匹配的θ角为')
disp(d)%负单轴晶体二类相位匹配
%相位匹配条件：     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;
y=solve(f1);
%disp(y)                  %输出所有解，包括实数和复数
yy=double(y);             %将符号型转换为数值型
%disp(yy)                 %输出数值型的解  
h=asin(sqrt(yy));
g=h*180/pi;               %转换为角度显示
% for i=1:length(g)               %另外一种选取实数的方法
%     if isreal(g(i))
%         m(length(g)+1)=g(i);
%     end
% end
% disp(m)
A=g==real(g);             %选取实数，如果是实数为1，不是就为0
fprintf('Ⅱ类相位匹配的θ角为        ')
disp(g(A))                %输出为1的实数
%输出的为二类相位匹配的θ角（角度），此时φ为0度。%求KDP晶体一类和二类相位匹配的有效非线性光学系数deff
d14=0.46;                 %d14=d36=0.46 (pm/V)
d36=0.46;
fai1=45;                  %一类相位匹配φ=45度
fai2=0;                   %二类相位匹配φ=0度
deff1=d14*sind(theta)*sind(2*fai1);
deff2=d14*sind(2*g(A))*cosd(2*fai2);
fprintf('Ⅰ类匹配的有效非线性系数deff为')
disp(deff1)
fprintf('Ⅱ类匹配的有效非线性系数deff为')
disp(deff2)%求离散角（离散效应）
%在晶体光学中，把光波波法线方向与光线方向的夹角称为离散角
%对于单轴晶体，O光的离散角恒等于零。e光的离散角为%倍频过程只有532产生离散角
B1=(1/2)*n2e^2*(n2o^-2-n2e^-2)*sind(2*theta);
B2=(1/2)*n2e^2*(n2o^-2-n2e^-2)*sind(2*g(A));
C1=B1*180/pi;
C2=B2*180/pi;
fprintf('Ⅰ类相位匹配的离散角')
disp(C1)
fprintf('Ⅱ类相位匹配的离散角')
disp(C2)
%对于正单轴晶体，ne>no,离散角C>0,e光的光线较其波法线靠近光轴。
%对于负单轴晶体，ne<no,离散角C<0,e光的光线较其波法线远离光轴。