帮我复现一下图1，用mathematica

要使用Mathematica复现图1中的数据，首先需要理解图1的数据生成原理。根据文中提供的信息，图1展示了BBO晶体在Ⅰ型非共线相位匹配下的光学参量放大特性，特别是信号光的增益谱和相位谱随着抽运光强度和非共线夹角的变化。

### 步骤1: 定义物理常数和公式

首先定义一些必要的物理常数和公式：

(* 物理常数 *)
\[Lambda]p = 532*10^-9; (* 抽运光波长, 单位: 米 *)
\[Lambda]s0 = 800*10^-9; (* 信号光中心波长, 单位: 米 *)
L = 13.5*10^-3; (* 晶体长度, 单位: 米 *)
Ip = 350*10^6; (* 抽运光强度, 单位: W/cm^2 *)

(* 折射率公式 *)
no[\[Lambda]_] := Sqrt[2.7359 + 0.01878/(\[Lambda]^2 - 0.01822) - 0.01354*\[Lambda]^2];
ne[\[Lambda]_] := Sqrt[2.3753 + 0.01224/(\[Lambda]^2 - 0.01667) - 0.01516*\[Lambda]^2];

(* 计算相位失配 *)
\[Delta]k[\[Lambda]s_, \[Alpha]_] := 
  (2 \[Pi]/\[Lambda]i[\[Lambda]s]) (ni[\[Lambda]s] - no[\[Lambda]s0]);

(* 信号光和闲频光的波长关系 *)
\[Lambda]i[\[Lambda]s_] := 1/(1/\[Lambda]p - 1/\[Lambda]s);

(* 信号光的折射率 *)
ni[\[Lambda]s_] := no[\[Lambda]s];

(* 有效增益系数 *)
gamma[\[Alpha]_] := 
  4 \[Pi] deff[\[Alpha]] Sqrt[Ip/(2 \[Epsilon]0 np ns ni[\[Lambda]s])];
deff[\[Alpha]_] := d22 Cos[\[Theta] + \[Rho]];
d22 = 1.87*10^-10; (* BBO晶体的有效非线性系数, 单位: m/V *)
\[Epsilon]0 = 8.854187817*10^-12; (* 真空介电常数, 单位: F/m *)
np = ne[\[Lambda]p]; (* 抽运光的折射率 *)
ns = no[\[Lambda]s0]; (* 信号光的折射率 *)
\[Theta] = ArcCos[(no[\[Lambda]s0]^2 - ne[\[Lambda]p]^2 Sin[\[Alpha]]^2)/(no[\[Lambda]s0]^2 - ne[\[Lambda]p]^2)];
\[Rho] = ArcSin[(ne[\[Lambda]p] - no[\[Lambda]s0]) Sin[\[Alpha]]/ne[\[Lambda]p]];

(* 计算增益和相位 *)
A[\[Lambda]s_, \[Alpha]_] := (\[Delta]k[\[Lambda]s, \[Alpha]] L)/2;
B[\[Lambda]s_, \[Alpha]_] := Sqrt[(gamma[\[Alpha]] L)^2 - A[\[Lambda]s, \[Alpha]]^2];
G[\[Lambda]s_, \[Alpha]_] := 1 + (gamma[\[Alpha]] L)^2 (Sinh[B[\[Lambda]s, \[Alpha]]]/B[\[Lambda]s, \[Alpha]])^2;
\[Phi][\[Lambda]s_, \[Alpha]_] := ArcTan[
   B[\[Lambda]s, \[Alpha]] Sin[A[\[Lambda]s, \[Alpha]]] Cosh[B[\[Lambda]s, \[Alpha]]] - 
    A[\[Lambda]s, \[Alpha]] Cos[A[\[Lambda]s, \[Alpha]]] Sinh[B[\[Lambda]s, \[Alpha]]],
   B[\[Lambda]s, \[Alpha]] Cos[A[\[Lambda]s, \[Alpha]]] Cosh[B[\[Lambda]s, \[Alpha]]] + 
    A[\[Lambda]s, \[Alpha]] Sin[A[\[Lambda]s, \[Alpha]]] Sinh[B[\[Lambda]s, \[Alpha]]]
 ];

### 步骤2: 绘制增益谱和相位谱

接下来绘制不同非共线夹角下的增益谱和相位谱：

(* 绘制增益谱 *)
Plot[Evaluate@Table[G[\[Lambda]s, \[Alpha]], {\[Alpha], {0.041, 0.0415, 0.042, 0.0425}}], 
 {\[Lambda]s, 750*10^-9, 850*10^-9}, 
 PlotRange -> All, 
 Frame -> True, 
 FrameLabel -> {"Signal Wavelength (m)", "Gain"}, 
 PlotLegends -> Placed[{"\[Alpha] = 0.041", "\[Alpha] = 0.0415", "\[Alpha] = 0.042", "\[Alpha] = 0.0425"}, Below]]

(* 绘制相位谱 *)
Plot[Evaluate@Table[\[Phi][\[Lambda]s, \[Alpha]], {\[Alpha], {0.041, 0.0415, 0.042, 0.0425}}], 
 {\[Lambda]s, 750*10^-9, 850*10^-9}, 
 PlotRange -> All, 
 Frame -> True, 
 FrameLabel -> {"Signal Wavelength (m)", "Phase Shift"}, 
 PlotLegends -> Placed[{"\[Alpha] = 0.041", "\[Alpha] = 0.0415", "\[Alpha] = 0.042", "\[Alpha] = 0.0425"}, Below]]

### 步骤3: 绘制不同抽运光强度下的增益谱和相位谱

接下来绘制不同抽运光强度下的增益谱和相位谱：

(* 绘制增益谱 *)
Plot[Evaluate@Table[G[\[Lambda]s, 0.0415] /. Ip -> ip, {ip, {250*10^6, 350*10^6, 450*10^6}}], 
 {\[Lambda]s, 750*10^-9, 850*10^-9}, 
 PlotRange -> All, 
 Frame -> True, 
 FrameLabel -> {"Signal Wavelength (m)", "Gain"}, 
 PlotLegends -> Placed[{"I_p = 250 MW/cm^2", "I_p = 350 MW/cm^2", "I_p = 450 MW/cm^2"}, Below]]

(* 绘制相位谱 *)
Plot[Evaluate@Table[\[Phi][\[Lambda]s, 0.0415] /. Ip -> ip, {ip, {250*10^6, 350*10^6, 450*10^6}}], 
 {\[Lambda]s, 750*10^-9, 850*10^-9}, 
 PlotRange -> All, 
 Frame -> True, 
 FrameLabel -> {"Signal Wavelength (m)", "Phase Shift"}, 
 PlotLegends -> Placed[{"I_p = 250 MW/cm^2", "I_p = 350 MW/cm^2", "I_p = 450 MW/cm^2"}, Below]]

这些代码将生成类似于图1所示的增益谱和相位谱。你可以根据需要调整参数和绘图范围，以便更好地复现原始图像。