Unable to plot array inside a running while loop

2 次查看(过去 30 天)
This is the code
% snyder paper 1981
% parallel waveguides Morishita K.,Optics Letters,Vol 13. No.2
% 15 th April 2020
clearvars
nc1 = 1.449734; % core
nclad1 = 1.44479; % cladding - surrounding
delta1 = 1 - (nclad1/nc1)^2;
a1 = 4.75;% core radius
nc2 = 1.449686; % core
nclad2 = 1.4428; % cladding - surrounding
delta2 = 1 - (nclad2/nc2)^2;
a2 = 4.75; %core radius
d =11.5;% core to core center-distance
lambda = 1.25; % start wavelength in microns
inc = 0.005; % increment for wavelength
%% code to store data
kk = (1.605-lambda)/inc + 1;
kk = round(kk,0);
pp1 = zeros(1,kk);
pp2 = zeros(1,kk);
wave = zeros(1,kk);
coup = zeros(1,kk);
length = zeros(1,kk);
betas = zeros(1,kk);
i=0;
for z = 500:250:1250
z % check value of z
while lambda<=1.605
k0 = 2*3.14/lambda;
V1 = sqrt((a1^2)*(k0^2)*(nc1^2-nclad1^2));
U1 = -0.0024*(V1^4)+0.0493*(V1^3)-0.3757*(V1^2)+1.3444*V1-0.0095;
beta1 = sqrt((1/(a1^2))*((V1^2)/delta1)*(1-delta1*((U1/V1)^2)));
W1 = sqrt((a1^2)*(beta1^2-(k0^2)*(nclad1^2)));
V2 = sqrt((a2^2)*(k0^2)*(nc2^2-nclad2^2));
U2 = -0.0024*(V2^4)+0.0493*(V2^3)-0.3757*(V2^2)+1.3444*V2-0.0095;
beta2 = sqrt((1/(a2^2))*((V2^2)/delta2)*(1-delta2*((U2/V2)^2)));
W2 = sqrt((a2^2)*(beta2^2-(k0^2)*(nclad2^2)));
deltabeta = beta1-beta2;
A1 = ((delta1*delta2)^(1/4))/(sqrt(a1*a2));
A2 = (U1*U2)/((V1*V2)^(3/2));
A3 = W1*d/a1;
A4 = besselk(0,A3,1);
A5 = (besselk(1,W1,1))*(besselk(1,W2,1));
C = A1*A2*A4/A5;
L = 3.14/(2*C);
betad = deltabeta/2;
betab = sqrt(betad^2+C^2);
F = 1/(1+(betad/C)^2);
P1= 1-F*((sin(betab*z))^2);
P2 = F*((sin(betab*z))^2);
i = i +1;
pp1(:,i) = P1*100;
pp2(:,i) = P2*100;
wave(:,i) = lambda*1000;
coup(:,i) = C;
length(:,i) = L;
betas(:,i) = deltabeta*10^3;
figure(1)
plot(wave(:,i),pp1(:,i),'r--') % plot the current value of pp1 against wave
hold on
lambda = lambda + inc;
end
i = 0 % reset i for next iteration
lambda = 1.25 % reset lambda for next iteration
end
hold off

采纳的回答

