Hi, everyone, why is the laser output power not shown in the figure?

33 次查看(过去 30 天)
function Pout = dualwavetest4(params,initialConditions, dt, dz, T, L,numSteps)
params = [2825 * 1e-9, 982 * 1e-9, 1620 * 1e-9, 2.136e-21*1e-4, 1.688e-21*1e-4,...
0.4e-21*1e-4,0.21e-21*1e-4,0.8e-3,6.9e-3,10e-6,120e-6,570e-6,0.37,...
0.63,0.99,0.00,0.001,0.85,0.006,0.004,0.14,0.34,0.012,0.015,0.18,...
0.44,1.0,2,2,0.1,0.16,2.8918e-10,1.28,0.362,0.542,1.55,1.12e-21*1e-4,...
2.31e-21*1e-4,3.02e-21*1e-4,3.2e+026,0.03,0.01625,0.0325,0.7511,0.0024,...
0.052,0,0.04,0,0.04,0.994,0.04,16,3e8,6.626e-34,10,5];
%params(1) = lambda_s = 2825 * 1e-9;
%params(2) = lambda_p1 = 982 * 1e-9;
%params(3) = lambda_p2 = 1620 * 1e-9;
%params(4) = sigma_ap1 = 2.136e-21*1e-4;
%params(5) = sigma_ep1 = 1.688e-21*1e-4;
%params(6) = sigma_ap2 = 0.4e-21*1e-4;
%params(7) = sigma_ep2 = 0.21e-21*1e-4;
%params(8) = tau1 = 0.8e-3;
%params(9) = tau2 = 6.9e-3;
%params(10)= tau3 = 10e-6;
%params(11)= tau4 = 120e-6;
%params(12)= tau5 = 570e-6;
%params(13)= beta21 = 0.37;
%params(14)= beta20 = 0.63;
%params(15)= beta32 = 0.99;
%params(16)= beta31 = 0.00;
%params(17)= beta30 = 0.001;
%params(18)= beta43 = 0.85;
%params(19)= beta42 = 0.006;
%params(20)= beta41 = 0.004;
%params(21)= beta40 = 0.14;
%params(22)= beta54 = 0.34;
%params(23)= beta53 = 0.012;
%params(24)= beta52 = 0.015;
%params(25)= beta51 = 0.18;
%params(26)= beta50 = 0.44;
%params(27)= beta10 = 1.0;
%params(28)= g1 = 2;
%params(29)= g2 = 2;
%params(30)= b1 = 0.1;
%params(31)= b2 = 0.16;
%params(32)= Aeff = 2.8918e-10;
%params(33)= W11 = 1.28;
%params(34)= W22 = 0.362;
%params(35)= W50 = 0.542;
%params(36)= W42 = 1.55;
%params(37)= sigma_25 = 1.12e-21*1e-4;
%params(38)= sigma_52 = 2.31e-21*1e-4;
%params(39)= sigma_es = 3.02e-21*1e-4;%
%
%
%params(40)= NEr = 3.2e+026;
%params(41)= alpha_p1 = 0.03;
%params(42)= alpha_p2 = 0.01625;
%params(43)= alpha_s = 0.0325;
%params(44)= gamma_s = 0.7511;
%params(45)= gamma_p1 = 0.0024; %
%params(46)= gamma_p2 = 0.052; %
%params(47)= Rp1 = 0;
%params(48)= Rp2 = 0.04;
%params(49)= Rp3 = 0;
%params(50)= Rp4 = 0.04;
%params(51)= Rs1 = 0.994;
%params(52)= Rs2 = 0.04;
%params(53)= L = 16;
%params(54)= c = 3e8;
%params(55)= h = 6.626e-34;
%Ppl = [Pinput_p1; Pinput_p2];
%params(56)= Pinput_p1 = 10;
%params(57)= Pinput_p2 = 5;
%params(58)= Ppr = 0;
N = zeros(6, 1); %
y = [params(56);0;params(57);0;0;0]; %
Z = linspace(0,params(53),100); %
t = linspace(0, 90, 100); %
numSteps = length(t);
dt = 0.01;
dz = 0.01;
y_history = zeros(length(y), numSteps);
N_history = zeros(6, numSteps);
for i = 1:numSteps
for j = 1:length(Z)-1%
% y = apply_boundary_conditions(y, params, Z(j));
% fprintf('Boundary conditions at Z=%.2f applied. Before: Pump powers: %.2e, %.2e; Laser power: %.2e\n', ...
% Z,y(1), y(3), y(5));
[N, y] = rk4_step(params, N, y, Z(j), dt, dz);
end
y_history(:, i) = y;
N_history(:, i) = N;
% y = apply_boundary_conditions(y, params, 0);
% y = apply_boundary_conditions(y, params, params(53));
% fprintf('After: Pump powers: %.2e, %.2e; Laser power: %.2e\n', ...
% y_history(1,:), y_history(3,:), y_history(5,:));
end
figure;
subplot(2,1,1);
plot(Z, y_history(1,:), 'b.-', 'DisplayName', 'Pp1+');
hold on;
plot(Z, y_history(2, :), 'g*-', 'DisplayName', 'Pp1-');
plot(Z, y_history(3, :), 'p-', 'DisplayName', 'Pp2+');
plot(Z, y_history(4, :), 'y>-', 'DisplayName', 'Pp2-');
plot(Z, y_history(5, :), 'r', 'DisplayName', 'Ps+');
plot(Z, y_history(6, :), 'k--', 'DisplayName', 'Ps-');
grid on;
title('Pump and laser powers');
legend('Location', 'best');
xlabel('Position z (m)');
ylabel('Power (W)');
subplot(2,1,2);
for k = 1:6
plot(Z, N_history(k, :), 'DisplayName', strcat('energy level ', num2str(k)));
hold on;
end
grid on;
title('Relative population density');
legend('Location', 'best');
xlabel('Position z (m)');
ylabel('N/N');
end
function y = apply_boundary_conditions(y, params, Z)
% if Z == 0 || Z == params(53)
% % Z = 0
% y(1) = params(47) * y(2) + params(56); % Pp1
% y(3) = params(49) * y(4) + params(57); % Pp2
%
% % Z = L
% y(2) = params(48) * y(1);
% y(4) = params(50) * y(3);
%
% % Boundary conditions
% y(5) = params(51) * y(6); % Emission conditions of laser power
% y(6) = params(52) * y(5);
% end
end
function [N, y] = rk4_step(params, N, y, Z, dt, dz)
params = [2825 * 1e-9, 982 * 1e-9, 1620 * 1e-9, 2.136e-21*1e-4, 1.688e-21*1e-4,...
0.4e-21*1e-4,0.21e-21*1e-4,0.8e-3,6.9e-3,10e-6,120e-6,570e-6,0.37,...
0.63,0.99,0.00,0.001,0.85,0.006,0.004,0.14,0.34,0.012,0.015,0.18,...
0.44,1.0,2,2,0.1,0.16,2.8918e-10,1.28,0.362,0.542,1.55,1.12e-21*1e-4,...
2.31e-21*1e-4,3.02e-21*1e-4,3.2e+026,0.03,0.01625,0.0325,0.7511,0.0024,...
0.052,0,0.04,0,0.04,0.994,0.04,16,3e8,6.626e-34,10,5];
k1_N = dt * rate_equations(N, y, params);
k1_y = dz * transmission_equations(N, y, params, Z);
k2_N = dt * rate_equations(N + 0.5 * k1_N, y + 0.5 * k1_y, params);
k2_y = dz * transmission_equations(N + 0.5 * k1_N, y + 0.5 * k1_y, params, Z);
k3_N = dt * rate_equations(N + 0.5 * k2_N, y + 0.5 * k2_y, params);
k3_y = dz * transmission_equations(N + 0.5 * k2_N, y + 0.5 * k2_y, params, Z);
k4_N = dt * rate_equations(N + k3_N, y + k3_y, params);
k4_y = dz * transmission_equations(N + k3_N, y + k3_y, params, Z);
N = N + (1/6) * (k1_N + 2*k2_N + 2*k3_N + k4_N);
y = y + (1/6) * (k1_y + 2*k2_y + 2*k3_y + k4_y);
end
function dNdt = rate_equations(N, y, params)
params = [2825 * 1e-9, 982 * 1e-9, 1620 * 1e-9, 2.136e-21*1e-4, 1.688e-21*1e-4,...
0.4e-21*1e-4,0.21e-21*1e-4,0.8e-3,6.9e-3,10e-6,120e-6,570e-6,0.37,...
0.63,0.99,0.00,0.001,0.85,0.006,0.004,0.14,0.34,0.012,0.015,0.18,...
0.44,1.0,2,2,0.1,0.16,2.8918e-10,1.28,0.362,0.542,1.55,1.12e-21*1e-4,...
2.31e-21*1e-4,3.02e-21*1e-4,3.2e+026,0.03,0.01625,0.0325,0.7511,0.0024,...
0.052,0,0.04,0,0.04,0.994,0.04,16,3e8,6.626e-34,10,5];
% RGSA = params(2)*params(45)*(params(4)*N(0)-params(5)*N(2))*(y(1)+y(2))/(params(55)*params(54)*params(32));
% RESA1 = params(2)*params(45)*(params(37)*N(2)-params(38)*N(5))*(y(1)+y(2))/(params(55)*params(54)*params(32));
% RESA2 = params(3)*params(46)*(params(6)*N(1)-params(7)*N(3))*(y(3)+y(4))/(params(55)*params(54)*params(32));
% RSE = params(1)*params(44)*params(39)*(b2*N(2)-g2/g1*b1*N(1))*(y(5)+y(6))/(params(55)*params(54)*params(32));
dNdt = zeros(6,1);
dNdt(6) = -N(6)/params(6)+params(6)*N(3)*N(3)-params(35)*N(6)*N(1)+params(36)*N(5)*N(3)+params(2)*params(45)*(params(37)*N(3)-params(38)*N(6))*(y(1)+y(2))/(params(55)*params(54)*params(32));
dNdt(5) = -N(5)/params(11)-params(36)*N(5)*N(3)+params(22)*N(6)/params(12);
dNdt(4) = -N(4)/params(10)+params(33)*N(2)*N(2)+params(35)*N(6)*N(1)+params(18)/params(11)*N(5)+params(23)/params(12)*N(6)+params(3)*params(46)*(params(6)*N(2)-params(7)*N(4))*(y(3)+y(4))/(params(55)*params(54)*params(32));
dNdt(3) = -N(3)/params(9)-2*params(34)*N(3)*N(3)-params(36)*N(5)*N(3)+params(15)/params(10)*N(4)+params(19)/params(11)*N(5)+params(24)/params(12)*N(6)+params(2)*params(45)*(params(4)*N(1)-params(5)*N(3))*(y(1)+y(2))/(params(55)*params(54)*params(32))-(params(2)*params(45)*(params(37)*N(3)-params(38)*N(6))*(y(1)+y(2))/(params(55)*params(54)*params(32)))-(params(1)*params(44)*params(39)*(params(31)*N(3)-params(29)/params(28)*params(30)*N(2))*(y(5)+y(6))/(params(55)*params(54)*params(32)));
dNdt(2) = -N(2)/params(8)-2*params(33)*N(2)*N(2)+params(35)*N(6)*N(1)+params(36)*N(5)*N(3)+params(13)/params(9)*N(3)+params(16)/params(10)*N(4)+params(20)/params(11)*N(5)+params(25)/params(12)*N(6)+params(1)*params(44)*params(39)*(params(31)*N(3)-params(29)/params(28)*params(30)*N(2))*(y(5)+y(6))/(params(55)*params(54)*params(32))-(params(3)*params(46)*(params(6)*N(2)-params(7)*N(4))*(y(3)+y(4))/(params(55)*params(54)*params(32)));
dNdt(1) = params(33)*N(2)*N(2)+params(34)*N(3)*N(3)-params(35)*N(6)*N(1)+params(27)/params(8)*N(2)+params(14)/params(9)*N(4)+params(17)/params(10)*N(3)+params(21)/params(11)*N(5)+params(26)/params(12)*N(6)-(params(2)*params(45)*(params(4)*N(1)-params(5)*N(3))*(y(1)+y(2))/(params(55)*params(54)*params(32)));
%N(6)+N(1)+N(2)+N(3)+N(4)+N(5)= params(40);
end
function dydz = transmission_equations(N, y, params, Z)
params = [2825 * 1e-9, 982 * 1e-9, 1620 * 1e-9, 2.136e-21*1e-4, 1.688e-21*1e-4,...
0.4e-21*1e-4,0.21e-21*1e-4,0.8e-3,6.9e-3,10e-6,120e-6,570e-6,0.37,...
0.63,0.99,0.00,0.001,0.85,0.006,0.004,0.14,0.34,0.012,0.015,0.18,...
0.44,1.0,2,2,0.1,0.16,2.8918e-10,1.28,0.362,0.542,1.55,1.12e-21*1e-4,...
2.31e-21*1e-4,3.02e-21*1e-4,3.2e+026,0.03,0.01625,0.0325,0.7511,0.0024,...
0.052,0,0.04,0,0.04,0.994,0.04,16,3e8,6.626e-34,10,5];
dydz = zeros(6,1);
dydz(1) = (-params(45)*(params(4)*N(1)-params(5)*N(3)+params(37)*N(3)-params(38)*N(6))-params(41))*y(1);
dydz(2) = (params(45)*(params(4)*N(1)-params(5)*N(3)+params(37)*N(3)-params(38)*N(6))+params(41))*y(2);
dydz(3) = (-params(46)*(params(6)*N(2)-params(7)*N(4))-params(42))*y(3);
dydz(4) = (params(46)*(params(6)*N(2)-params(7)*N(4))+params(42))*y(4);
dydz(5) = (params(44)*params(39)*(params(31)*N(3)-params(29)/params(28)*params(30)*N(2))-params(43))*y(5);
dydz(6) = (-params(44)*params(39)*(params(31)*N(3)-params(29)/params(28)*params(30)*N(2))+params(43))*y(6);
fprintf('At Z=%.2f, Rate of change for pump powers: %.2e, %.2e; Laser power: %.2e\n', ...
Z, dydz(1), dydz(3), dydz(5));
if Z == 0
% Z = 0
y(1) = params(47)*y(2)+params(56);%y(1) = Rp1*y(2)+P_input_p1;
y(3) = params(49)*y(4)+params(57);%y(3) = Rp3*y(4)+P_input_p2;
y(5) = params(51)*y(6);%y(5) = Rs1*y(6);
elseif Z == params(53)
y(2) = params(48)*y(1);%y(2) = Rp2*y(1);
y(4) = params(50)*y(3);%y(4) = Rp4*y(3);
y(6) = params(52)*y(5);%y(6) = Rs2*y(5);
end
fprintf('At Z=%.2f, Rate of change for pump powers: %.2e, %.2e; Laser power: %.2e\n', ...
Z, dydz(1), dydz(3), dydz(5));
end

