利用numpy库实现激光器领域的时域行波模型 该模型输出激光器信号，频谱等.zip

#!/usr/bin/python # -*- coding: UTF-8 -*- import time import random import math import numpy import matplotlib.pyplot as plt # 数学常量 PI = 3.1415926 # 物理常量： SPEED_LIGHT = 2.997e10 # cm/s ELECTRON_CHARGE = 1.602e-19 # C H = 6.63e-34 # 材料参数： KACHA = 0.08 A = 1.846*1e9 # 1/s B = 1.26*1e-10 # m^3/s C = 3.553*1e-29 # m^6/s N_EFF = 3.2263 N_G = 3.4034 G_N = 2000 # cm N0 = 0.5e18 # cm-3 BATA = 0.5e-4 AFA = 12 # cm-1 EPSILON = 15e-17 # cm3 AFA_M = 2.5 # 结构参数 R_RIGHT = math.sqrt(0.005) R_LEFT = math.sqrt(0.95) L_CAV = 250e-4 # cm W_CAV = 1.8e-4 # cm D_ACT = 125*1e-9*100 # cm KAPPA0 = 70 # cm-1 LAMBDA_BRAG = 1577 # nm # 仿真参数 SUB_NUM = 50 # 输入参数 MOD_RATE = 10*(2**30) # Hz I_BIAS = 36e-3 # A I_MOD = 33.8e-3 BIT_NUM = 200 TIME_START = 1.5e-9 # s # 初始化参数 V_EFF = SPEED_LIGHT/N_EFF V_G = SPEED_LIGHT/N_G DELTA_Z = L_CAV/SUB_NUM DELTA_T = DELTA_Z/V_G I_INJECT = I_BIAS KAPPA_RIGHT = KAPPA0 * numpy.ones(int(SUB_NUM * 0.6), dtype=complex, order='C') KAPPA_LEFT = KAPPA0 * numpy.ones(SUB_NUM - int(SUB_NUM*0.6), dtype=complex, order='C') KAPPA = numpy.hstack((KAPPA_RIGHT, KAPPA_LEFT)) L_SUB = numpy.linspace(0.0, L_CAV * 1e4, SUB_NUM + 1) N_CARRIER = N0*numpy.linspace(1, 1, SUB_NUM) PHOTON_DENSITY = numpy.zeros(SUB_NUM, dtype=float, order='C') F_LIGHT = numpy.zeros(SUB_NUM + 1, dtype=complex, order='C') R_LIGHT = numpy.zeros(SUB_NUM + 1, dtype=complex, order='C') BIT_TIME = 1.0/MOD_RATE BIT_DOTS_NUM = int(BIT_TIME//DELTA_T) BIT_IDX_START = int(TIME_START//DELTA_T) MOD_DOTS_NUM = int((BIT_NUM*BIT_DOTS_NUM)/2)*2 + 1 ALL_DOTS = BIT_IDX_START + MOD_DOTS_NUM SIG_ARRAY = numpy.zeros(ALL_DOTS, dtype=complex, order='C') DENSITY_TO_POWER = numpy.sqrt(H*SPEED_LIGHT*SPEED_LIGHT/N_G*W_CAV*D_ACT*(1 - R_RIGHT**2)\ /LAMBDA_BRAG/KACHA*1e10) def __outer_circulation(): global I_INJECT, N_CARRIER, PHOTON_DENSITY bit_rem = 0 for idx_outer in range(0, ALL_DOTS): if idx_outer > BIT_IDX_START: if 0 == bit_rem: k = random.randint(0, 1) I_INJECT = I_BIAS + I_MOD/2*(-1)**k # I_INJECT = I0 + mod_I / 2; # I_INJECT = I0 - mod_I / 2; bit_rem = (idx_outer-BIT_IDX_START) % BIT_DOTS_NUM J = I_INJECT/L_CAV/W_CAV gain0 = G_N * numpy.log(N_CARRIER / N0) gain1 = G_N * numpy.log(N_CARRIER / N0)/(1 + EPSILON * PHOTON_DENSITY) delta_n = -((LAMBDA_BRAG*1e-7/4/PI)*AFA_M*KACHA)*gain0 delta = 2*PI*1e7/LAMBDA_BRAG*delta_n G = 0.5*KACHA*gain1 - AFA/2 - 1j*delta a11 = 1/numpy.cosh(KAPPA * DELTA_Z) a12 = 1j*numpy.tanh(KAPPA * DELTA_Z) a21 = a12 a22 = a11 F_LIGHT[0] = R_RIGHT*R_LIGHT[0] F0 = F_LIGHT F_LIGHT[1:]=numpy.exp(DELTA_Z*G)*F0[: -1] R0 = R_LIGHT Fz = F_LIGHT[1:] R_LIGHT[-1] = R_LEFT*F0[-1] R_LIGHT[:-1]=numpy.exp(DELTA_Z*G)* R0[1:] tao = 15 * numpy.sqrt(BATA*B/L_CAV/V_EFF*(N_CARRIER**2)) # if idx_outer > BIT_IDX_START: # tao = 0 SF = tao*numpy.random.randn(1, SUB_NUM)* numpy.exp(1j*2*PI*numpy.random.rand(1,SUB_NUM)) F_LIGHT[1:]=a11*F_LIGHT[1:] + a12*R_LIGHT[:-1]+DELTA_Z*SF tao = 15 * numpy.sqrt(BATA * B / L_CAV / V_EFF * (N_CARRIER ** 2)) # if idx_outer > BIT_IDX_START: # tao = 0 SR = tao*numpy.random.randn(1, SUB_NUM)* numpy.exp(1j*2*PI*numpy.random.rand(1,SUB_NUM)) R_LIGHT[:-1]=a21*Fz + a22*R_LIGHT[:-1]+DELTA_Z * SR PHOTON_DENSITY = 0.5*((numpy.abs(F_LIGHT[:-1]))**2 + (numpy.abs(R_LIGHT[:-1]))**2 + (numpy.abs(F_LIGHT[1:]))**2 + (numpy.abs(R_LIGHT[1:]))**2) N_CARRIER = N_CARRIER + (J/ELECTRON_CHARGE/D_ACT - A*N_CARRIER - B*(N_CARRIER**2) - C*(N_CARRIER**3) - V_G*gain1*PHOTON_DENSITY)*DELTA_T SIG_ARRAY[idx_outer] = DENSITY_TO_POWER*R_LIGHT[0] if 0 == idx_outer%(50*BIT_DOTS_NUM): print(idx_outer//BIT_DOTS_NUM) def __fft(time_single): NFFT = 2**int((math.ceil(numpy.log2(numpy.shape(time_single)[0])))) fre_single = numpy.fft.fft(time_single, NFFT) return NFFT, 10*numpy.log10(numpy.abs(fre_single/numpy.max(numpy.abs(fre_single)))**2) def __get_time(second): second = int(second) seconds = second %60 all_mins = math.floor(second/60) mins = all_mins%60 hours = math.floor(all_mins/60) return "%2d:%2d:%2d" %(hours, mins, seconds) if __name__ == '__main__': time_start = time.time() __outer_circulation() plot_sig = numpy.abs(SIG_ARRAY)**2 Time = 1e12*DELTA_T*numpy.linspace(0, ALL_DOTS-1, ALL_DOTS) # ps plt.close() plt.figure(figsize=(32, 18)) plt.plot(Time,plot_sig, linestyle="-", color="blue") plt.xlim(Time[0], Time[-1]) max_sig = numpy.max(plot_sig) min_sig = numpy.min(plot_sig) plt.ylim(min_sig, max_sig) plt.xticks(numpy.linspace(Time[0], Time[-1], 15, endpoint=True), fontsize=30) plt.yticks(numpy.linspace(min_sig, max_sig, 10, endpoint=True),fontsize=30) plt.xlabel("Time(ps)", fontsize=40) plt.ylabel("Power(mW)", fontsize=40) plt.savefig('./功率图.png') # plt.show() if BIT_NUM>20: fft_dots = 20*BIT_DOTS_NUM else: fft_dots = BIT_NUM*BIT_DOTS_NUM NFFT, fre_single = __fft(SIG_ARRAY[-fft_dots:]) delta_v_fft = 1/DELTA_T*numpy.linspace(-0.5, 0.5, NFFT) delta_lambda_fft = -((10**-9*LAMBDA_BRAG)**2)*delta_v_fft/(SPEED_LIGHT*0.01)*10**9 + LAMBDA_BRAG plt.figure(figsize=(32, 18)) plt.plot(delta_lambda_fft, fre_single, linestyle="-", color="red") plt.xticks(fontsize=30) plt.yticks(fontsize=30) plt.xlabel("lambda(nm)", fontsize=40) plt.ylabel("Power(dB)", fontsize=40) # plt.savefig('./频谱.png') plt.show() sig_start = BIT_IDX_START+4*BIT_DOTS_NUM numpy.save("single_out.npy", plot_sig[sig_start:]) numpy.save("DELTA_T.npy", DELTA_T) numpy.save("BIT_DOTS_NUM.npy", BIT_DOTS_NUM) time_consume = __get_time(time.time() - time_start) print("耗时：%s" % time_consume)