Sriram Tadavarty
Sriram Tadavarty 2020-4-21
编辑:Sriram Tadavarty 2020-4-21
Hi Agnivo,
A minor update to your code in placing the plot after the while loops solves this. I commented out the lines that are not required.
clearvars
nc1 = 1.449734; % core
nclad1 = 1.44479; % cladding - surrounding
delta1 = 1 - (nclad1/nc1)^2;
a1 = 4.75;% core radius
nc2 = 1.449686; % core
nclad2 = 1.4428; % cladding - surrounding
delta2 = 1 - (nclad2/nc2)^2;
a2 = 4.75; %core radius
d =11.5;% core to core center-distance
lambda = 1.25; % start wavelength in microns
inc = 0.005; % increment for wavelength
%% code to store data
kk = (1.605-lambda)/inc + 1;
kk = round(kk,0);
pp1 = zeros(1,kk);
pp2 = zeros(1,kk);
wave = zeros(1,kk);
coup = zeros(1,kk);
length = zeros(1,kk);
betas = zeros(1,kk);
i=0;
for z = 500:250:1250
z % check value of z
while lambda<=1.605
k0 = 2*3.14/lambda;
V1 = sqrt((a1^2)*(k0^2)*(nc1^2-nclad1^2));
U1 = -0.0024*(V1^4)+0.0493*(V1^3)-0.3757*(V1^2)+1.3444*V1-0.0095;
beta1 = sqrt((1/(a1^2))*((V1^2)/delta1)*(1-delta1*((U1/V1)^2)));
W1 = sqrt((a1^2)*(beta1^2-(k0^2)*(nclad1^2)));
V2 = sqrt((a2^2)*(k0^2)*(nc2^2-nclad2^2));
U2 = -0.0024*(V2^4)+0.0493*(V2^3)-0.3757*(V2^2)+1.3444*V2-0.0095;
beta2 = sqrt((1/(a2^2))*((V2^2)/delta2)*(1-delta2*((U2/V2)^2)));
W2 = sqrt((a2^2)*(beta2^2-(k0^2)*(nclad2^2)));
deltabeta = beta1-beta2;
A1 = ((delta1*delta2)^(1/4))/(sqrt(a1*a2));
A2 = (U1*U2)/((V1*V2)^(3/2));
A3 = W1*d/a1;
A4 = besselk(0,A3,1);
A5 = (besselk(1,W1,1))*(besselk(1,W2,1));
C = A1*A2*A4/A5;
L = 3.14/(2*C);
betad = deltabeta/2;
betab = sqrt(betad^2+C^2);
F = 1/(1+(betad/C)^2);
P1= 1-F*((sin(betab*z))^2);
P2 = F*((sin(betab*z))^2);
i = i +1;
pp1(:,i) = P1*100;
pp2(:,i) = P2*100;
wave(:,i) = lambda*1000;
coup(:,i) = C;
length(:,i) = L;
betas(:,i) = deltabeta*10^3;
% figure(1)
% plot(wave(:,i),pp1(:,i),'r--') % plot the current value of pp1 against wave
% hold on
lambda = lambda + inc;
end
figure
plot(wave,pp1,'r--') % plot the complete wave for the iteration
i = 0 % reset i for next iteration
lambda = 1.25 % reset lambda for next iteration
end
Hope this helps.
Regards,
Sriram

更多回答(1 个)

David Hill
David Hill 2020-4-21
If you want everything stored:
nc1 = 1.449734; % core
nclad1 = 1.44479; % cladding - surrounding
delta1 = 1 - (nclad1/nc1)^2;
a1 = 4.75;% core radius
nc2 = 1.449686; % core
nclad2 = 1.4428; % cladding - surrounding
delta2 = 1 - (nclad2/nc2)^2;
a2 = 4.75; %core radius
d =11.5;% core to core center-distance
lambda = 1.25; % start wavelength in microns
inc = 0.005; % increment for wavelength
%% code to store data
kk = (1.605-lambda)/inc + 1;
kk = round(kk,0);
pp1 = zeros(1,kk);
pp2 = zeros(1,kk);
wave = zeros(1,kk);
coup = zeros(1,kk);
length = zeros(1,kk);
betas = zeros(1,kk);
i=0;
count=0;
for z = 500:250:1250
%z % check value of z
count=count+1;
while lambda<=1.605
k0 = 2*3.14/lambda;
V1 = sqrt((a1^2)*(k0^2)*(nc1^2-nclad1^2));
U1 = -0.0024*(V1^4)+0.0493*(V1^3)-0.3757*(V1^2)+1.3444*V1-0.0095;
beta1 = sqrt((1/(a1^2))*((V1^2)/delta1)*(1-delta1*((U1/V1)^2)));
W1 = sqrt((a1^2)*(beta1^2-(k0^2)*(nclad1^2)));
V2 = sqrt((a2^2)*(k0^2)*(nc2^2-nclad2^2));
U2 = -0.0024*(V2^4)+0.0493*(V2^3)-0.3757*(V2^2)+1.3444*V2-0.0095;
beta2 = sqrt((1/(a2^2))*((V2^2)/delta2)*(1-delta2*((U2/V2)^2)));
W2 = sqrt((a2^2)*(beta2^2-(k0^2)*(nclad2^2)));
deltabeta = beta1-beta2;
A1 = ((delta1*delta2)^(1/4))/(sqrt(a1*a2));
A2 = (U1*U2)/((V1*V2)^(3/2));
A3 = W1*d/a1;
A4 = besselk(0,A3,1);
A5 = (besselk(1,W1,1))*(besselk(1,W2,1));
C = A1*A2*A4/A5;
L = 3.14/(2*C);
betad = deltabeta/2;
betab = sqrt(betad^2+C^2);
F = 1/(1+(betad/C)^2);
P1= 1-F*((sin(betab*z))^2);
P2 = F*((sin(betab*z))^2);
i = i +1;
pp1(count,i) = P1*100;
pp2(count,i) = P2*100;
wave(count,i) = lambda*1000;
coup(count,i) = C;
length(count,i) = L;
betas(count,i) = deltabeta*10^3;
lambda = lambda + inc;
end
hold on;
plot(wave(count,:),pp1(count,:),'r--')
i = 0; % reset i for next iteration
lambda = 1.25; % reset lambda for next iteration
end

类别

Help CenterFile Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by