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
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.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!