I want to add only values of Q_sun in each iteration answer, should not consider NaN

clear
clc
n = 172;                       %no. of day
H = 20000;                     % altitude in meter 
CF = 0.5;                      %cloud fraction
V_wind =10;                    %wind velocity   
%%Constants taken
u = [7.447323    2.071808    9.010095    7.981491 194];
a = u(1); % 7.447323;
b = u(2); % 2.071808;
c = u(3); % 9.010095;
d = u(4); % 7.981491;
l = u(5); % 194;

%%constants
alpha = 0.15;                                              %absorptivity of film
albedo_g = 0.35;                                           %ground albedo                                     
albedo_c = 0.5;                                            %cloud albedo
P_sea = 101325;                                            %pressure at sea level
e = 0.0167;                                                %eccentricity of earth
T_g = 280;                                                 %temperature of ground in k
T_cl = 262.15;                                             % assuming cloud at altitude 4km in k
R = 6371000;                                               %redius of earth
epsilon_e = 0.95;                                          %earth emissivity
alpha_ir = 0.85;                                           %infrared absorbtivity of film 
sigma = 5.67.*10.^-8;                                      %stefan boltzman constant  
g = 9.81;
% [T_oa, P_oa, rho_oa]  = atmos(H);
syms x t
Fun = 1/8*(sqrt(u(1)*(u(5)-x).*(u(2).*x-u(5)*sqrt(u(3))+sqrt(u(3)*u(5).^2-u(4)*u(5).*x))));
h=1;
k = 1;
zi = k:k:360;
x=0:h:l;                                                             
Fun1 = matlabFunction(Fun);
area = 0;
Q_sun = 0;
Q_earth = 0;
Q_sky = 0;
Q_IR = 0;
Q_aIR = 0;
Q_oc = 0;
for j=1:length(zi)
    for i=1:(length(x)-1)
        y1 = Fun1(x(i));
        z1 = Fun1(x(i));
        x2 = x(i+1);
        y2 = Fun1(x(i+1));
        z2 = Fun1(x(i+1));
        x3 = x(i+1);
        y3 = y2*sin((pi/180)*zi(j));
        z3 = z2*cos((pi/180)*zi(j));
        c1 = y2*cos((pi/180)*k);
        c2 = z2*sin((pi/180)*k);
        T1 = [0, y3-y2, z3-z2];
        T2 = [x(i)-x(i+1), y1-y2, z1-z2];
        vec = cross(T1,T2);                                                % normal vector
        deltal = sqrt((y2-y1)^2+h^2);
        deltazi = sqrt(c2^2+(y2-c1)^2);
        dA11 = deltal.*deltazi;                                            
        area = area +dA11;
        %%angles 
        w = 15.*(12-10);                                                    % hour angle   
        phi = 18.975;                                                      % lattitude 
        delta = 23.45.*sin((pi/180).*(360/365).*(284+n));%declination angle
        omega = asin(sin((pi/180).*delta).*sin((pi/180).*phi)+cos((pi/180).*delta).*cos((pi/180).*phi).*cos((pi/180).*w)); %elevation angle
        psi = asin((sin(pi/180).*w).*cos((pi/180).*delta)/cos(omega));     % azimuth angle 
        thetai = (pi/2)+asin(vec(2)/(sqrt(vec(1)^2+vec(3)^2)));         
        E = [cos(omega)*cos(psi), -sin(psi)*cos(omega), sin(omega)];       % solar vector 
        beta = acos((dot(E,vec)/(norm(E).*norm(vec))));                    % angle between two vector
        MA = 2.*pi.*n/365;
        M = (5466/P_sea).*(sqrt(1229+(614.*sin(omega)).^2)-614.*sin(omega));
        Tau_atm = 0.5.*(exp(-0.65.*M)+exp(-0.95.*M));
        zita = MA+2.*e.*sin(M)+1.25.*e.^2.*sin(2.*M); 
        ID = (1367.*Tau_atm.^M).*((1+e.*cos(zita))./(1-e.^2)).^2;
        Q_D = alpha.*dA11.*ID.*sin(omega).*cos(pi-beta);                   % direct radiation in Watt

          I_d = (1/2).*ID.*sin(omega).*(1.-Tau_atm.^M)./(1.-1.4.*log(Tau_atm));
          Q_d = alpha.*I_d.*dA11.*(1+cos(thetai))/2;                         % diffuse radiation

          %Reflected solar radiation model
          albedo = albedo_g.*(1-CF).^2+albedo_c.*CF;
          Ir = albedo.*(ID.*sin(omega)+I_d);
          Q_r = alpha.*Ir.*dA11.*(1-cos(thetai))/2;                          % reflected radiation
          if beta <= pi/2                                                    % this loop is for consideration of direct radiation,  
              Q_sun=Q_sun+Q_D+Q_d+Q_r;                                       
          else Q_sun = Q_sun+Q_d+Q_r;
          end
      end 
  end
  disp(Q_sun)

code gives Ans : NaN

Answer 1

[Walter Roberson]

          if beta <= pi/2                                                    % this loop is for consideration of direct radiation,  
             to_add = [Q_d, Q_d, Q_r];
          else
             to_add = [Q_d, Q_r];
          end
          Q_sun = Q_sun + to_add(~isnan(to_add));