Unable to plot array inside a running while loop

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

Answer 1

[Sriram Tadavarty]

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

[Agnivo Gosai]

Thanks

Answer 2

[David Hill]

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

[Agnivo Gosai]

Thanks, this is very helpful