Unable to plot array inside a running while loop
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Agnivo Gosai
le 21 Avr 2020
Modifié(e) : Agnivo Gosai
le 21 Avr 2020
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
0 commentaires
Réponse acceptée
Sriram Tadavarty
le 21 Avr 2020
Modifié(e) : Sriram Tadavarty
le 21 Avr 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
Plus de réponses (1)
David Hill
le 21 Avr 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 commentaire
Voir également
Catégories
En savoir plus sur Solver Outputs and Iterative Display dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!