How to skip NaN arrays outputs of a for loop?
    5 views (last 30 days)
  
       Show older comments
    
clc; clear all; close all;
syms m1 m2 s1 s2 j1 j2 
lambda = 1060*10^-9;
M=1
z=linspace(0.00001,8000);
wo = 0.02;
C = 10^(-7);
k=2*pi/lambda;
po=(0.545*C^2*k^2*z).^(-(3/5));
b=0.1;
T1=0;
T2=0;
T3=0;
T4=0;
T5=0;
T6=0;
delta= ((1i*k)./(2*z))+(1./(wo.^2))+(1./(po.^2));
delta_star= subs(delta, 1i, -1i);
etta = delta_star - (1./(delta.*po.^4));
X=((k./(2.*z)).^2).*((1./(2.*1i.*sqrt(delta))).^M).*(1./(16.*delta.*etta));
alpha = (1i.*k./(2.*z)).*(1./(delta.*po.^2)-1);
beta_p = (b./(2.*wo)).*(1./(delta.*po.^2)+1);
beta_n = (b./(2.*wo)).*(1./(delta.*po.^2)-1);
A = (alpha.^2./etta)-((k.^2)-(4.*z.^2.*delta));
B1_p = ((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_p)./etta);
B1_n = ((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_n)./etta);
B2_p = -(((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_p)./etta));
B2_n = -(((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_n)./etta));
C1_p = ((b.^2)./(4.*wo.^2.*delta))+(((beta_p).^2)./etta);
C1_n = ((b.^2)./(4.*wo.^2.*delta))+(((beta_n).^2)./etta);
C2_p = ((b.^2)./(4.*wo.^2.*delta))+(((beta_n).^2)./etta);
C2_n = ((b.^2)./(4.*wo.^2.*delta))+(((beta_p).^2)./etta);
m1 = 0:M;
m2=0:M;
s1 = m1./2;
s2=(M-m1)./2;
%%%% j1=m1-2.*s1 %%%%
allVecs_1 = {-2*s1,m1}; 
sub = cell(1,numel(allVecs_1));
[sub{:}] = ndgrid(allVecs_1{:});
sub = cellfun(@(x)x(:),sub,'UniformOutput', false);
j= cell2mat(sub);
j11 = [ sum(j(:,[1 2]),2) ];
j1=j11(j11>=0);
%%%% j2=(M-m1-(2.*s2)) %%%%
allVecs_1 = {-1*m1,-2*s2};
sub = cell(1,numel(allVecs_1));
[sub{:}] = ndgrid(allVecs_1{:});
sub = cellfun(@(x)x(:),sub,'UniformOutput', false);
j21= cell2mat(sub);
j2212 =M + [ sum(j21(:,[1 2]),2) ];
j2=j2212(j2212>=0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%
[rm1,cm1] = size(m1);
[rm2,cm2] = size(m2);
[rs1,cs1] = size(s1);
[rs2,cs2] = size(s2);
[rj1,cj1] = size(j1);
[rj2,cj2] = size(j2);
for Nm1 = 1:cm1
    for Mm1 = 1:rm1
        M1= m1(Mm1,Nm1);
        T2=0;
    for Nm2 = 1:cm2
        for Mm2 = 1:rm2
            M2= m2(Mm2,Nm2);
            T3=0;
        for Ns1 = 1:cs1
            for Ms1 = 1:rs1
               S1= s1(Ms1,Ns1);
               T4=0;
               for Nj1 = 1:cj1
                    for Mj1 = 1:rj1
                        J1= j1(Mj1,Nj1);
                        T5=0;
                        for Ns2 = 1:cs1
                            for Ms2 = 1:rs2
                                S2= s2(Ms2,Ns2);
                                T6=0;
                                for Nj2 = 1:cj2
                                    for Mj2 = 1:rj2
                                        J2= j2(Mj2,Nj2);
                                        e3=(exp(C1_p).*hermiteH(M-M2+J2,((1i.*beta_p)./sqrt(etta))))+(exp(C2_n).*hermiteH(M-M2+J2,((-1i.*beta_n)./sqrt(etta))));
                                        e4=(exp(C1_n).*hermiteH(M-M2+J2,((1i.*beta_n)./sqrt(etta))))+(exp(C2_p).*hermiteH(M-M2+J2,((-1i.*beta_p)./sqrt(etta))));
                                        T6=T6+(gamma(1+M-M1-(2.*S2))/(gamma(J2+1)*gamma(M-M1-(2.*S2)-J2+1))).*((1./(2.*1i.*sqrt(etta))).^(M-M2+J2)).*((1./(po.^2)).^J2).*((b./2.*wo).^(M-M1-2*S2-J2)).*(e4+e3);
                                    end
                                end  
                                e1=(exp(C1_n).*hermiteH(M2+J1,((1i.*beta_n)./sqrt(etta))))+(exp(C2_p).*hermiteH(M2+J1,((-1i.*beta_p)./sqrt(etta))));
                                e2=(exp(C1_p).*hermiteH(M2+J1,((1i.*beta_p)./sqrt(etta))))+(exp(C2_n).*hermiteH(M2+J1,((-1i.*beta_n)./sqrt(etta))));
                                T5=T5+T6.*(((gamma(M-M1+1).*(-1).^S2)./(gamma(S2+1).*gamma(M-M1-2.*S2+1)))).*(((2.*1i)./(delta.^(0.5))).^(M-M1-2.*S2)).*(e1+e2);
                                T5 = rmmissing(vpa(T5))
                            end
                        end
                        T4=T4+T5.*((gamma(M1-(2.*S2)+1))/(gamma(J1+1)*gamma(1+M1-(2.*S2)-J1))).*((1./(2.*1i.*sqrt(etta))).^(M2+J2)).*((1./(po.^2)).^J1).*((b./2.*wo).^(M1-2.*S1-J1));
                         T4 = rmmissing(vpa(T4))
                    end
               end
            end
            T3=T3+T4.*(((gamma(M1+1).*(-1).^S1)./(gamma(S1+1).*gamma(1+M1-2.*S1))).*(((2.*1i)./(delta.^(0.5))).^(M1-2.*S1)));
            T3 = rmmissing(vpa(T3))
        end
        T2=T2+T3*nchoosek(M,M2)*((-1i)^(M-M2));
        T2 = rmmissing(vpa(T2))
    end
    T1=T1+T2*nchoosek(M,M1)*((1i)^(M-M1));
    T1 = rmmissing(vpa(T1))
end
end
end
I0 = abs(real(X.*T1))
2 Comments
  Matt J
      
      
 on 18 Jul 2022
				What is your question? The title of your post doesn't convey what problems you're having with the posted code.
Accepted Answer
  Jan
      
      
 on 18 Jul 2022
        doBreak = false;
for Nj2 = 1:cj2
    for Mj2 = 1:rj2
        J2= j2(Mj2,Nj2);
        e3=(exp(C1_p).*hermiteH(M-M2+J2,((1i.*beta_p)./sqrt(etta))))+(exp(C2_n).*hermiteH(M-M2+J2,((-1i.*beta_n)./sqrt(etta))));
        e4=(exp(C1_n).*hermiteH(M-M2+J2,((1i.*beta_n)./sqrt(etta))))+(exp(C2_p).*hermiteH(M-M2+J2,((-1i.*beta_p)./sqrt(etta))));
        tmp = (gamma(1+M-M1-(2.*S2))/(gamma(J2+1)*gamma(M-M1-(2.*S2)-J2+1))).*((1./(2.*1i.*sqrt(etta))).^(M-M2+J2)).*((1./(po.^2)).^J2).*((b./2.*wo).^(M-M1-2*S2-J2)).*(e4+e3);
        if isnan(tmp)
            doBreak = true;
            break;
        end
        T6 = T6 + tmp;
    end
    if doBreak  % The loop over Nj2 also?
        break;
    end
end  
0 Comments
More Answers (0)
See Also
Categories
				Find more on Loops and Conditional Statements 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!

