Unable to plot array inside a running while loop
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Agnivo Gosai
el 21 de Abr. de 2020
Editada: Agnivo Gosai
el 21 de Abr. de 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 comentarios
Respuesta aceptada
Sriram Tadavarty
el 21 de Abr. de 2020
Editada: Sriram Tadavarty
el 21 de Abr. de 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
Más respuestas (1)
David Hill
el 21 de Abr. de 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 comentario
Ver también
Categorías
Más información sobre Solver Outputs and Iterative Display en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!