clear; clc;
%% 1. 固定参数定义
eta = 0.5;
P_c = 10e-6; % 控制场泵浦功率
hbar = 1.0545718e-34;% 约化普朗克常数
omega_c = 2 * pi * 3e8 / 1573e-9; % 控制场角频率
omega_m = 2 * pi * 3.68e9; % 机械振子角频率
epsilon_c = sqrt(P_c / (hbar * omega_c)); % 控制场振幅
p_1 = 0.05 * P_c; % 探测场1泵浦功率
p_2 = 0.05 * P_c; % 探测场2泵浦功率
Delta = omega_m; % 失谐量（示例值）
kappa = 0.1 * omega_m;% 腔耗散率
g = 2*pi*1e6; % 光-力耦合强度
gamma_1 = 0.005*omega_m; % 原子耗散率1
gamma_2 = 0.005*omega_m; % 原子耗散率2
J = 0.5* (gamma_1 + gamma_2); % 固定J为 0.5(γ₁+γ₂)
i = 1i; % 复数单位
tolerance = 1e-20; % 收敛容差
max_iterations = 100000; % 最大迭代次数
%% 2. 坐标轴范围与网格生成
delta1_ratios = linspace(-2, 2, 500); % X轴：δ₁/ωₘ 范围（匹配目标图刻度）
delta2_ratios = linspace(-2, 2, 500); % Y轴：δ₂/ωₘ 范围（匹配目标图刻度）
[X, Y] = meshgrid(delta1_ratios, delta2_ratios); % 生成网格
Z = zeros(size(X)); % 存储log10(ηₛ⁻)的矩阵
%% 3. 遍历网格计算ηₛ⁻
for i_row = 1:size(X, 1)
for j_col = 1:size(X, 2)
delta1_ratio = X(i_row, j_col); % 当前δ₁/ωₘ
delta2_ratio = Y(i_row, j_col); % 当前δ₂/ωₘ
delta1 = delta1_ratio * omega_m; % δ₁ = (δ₁/ωₘ) * ωₘ
delta2 = delta2_ratio * omega_m; % δ₂ = (δ₂/ωₘ) * ωₘ（原delta2改为变量）
% --- 迭代求解a和b1（δ₁、δ₂变化时，重新迭代保证收敛）---
a = 1; b1 = 1;
for iter = 1:max_iterations
a_prev = a; b1_prev = b1;
D = (i * omega_m + gamma_1/2) * (i * omega_m - gamma_2/2) + J^2;
b1 = (-i * g * abs(a)^2 * (-i * omega_m + gamma_2/2)) / D;
denominator = i * Delta + kappa/2 - i * g * (b1 + conj(b1));
a = (sqrt(eta*kappa) * epsilon_c) / denominator;
if abs(a - a_prev) < tolerance && abs(b1 - b1_prev) < tolerance
break; % 收敛则退出迭代
end
end
% --- 计算δ₂相关的辅助量（δ₂为变量，每次重新计算）---
D_plus_delta2 = -i * delta2 + i * omega_m + gamma_1 / 2;
D_minus_delta2 = -i * delta2 - i * omega_m + gamma_1 / 2;
E_plus_delta2 = -i * delta2 + i * omega_m - gamma_2 / 2;
E_minus_delta2 = -i * delta2 - i * omega_m - gamma_2 / 2;
F_plus_delta2 = -i * delta2 + (i*Delta - i*g*(conj(b1)+b1)) + kappa / 2;
F_minus_delta2 = -i * delta2 - (i*Delta - i*g*(conj(b1)+b1)) + kappa / 2;
L_plus_delta2 = D_plus_delta2 * E_plus_delta2 + J^2;
L_minus_delta2 = D_minus_delta2 * E_minus_delta2 + J^2;
alpha_delta2 = -E_minus_delta2 * L_plus_delta2 + E_plus_delta2 * L_minus_delta2;
beta_delta2 = L_plus_delta2 * L_minus_delta2 * F_minus_delta2 - g^2 * a * conj(a) * alpha_delta2;
% --- 计算y_omega_ratio（对应δ₁/ωₘ，因 y = -δ₁ = -delta1_ratio*omega_m）---
y_omega_ratio = delta1_ratio;
% --- 调用函数计算ηₛ⁻ ---
[~, ~, ~, ~, ~, eta_s_minus, ~, ~] = calculate_variables(...
y_omega_ratio, a, b1, beta_delta2, gamma_1, gamma_2, J, g, kappa, omega_m, Delta, eta, hbar, p_1, p_2, delta2, ...
D_plus_delta2, D_minus_delta2, E_plus_delta2, E_minus_delta2, F_plus_delta2, F_minus_delta2, L_plus_delta2, L_minus_delta2, alpha_delta2, omega_c);
Z(i_row, j_col) = log10(eta_s_minus); % 存储对数值
end
end
%% 4. 三维绘图与美化
figure('Position', [100, 100, 800, 600]);
surf(X, Y, Z); % 绘制曲面
shading interp; % 平滑着色（使颜色过渡自然）
colormap(jet); % 颜色映射（蓝→红，匹配目标图）
colorbar; % 显示颜色条
clim([-14, -2]); % 颜色范围（匹配目标图的-14到-2）
xlabel('$\delta_1/\omega_m$', 'Interpreter', 'latex', 'FontSize', 12);
ylabel('$\delta_2/\omega_m$', 'Interpreter', 'latex', 'FontSize', 12);
zlabel('$\log_{10}\eta_s^-$', 'Interpreter', 'latex', 'FontSize', 12);
title('(b)', 'FontSize', 16);
view(30, 30); % 调整视角（匹配目标图视角）
grid on; % 显示网格
%% 5. 辅助函数：计算ηₛ⁻（需放在脚本末尾）
function [sigma, psi, phi1, phi2, phi3, eta_s_minus, epsilon_1, y_omega_ratio] = calculate_variables(...
y_omega_ratio, a, b1, beta_delta2, gamma_1, gamma_2, J, g, kappa, omega_m, Delta, eta, hbar, p_1, p_2, delta2, ...
D_plus_delta2, D_minus_delta2, E_plus_delta2, E_minus_delta2, F_plus_delta2, F_minus_delta2, ~, ~, alpha_delta2, omega_c)
y = -y_omega_ratio * omega_m; % y = -δ₁（因 y_omega_ratio = δ₁/ωₘ）
% --- 计算D_+(y), D_-(y), E_+(y), E_-(y), F_+(y), F_-(y) ---
D_plus_y = -1i * y + 1i * omega_m + gamma_1/2;
D_minus_y = -1i * y - 1i * omega_m + gamma_1/2;
E_plus_y = -1i * y + 1i * omega_m - gamma_2/2;
E_minus_y = -1i * y - 1i * omega_m - gamma_2/2;
F_plus_y = -1i * y + (1i*Delta - 1i*g*(conj(b1)+b1)) + kappa/2;
F_minus_y = -1i * y - (1i*Delta - 1i*g*(conj(b1)+b1)) + kappa/2;
% --- 计算L_+(y), L_-(y), alpha(y), beta(y) ---
L_plus_y = D_plus_y * E_plus_y + J^2;
L_minus_y = D_minus_y * E_minus_y + J^2;
alpha_y = -E_minus_y * L_plus_y + E_plus_y * L_minus_y;
beta_y = L_plus_y * L_minus_y * F_minus_y - g^2 * a * conj(a) * alpha_y;
% --- 计算探测场1振幅epsilon_1 ---
epsilon_1 = sqrt(p_1 / (hbar * (y + omega_c)));
% --- 计算A1_plus, A1_minus_conj ---
numerator_A1_plus = sqrt(eta * kappa) * epsilon_1 * beta_y;
denominator_A1_plus = F_plus_y * beta_y + (g^2) * a * conj(a) * alpha_y * F_minus_y;
A1_plus = numerator_A1_plus / denominator_A1_plus * (abs(denominator_A1_plus)>=1e-20); % 避免除零
numerator_A1_minus_conj = g^2 * conj(a)^2 * alpha_y * A1_plus;
denominator_A1_minus_conj = beta_y;
A1_minus_conj = numerator_A1_minus_conj / denominator_A1_minus_conj * (abs(denominator_A1_minus_conj)>=1e-20);
% --- 计算B1_plus, B1_minus_conj ---
numerator_B1_plus = 1i * g * E_plus_y * (a * A1_minus_conj + conj(a) * A1_plus);
denominator_B1_plus = D_plus_y * E_plus_y + J^2;
B1_plus = numerator_B1_plus / denominator_B1_plus * (abs(denominator_B1_plus)>=1e-30);
numerator_B1_minus_conj = -1i * g * E_minus_y * (a * A1_minus_conj + conj(a) * A1_plus);
denominator_B1_minus_conj = D_minus_y * E_minus_y + J^2;
B1_minus_conj = numerator_B1_minus_conj / denominator_B1_minus_conj * (abs(denominator_B1_minus_conj)>=1e-20);
% --- 计算A2_plus, A2_minus_conj ---
numerator_A2_plus = sqrt(eta * kappa) * sqrt(p_2 / (hbar * (delta2 + omega_c))) * beta_delta2;
denominator_A2_plus = F_plus_delta2 * beta_delta2 + (g^2) * a * conj(a) * alpha_delta2 * F_minus_delta2;
A2_plus = numerator_A2_plus / denominator_A2_plus * (abs(denominator_A2_plus)>=1e-20);
numerator_A2_minus_conj = g^2 * conj(a)^2 * alpha_delta2 * A2_plus;
denominator_A2_minus_conj = beta_delta2;
A2_minus_conj = numerator_A2_minus_conj / denominator_A2_minus_conj * (abs(denominator_A2_minus_conj)>=1e-20);
% --- 计算B2_plus, B2_minus_conj ---
numerator_B2_plus = 1i * g * E_plus_delta2 * (a * A2_minus_conj + conj(a) * A2_plus);
denominator_B2_plus = D_plus_delta2 * E_plus_delta2 + J^2;
B2_plus = numerator_B2_plus / denominator_B2_plus * (abs(denominator_B2_plus)>=1e-20);
numerator_B2_minus_conj = -1i * g * E_minus_delta2 * (a * A2_minus_conj + conj(a) * A2_plus);
denominator_B2_minus_conj = D_minus_delta2 * E_minus_delta2 + J^2;
B2_minus_conj = numerator_B2_minus_conj / denominator_B2_minus_conj * (abs(denominator_B2_minus_conj)>=1e-20);
% --- 和频部分：计算sigma, phi1, phi2, phi3, psi, As_minus_conj ---
omega = y - delta2; % 和频角频率
D_plus_omega = -1i * omega + 1i * omega_m + gamma_1/2;
D_minus_omega = -1i * omega - 1i * omega_m + gamma_1/2;
E_plus_omega = -1i * omega + 1i * omega_m - gamma_2/2;
E_minus_omega = -1i * omega - 1i * omega_m - gamma_2/2;
F_plus_omega = -1i * omega + (1i*Delta - 1i*g*(conj(b1)+b1)) + kappa/2;
F_minus_omega = -1i * omega - (1i*Delta - 1i*g*(conj(b1)+b1)) + kappa/2;
L_plus_omega = D_plus_omega * E_plus_omega + J^2;
L_minus_omega = D_minus_omega * E_minus_omega + J^2;
alpha_omega = -E_minus_omega * L_plus_omega + E_plus_omega * L_minus_omega;
beta_omega = L_plus_omega * L_minus_omega * F_minus_omega - g^2 * abs(a)^2 * alpha_omega;
sigma = (D_minus_omega * E_minus_omega + J^2) * (D_plus_omega * E_plus_omega + J^2);
phi1 = A1_minus_conj*(B2_plus + B2_minus_conj) + A2_minus_conj*(B1_plus + B1_minus_conj);
phi2 = A1_plus*(B2_plus + B2_minus_conj) + A2_plus*(B1_plus + B1_minus_conj);
phi3 = A1_minus_conj*A2_plus + A2_minus_conj*A1_plus;
psi = 1i * g * (F_minus_omega * a * sigma * phi1 + beta_omega * conj(a) * phi2 - beta_omega * a * phi1);
% --- 计算As_minus_conj ---
numerator_As_minus_conj = g^2 * conj(a)^2 * alpha_omega * A1_plus + g^2 * conj(a) * alpha_omega*phi3 - 1i * g * sigma * phi1;
denominator_As_minus_conj = beta_omega;
As_minus_conj = numerator_As_minus_conj / denominator_As_minus_conj * (abs(denominator_As_minus_conj)>=1e-40);
% --- 计算ηₛ⁻ ---
eta_s_minus = abs(-sqrt(eta * kappa) * conj(As_minus_conj) / epsilon_1);
% 返回变量（未使用的用~占位，需保持数量一致）
sigma = sigma;
psi = psi;
phi1 = phi1;
phi2 = phi2;
phi3 = phi3;
y_omega_ratio = y_omega_ratio;
epsilon_1 = epsilon_1;
end优化我的代码
最新发布
扫一扫
8万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
