Unable to plot array inside a running while loop
Show older comments
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
Accepted Answer
Sriram Tadavarty
on 21 Apr 2020
Edited: Sriram Tadavarty
on 21 Apr 2020
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 Comment
Agnivo Gosai
on 21 Apr 2020
More Answers (1)
David Hill
on 21 Apr 2020
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
1 Comment
Agnivo Gosai
on 21 Apr 2020
Edited: Agnivo Gosai
on 21 Apr 2020
Categories
Find more on Deep Learning Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
