複数のfor及びif​文を含むプログラムを​実行し​たが出てこな​い答え​がある。

下記のコードを実行しました。

本当は、（Zs、Zp1、Zr1、Zp2、Zr2、mn1、mn2）＝（17、25、67、28、76、1.2、1.05）も答えになるはずなのですが、答えとして出力されません。（これだけではありませんが…）

傾向を見ていると、mn1とmn2が1.2と1.05の組み合わせはそもそも結果で出てこないです。

%モジュール%
for mn1=1:0.05:1.2
   for mn2=1:0.05:1.2

の範囲を変える（0.8-1.2など）とこの現象が起きます。

原因は何でしょうか。浮動小数点誤差などが影響しているのでしょうか？また、回避の方法はありますでしょうか。

▼コード

clear

 K=1;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%固定パラメータ%%%%%%%%%%%%%%%%%%%%%%%%%% 
    B1=0;%%%1段目ねじれ角(°）
    B2=0;%%%2段目ねじれ角(°）
    An1=20;%%%1段目歯直角圧力角An(°)
    An2=20;%%%2段目歯直角圧力角An(°)
    mu=0.1;%%%歯面摩擦係数
    N=100;%%歯数MAX
    n=3;%%プラネタリ個数
    gl=380;%%減速比下限
    gh=410;%%減速比上限
    epsilon_l=1.2;%%かみ合い率下限値
    D=98;%%外径制約
    
    Xs=0;
    Xp1=0;
    Xr1=0;
    Xp2=0;
    Xr2=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%求めたい解%%%%%%%%%%%%%%%%%
%モジュール%
for mn1=1:0.05:1.2
   for mn2=1:0.05:1.2
       
%歯数%
for Zs=5:N
    for Zp1=5:N
        for Zr1=5:N
            for Zp2=5:N
                for Zr2=5:N

 
                    
 %%%%%%%%%%%%%%%%%%%%%%%%%%条件式%%%%%%%%%%%%%%%%%%%%%%%%%%

 I1=Zr1/Zs;
 I2=(Zr1*Zp2)/(Zr2*Zp1);
 
 
 %%%条件1%%%
     if Zr1<=Zp1
       continue
     end
    
     if Zr2<=Zp2
       continue
     end
 
 %%%条件2%%%
  g=1/((1-I2)/(1+I1)); 
    if g<=gl
       continue
    end
                    
    if g>=gh
       continue
    end 
    
    
 %%%条件3%%%
  A=(Zs+Zr1)/n;
    if  A~=floor(A)
        continue
    end
    
 %%%条件4%%%
 mt1=mn1/cos(deg2rad(B1));
 mt2=mn2/cos(deg2rad(B2));
 
    if mt1*(Zs+Zp1)~=mt1*(Zr1-Zp1)
        continue
    end
    
    if mt1*(Zr1-Zp1)~=mt2*(Zr2-Zp2)
        continue
    end
    
 %%%条件5%%%
  P1= rem(Zr1,n)/n;
  P2= rem(Zr2,n)/n;
  np= n*gcd(Zp1,Zp2);
 
  P11=zeros(np-1,1);
  P22=zeros(np-1,1);
    for i0=0:(np-1)
        P11(i0+1)=(rem(i0*Zp1,np))/np;
        P22(i0+1)=(rem(i0*Zp2,np))/np;
    end
                        
   P_1= P1==P11;
   P_2= P2==P22;
                        
   A_P1=find(P_1);
   A_P2=find(P_2);
                       
    if isempty(A_P1)==1 || isempty(A_P2)==1
       continue
    end
    
    A_P=zeros(length(A_P1),length(A_P2));
    for i=1:length(A_P1)
        for ii=1:length(A_P2)
            A_P(i,ii)=A_P1(i)==A_P2(ii);
        end
    end
    
    AA_P=find(A_P);
    
    if isempty(AA_P)==1
        continue
    end
    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%サン-プラ①%%%%%%%%%%%%%%%%%%%%%%%%%% 

 At1=atan(tan(deg2rad(An1))/cos(deg2rad(B1)));%%%正面圧力角At(rad)
 At11=rad2deg(At1);%%%正面圧力角At(°)
 I_At1=tan(At1)-(pi()*At11)/180;%%%インボリュートAt(rad)
 I_Awt1=((2*(Xs+Xp1))/(Zs+Zp1))*tan(deg2rad(An1))+I_At1;%%%インボリュートAwt(rad)
 
 err1=1;
 i1=0;
 while err1>0.000001
     i1=i1+1;
     
     if i1 == 1
          X(i1)=1+(I_Awt1-tan(1)+1)/tan(1)^2;
          continue
     end
     
 X(i1)=X(i1-1)+(I_Awt1-tan(X(i1-1))+X(i1-1))/tan(X(i1-1))^2;
 
  err1=abs(X(i1)-X(i1-1));
  
  if i1 > 15
      return
  end
  
 end
 
 Awt1=rad2deg(X(end));%%%%正面かみ合い圧力角(°)
 
 ds=(Zs*mn1)/cos(deg2rad(B1));%%%サン基準円直径(㎜)
 dp1=(Zp1*mn1)/cos(deg2rad(B1));%%%プラ①基準円直径(㎜)
 dbs=ds*cos(At1);%%%サン基礎円直径(㎜)
 dbp1=dp1*cos(At1);%%%プラ①基礎円直径(㎜)
 dws=dbs/cos(deg2rad(Awt1));%%%サンかみ合いピッチ円直径(㎜)
 dwp1=dbp1/cos(deg2rad(Awt1));%%%プラ①かみ合いピッチ円直径(㎜)
 rs=dws/2;%%%サンかみ合いピッチ円半径(㎜)
 rp1=dwp1/2;%%%プラ①かみ合いピッチ円半径(㎜)
 ysp1=((Zs+Zp1)/(2*cos(deg2rad(B1))))*((cos(At1)/cos(deg2rad(Awt1)))-1);%%%中心距離修正係数
 has=(1+ysp1-Xp1)*mn1;%%%サン歯末のたけ
 hap1=(1+ysp1-Xs)*mn1;%%%プラ①歯末のたけ
 das=ds+2*has;%%%サン歯先円直径
 dap1=dp1+2*hap1;%%%プラ①歯先円直径
 dfs=ds-2*has;%%%サン歯底円直径
 dfp1=dp1-2*hap1;%%%プラ①歯底円直径
 asp1=((Zs+Zp1)/(2*cos(deg2rad(B1)))+ysp1)*mn1;%%%中心距離
 alpha_s=acos(dbs/das);%%%サン歯先円圧力角(rad)
 alpha_p1=acos(dbp1/dap1);%%%プラ①歯先円圧力角(rad)
 epsilon_a1=(Zp1/(2*pi()))*(tan(alpha_p1)-tan(deg2rad(Awt1)));%%%かみ合い因子率εα1
 epsilon_a2=(Zs/(2*pi()))*(tan(alpha_s)-tan(deg2rad(Awt1)));%%%かみ合い因子率εα2
 epsilon_a=epsilon_a1^2+epsilon_a2^2-epsilon_a1-epsilon_a2+1;%%%かみ合い因子率εα
 
 eta_alpha=1-mu*pi()*((1/Zs)+(1/Zp1))*epsilon_a;%%%サン-プラ①基準効率
 
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%プラ①-リング①%%%%%%%%%%%%%%%%%%%%%%%%%% 

 At2=atan(tan(deg2rad(An1))/cos(deg2rad(B1)));%%%正面圧力角At(rad)
 At22=rad2deg(At2);%%%正面圧力角At(°)
 I_At2=tan(At2)-(pi()*At22)/180;%%%インボリュートAt(rad)
 I_Awt2=((2*(Xr1-Xp1))/(Zr1-Zp1))*tan(deg2rad(An1))+I_At2;%%%インボリュートAwt(rad)
 
 err2=1;
 i2=0;
 while err2>0.000001
     i2=i2+1;
     if i2 == 1
          Y(i2)=1+(I_Awt2-tan(1)+1)/tan(1)^2;
          continue
     end
 Y(i2)=Y(i2-1)+(I_Awt2-tan(Y(i2-1))+Y(i2-1))/tan(Y(i2-1))^2;
 
  err2=abs(Y(i2)-Y(i2-1));
  
  if i2 > 15
      return
  end
  
 end
 
 Awt2=rad2deg(Y(end));%%%%正面かみ合い圧力角(°)
 
 dp11=(Zp1*mn1)/cos(deg2rad(B1));%%%プラ①基準円直径(㎜)
 dr1=(Zr1*mn1)/cos(deg2rad(B1));%%%リング①基準円直径(㎜)
 dbp11=dp11*cos(At2);%%%プラ①基礎円直径(㎜)
 dbr1=dr1*cos(At2);%%%リング①基礎円直径(㎜)
 dwp11=dbp11/cos(deg2rad(Awt2));%%%プラ①かみ合いピッチ円直径(㎜)
 dwr1=dbr1/cos(deg2rad(Awt2));%%%リング①かみ合いピッチ円直径(㎜)
 rp11=dwp11/2;%%%プラ①かみ合いピッチ円半径(㎜)
 rr1=dwr1/2;%%%リング①かみ合いピッチ円半径(㎜)
 yp1r1=((Zr1-Zp1)/(2*cos(deg2rad(B1))))*((cos(At2)/cos(deg2rad(Awt2)))-1);%%%中心距離修正係数
 hap11=(1+Xp1)*mn1;%%%プラ①歯末のたけ
 har1=(1-Xr1)*mn1;%%%リング①歯末のたけ
 dap11=dp11+2*hap11;%%%プラ①歯先円直径
 dar1=dr1-2*har1;%%%リング①歯先円直径
 dfp11=dp11-2*hap11;%%%プラ①歯底円直径
 dfr1=dr1+2*har1;%%%リング①歯底円直径
 ap1r1=((Zr1-Zp1)/(2*cos(deg2rad(B1)))+yp1r1)*mn1;%%%中心距離
 alpha_p11=acos(dbp11/dap11);%%%プラ①歯先円圧力角(rad)
 alpha_r1=acos(dbr1/dar1);%%%リング①歯先円圧力角(rad)
 epsilon_b1=-(Zr1/(2*pi()))*(tan(alpha_r1)-tan(deg2rad(Awt2)));%%%かみ合い因子率εβ1
 epsilon_b2=(Zp1/(2*pi()))*(tan(alpha_p11)-tan(deg2rad(Awt2)));%%%かみ合い因子率εβ2
 epsilon_b=epsilon_b1^2+epsilon_b2^2-epsilon_b1-epsilon_b2+1;%%%かみ合い因子率εβ
 
 eta_beta=1-mu*pi()*((1/Zp1)-(1/Zr1))*epsilon_b;%%%プラ①-リング①基準効率
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%プラ②-リング②%%%%%%%%%%%%%%%%%%%%%%%%%% 

 At3=atan(tan(deg2rad(An2))/cos(deg2rad(B2)));%%%正面圧力角At(rad)
 At33=rad2deg(At3);%%%正面圧力角At(°)
 I_At3=tan(At3)-(pi()*At33)/180;%%%インボリュートAt(rad)
 I_Awt3=((2*(Xr2-Xp2))/(Zr2-Zp2))*tan(deg2rad(An2))+I_At3;%%%インボリュートAwt(rad)
 
 err3=1;
 i3=0;
 while err3>0.000001
     i3=i3+1;
     if i3 == 1
          Z(i3)=1+(I_Awt3-tan(1)+1)/tan(1)^2;
          continue
     end
 Z(i3)=Z(i3-1)+(I_Awt3-tan(Z(i3-1))+Z(i3-1))/tan(Z(i3-1))^2;
 
  err3=abs(Z(i3)-Z(i3-1));
  
  if i3 > 15
      return
  end
  
 end
 
 Awt3=rad2deg(Z(end));%%%%正面かみ合い圧力角(°)
 
 dp2=(Zp2*mn2)/cos(deg2rad(B2));%%%プラ②基準円直径(㎜)
 dr2=(Zr2*mn2)/cos(deg2rad(B2));%%%リング②基準円直径(㎜)
 dbp2=dp2*cos(At3);%%%プラ②基礎円直径(㎜)
 dbr2=dr2*cos(At3);%%%リング②基礎円直径(㎜)
 dwp2=dbp2/cos(deg2rad(Awt3));%%%プラ②かみ合いピッチ円直径(㎜)
 dwr2=dbr2/cos(deg2rad(Awt3));%%%リング②かみ合いピッチ円直径(㎜)
 rp2=dwp2/2;%%%プラ②かみ合いピッチ円半径(㎜)
 rr2=dwr2/2;%%%リング②かみ合いピッチ円半径(㎜)
 yp2r2=((Zr2-Zp2)/(2*cos(deg2rad(B2))))*((cos(At3)/cos(deg2rad(Awt3)))-1);%%%中心距離修正係数
 hap2=(1+Xp2)*mn2;%%%プラ②歯末のたけ
 har2=(1-Xr2)*mn2;%%%リング②歯末のたけ
 dap2=dp2+2*hap2;%%%プラ②歯先円直径
 dar2=dr2-2*har2;%%%リング②歯先円直径
 dfp2=dp2-2*hap2;%%%プラ②歯底円直径
 dfr2=dr2+2*har2;%%%リング②歯底円直径
 ap2r2=((Zr2-Zp2)/(2*cos(deg2rad(B2)))+yp2r2)*mn2;%%%中心距離
 alpha_p2=acos(dbp2/dap2);%%%プラ②歯先円圧力角(°)
 alpha_r2=acos(dbr2/dar2);%%%リング②歯先円圧力角(°)
 epsilon_g1=-(Zr2/(2*pi()))*(tan(alpha_r2)-tan(deg2rad(Awt3)));%%%かみ合い因子率εγ1
 epsilon_g2=(Zp2/(2*pi()))*(tan(alpha_p2)-tan(deg2rad(Awt3)));%%%かみ合い因子率εγ2
 epsilon_g=epsilon_g1^2+epsilon_g2^2-epsilon_g1-epsilon_g2+1;%%%かみ合い因子率εγ
 
 eta_gamma=1-mu*pi()*((1/Zp2)-(1/Zr2))*epsilon_g;%%%プラ②-リング②基準効率
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%効率算出%%%%%%%%%%%%%%%%%%%%%%%%%%
    if I2<1
      eta=((1-I2)/(1+I1))*((1+(eta_alpha*eta_beta*I1))/(1-(eta_beta*eta_gamma*I2))); %%%I2<1 順駆動
      eta_gyaku=((1+I1)/(1-I2))*((eta_alpha*(eta_beta*eta_gamma-I2))/(eta_gamma*(eta_alpha*eta_beta+I1)));%%%I2<1 逆駆動
   
    elseif I1>1
      eta=((1-I2)/(1+I1))*((eta_gamma*(eta_beta+eta_alpha*I1))/(eta_beta*eta_gamma-I2));%%%I1>1 順駆動
      eta_gyaku=((1+I1)/(1-I2))*((eta_alpha*(1-eta_beta*eta_gamma*I2))/(eta_alpha+eta_beta*I1));%%%I1>1 逆駆動
    end
    
    if isreal(eta) == 0 || isreal(eta_gyaku) == 0
        continue
    end
 

 %%%%%%%%%%%%%%%%%%%%%%%%%%条件式%%%%%%%%%%%%%%%%%%%%%%%%%%
 

%%条件6%%%
    if dfr1+2*mn1>=D
       continue
    end
    
    if dfr2+2*mn1>=D
    continue
    end  


%%%条件7%%%
    if dap1/(asp1*sin(pi()/n))>=2
        continue
    end

    if dap2/(asp1*sin(pi()/n))>=2
        continue
    end


%%%条件8%%%
    if epsilon_a1+epsilon_a2<epsilon_l
        continue
    end    
    
    if epsilon_b1+epsilon_b2<epsilon_l
        continue
    end    
   
    if epsilon_g1+epsilon_g2<epsilon_l
        continue
    end    

%%%条件9%%%
   inv1=Zp1/((1-(tan(alpha_r1)/tan(deg2rad(Awt2))))*Zr1);
   inv2=Zp2/((1-(tan(alpha_r2)/tan(deg2rad(Awt3))))*Zr2);
   
   if inv1<1
       continue
   end
   
   if inv2<1
       continue
   end

%%%条件10%%%
  sita_p1=acos(((dar1/2)^2-(dap11/2)^2-ap1r1^2)/2/ap1r1/(dap11/2))+tan(alpha_p11)-alpha_p11-I_Awt2;
  sita_r1=acos((ap1r1^2+(dar1/2)^2-(dap11/2)^2)/2/ap1r1/(dar1/2));
  sita_p2=acos(((dar2/2)^2-(dap2/2)^2-ap2r2^2)/2/ap2r2/(dap2/2))+tan(alpha_p2)-alpha_p2-I_Awt3;
  sita_r2=acos((ap2r2^2+(dar2/2)^2-(dap2/2)^2)/2/ap2r2/(dar2/2));
  
  tro1=sita_p1*(Zp1/Zr1)+I_Awt2-(tan(alpha_r1)-(pi()*rad2deg(alpha_r1))/180);
  tro2=sita_p2*(Zp2/Zr2)+I_Awt3-(tan(alpha_r2)-(pi()*rad2deg(alpha_r2))/180);
  
  if tro1<sita_r1
      continue
  end
  
  if tro2<sita_r2
      continue
  end
  


                 Answer(K,:)=[Zs Zp1 Zr1 Zp2 Zr2 mn1 mn2 g eta eta_gyaku];
                    K=K+1;   

                    
                end
            end
        end
    end
end
   end
end       

[実香]

一応、条件4は以下のように修正しました。

しかし、mn1とmn2の探索範囲を広げたりすると、本来出てくるべき答えが出てきません…

 %%%条件4%%%
 mt1=mn1/cos(deg2rad(B1));%%1段目正面モジュール
 mt2=mn2/cos(deg2rad(B2));%%2段目正面モジュール
 
    if abs(mt1*(Zs+Zp1)-mt1*(Zr1-Zp1))>1e-5
        continue
    end
    
    if abs(mt1*(Zr1-Zp1)-mt2*(Zr2-Zp2))>1e-5
        continue
    end