电光调Q设计仿真python代码

import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate import solve_ivp

import warnings

import time

# 屏蔽不必要的警告

warnings.filterwarnings('ignore')

# ============================================================

# 1. 自动中文字体配置 (解决乱码)

# ============================================================

def set_chinese_font():

    import matplotlib as mpl

    # 优先尝试黑体或微软雅黑

    fonts = ['SimHei', 'Microsoft YaHei', 'Arial Unicode MS', 'sans-serif']

    for font in fonts:

        try:

            plt.rcParams['font.sans-serif'] = [font]

            plt.rcParams['axes.unicode_minus'] = False

            # 简单测试绘图

            fig = plt.figure()

            plt.close(fig)

            break

        except:

            continue

set_chinese_font()

# ============================================================

# 2. 物理常数与激光参数 (Nd:YVO4 @ 1064nm)

# ============================================================

c        = 3.0e8

h        = 6.626e-34

lambda_  = 1064e-9

nu       = c / lambda_

lambda_p = 808e-9

# Nd:YVO4 晶体参数

sigma_e  = 25e-23       # 受激辐射截面 (m²)

tau_f    = 90e-6        # 上能级寿命 (s)

l_c      = 10e-3        # 晶体长度 (m)

w_c      = 200e-6       # 模场半径 (m)

A_c      = np.pi * w_c**2

# 腔参数

L_cav    = 100e-3       # 腔长 (m)

R_oc     = 0.85         # 输出镜反射率

T_oc     = 1.0 - R_oc

L_int    = 0.03         # 腔内损耗

t_r      = 2.0 * L_cav / c # 往返时间

# 预计算系数

delta    = L_int + 0.5 * T_oc

t_c_on   = t_r / (2.0 * delta)

eta_pump = (lambda_p / lambda_) * 0.9 * 0.8

E_ph_p   = h * c / lambda_p

R_coeff  = eta_pump / (E_ph_p * A_c * l_c)

Pout_coeff = T_oc / (2.0 * delta) * (A_c * L_cav) * (h * nu) / t_c_on

beta     = 1e-8

# ============================================================

# 3. 核心仿真引擎 (速率方程)

# ============================================================

def rhs_qswitch(t, y, R_pump):

    N, phi = y

    sc = sigma_e * c

    dN   = R_pump - N / tau_f - sc * N * phi

    dphi = (sc * N * (l_c/L_cav) - 1.0/t_c_on) * phi + beta * N / tau_f

    return [dN, dphi]

def simulate_single(P_pump, f_rep=10e3, duty=0.9):

    """单脉冲全过程仿真"""

    T_rep = 1.0 / f_rep

    t_pump = duty * T_rep

    # 阶段1：解析解计算积累后的初始反转粒子数

    N_ss = (R_coeff * P_pump) * tau_f

    N_init = N_ss * (1.0 - np.exp(-t_pump / tau_f))

    # 阶段2：自适应步长求解脉冲

    t_max = 3e-6 # 仿真 3us 窗口

    sol = solve_ivp(rhs_qswitch, (0, t_max), [N_init, 1e8], 

                    args=(R_coeff * P_pump,), method='RK45', rtol=1e-7)

    t = sol.t

    N = sol.y[0]

    phi = sol.y[1]

    P_out = np.maximum(phi * Pout_coeff, 0)