可以把共线相位匹配3.1这里讲解清楚，同时给我一个复现图2(b)的代码吗，matlab或者mathematica

### 共线相位匹配（3.1）

#### 3.1.1 双o类型

**定义**：双o类型是指在三波相互作用中，有两个光波为o光（寻常光），另一个为e光（非常光）。

**计算公式**：
对于负单轴晶体（例如ooe类型），可以根据以下公式计算相位匹配角：

\[ n(\theta) = n_o \sqrt{1 + \tan^2 \theta} / \sqrt{1 + \left(\frac{n_o}{n_e}\right)^2 \tan^2 \theta} \]

其中：
- \( n_o \) 和 \( n_e \) 分别是o光和e光的折射率。
- \( \theta \) 是e光与光轴的夹角。

根据相位匹配条件 \( n_1 \lambda_1 + n_2 \lambda_2 = n_3 \lambda_3 \)，可以得到：

\[ A + B = C \left(1 + \tan^2 \theta\right) / \left(1 + \left(\frac{F}{E}\right)^2 \tan^2 \theta\right) \]

其中：
- \( A = n_{o1} / \lambda_1 \)
- \( B = n_{o2} / \lambda_2 \)
- \( C = n_{o3} / \lambda_3 \)
- \( D = n_{e1} / \lambda_1 \)
- \( E = n_{e2} / \lambda_2 \)
- \( F = n_{e3} / \lambda_3 \)

求解上述方程，得到相位匹配角 \( \theta \)：

\[ \tan^2 \theta = \frac{W \left(1 - U/V\right)}{1 - W} \]

其中：
- \( U = \left(\frac{A + B}{C}\right)^2 \)
- \( V = \left(\frac{A + B}{F}\right)^2 \)
- \( W = 1 \)

#### 3.1.2 双e类型

**定义**：双e类型是指在三波相互作用中，有两个光波为e光，另一个为o光。

**计算公式**：
对于负单轴晶体（例如oee类型），可以根据以下公式计算相位匹配角：

\[ n_1 = n_{o1}, \quad n_2 = n_{o2} \sqrt{1 + \tan^2 \theta} / \sqrt{1 + \left(\frac{n_{o2}}{n_{e2}}\right)^2 \tan^2 \theta}, \quad n_3 = n_{o3} \sqrt{1 + \tan^2 \theta} / \sqrt{1 + \left(\frac{n_{o3}}{n_{e3}}\right)^2 \tan^2 \theta} \]

根据相位匹配条件 \( n_1 \lambda_1 + n_2 \lambda_2 = n_3 \lambda_3 \)，可以得到：

\[ A + B \left(1 + \tan^2 \theta\right) / \left(1 + \left(\frac{E}{D}\right)^2 \tan^2 \theta\right) = C \left(1 + \tan^2 \theta\right) / \left(1 + \left(\frac{F}{E}\right)^2 \tan^2 \theta\right) \]

求解上述方程，得到相位匹配角 \( \theta \)：

\[ \tan^2 \theta = \frac{W \left(1 - U/V\right)}{1 - W} \]

其中：
- \( U = \left(\frac{C}{A + B}\right)^2 \)
- \( V = \left(\frac{F}{A + E}\right)^2 \)
- \( W = \left(\frac{F}{C}\right)^2 \)

### 复现图2(b)的MATLAB代码

以下是用于复现图2(b)的MATLAB代码示例：

% 定义波长和折射率
lambda1 = 1588.1e-9; % 波长1
lambda2 = 800e-9;    % 波长2
lambda3 = 532e-9;    % 波长3

% 折射率数据（示例数据）
no1 = 1.67; % o光折射率1
no2 = 1.65; % o光折射率2
no3 = 1.64; % o光折射率3
ne1 = 1.69; % e光折射率1
ne2 = 1.66; % e光折射率2
ne3 = 1.63; % e光折射率3

% 计算等效折射率
A = no1 / lambda1;
B = no2 / lambda2;
C = no3 / lambda3;
D = ne1 / lambda1;
E = ne2 / lambda2;
F = ne3 / lambda3;

% 计算近似解的误差
eccentricity_values = linspace(0, 0.5, 100);
error_values = zeros(size(eccentricity_values));

for i = 1:length(eccentricity_values)
    ecc = eccentricity_values(i);
    % 修改折射率椭圆的长半轴
    ne2_mod = sqrt(ne2^2 * (1 - ecc^2));
    F_mod = ne3 / lambda3;
    
    % 计算相位匹配角
    U = (C / (A + B))^2;
    V = (F_mod / (A + E))^2;
    W = (F_mod / C)^2;
    tan_theta_squared = (W * (1 - U/V)) / (1 - W);
    theta = atan(sqrt(tan_theta_squared)) * 180 / pi;
    
    % 数值解
    f = @(theta) A + B * (1 + tan(theta)^2) / (1 + (E/D)^2 * tan(theta)^2) - C * (1 + tan(theta)^2) / (1 + (F_mod/E)^2 * tan(theta)^2);
    theta_numerical = fzero(f, [0, pi/2]) * 180 / pi;
    
    % 计算误差
    error_values(i) = abs(theta - theta_numerical);
end

% 绘制误差与离心率的关系
figure;
plot(eccentricity_values, error_values, 'LineWidth', 2);
xlabel('Refractive Index Ellipse Eccentricity');
ylabel('Approximate Solution Error (degrees)');
title('Relationship Between Approximate Solution Error and Refractive Index Ellipse Eccentricity');
grid on;

这段代码首先定义了波长和折射率，然后计算了等效折射率。接着，它通过改变折射率椭圆的离心率来计算近似解和数值解的误差，并绘制了误差与离心率的关系图。你可以根据实际的折射率数据调整 `no1`, `no2`, `no3`, `ne1`, `ne2`, `ne3` 的值。