采纳的回答

Walter Roberson
Walter Roberson 2024-7-15,23:35
plot(Z, y_history(4, :), 'y>-', 'DisplayName', 'Pp2-');
plot(Z, y_history(5, :), 'r', 'DisplayName', 'Ps+');
plot(Z, y_history(6, :), 'k--', 'DisplayName', 'Ps-');
All of those values are identical to zero.
dNdt(6) = -N(6)/params(6)+params(6)*N(3)*N(3)-params(35)*N(6)*N(1)+params(36)*N(5)*N(3)+params(2)*params(45)*(params(37)*N(3)-params(38)*N(6))*(y(1)+y(2))/(params(55)*params(54)*params(32));
Every term there depends upon a value from N. When N is initialized to all zero, the result must be identical to zero.
N = N + (1/6) * (k1_N + 2*k2_N + 2*k3_N + k4_N);
all-zero plus terms that are all-zero, gives all zero. So N stays all zero.
dydz(4) = (params(46)*(params(6)*N(2)-params(7)*N(4))+params(42))*y(4);
dydz(5) = (params(44)*params(39)*(params(31)*N(3)-params(29)/params(28)*params(30)*N(2))-params(43))*y(5);
dydz(6) = (-params(44)*params(39)*(params(31)*N(3)-params(29)/params(28)*params(30)*N(2))+params(43))*y(6);
y(4), y(5), y(6) happen to all be 0, and due to the multiplication by y(4), y(5), y(6), those derivatives all come out as 0. So the last 3 y values all come out as 0.
  1 个评论
LinHai
LinHai 2024-7-16,6:37
Thanks for your answer. Here, I set the initial value of the laser power to 2 W, logically speaking, the laser power will increase and stabilize to a fixed value, why the laser power keeps decreasing? And how should I set the initial value of N?

请先登录,再进行评论。

更多回答(0 个)

产品


版本

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by