I want values of Q_sun at different given time but got error: In an assignment A(I) = B, the number of elements in B and I must be the same.

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 
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=10;
k = 10;
zi = k:k:360;
x=0:h:l;                                                             
Fun1 = matlabFunction(Fun);
area = 0;
Q_D = 0;
Q_d=0;
Q_r=0;
Q_earth = 0;
Q_sky = 0;
Q_IR = 0;
Q_aIR = 0;
Q_oc = 0;
t = 10:1:14;
Q_sun = zeros(1,length(t));
for c = 1:(length(t));
    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-t(c));                                                    % 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;
              QD = alpha.*dA11.*ID.*sin(omega).*cos(pi-beta);                   % direct radiation in Watt
              Q_D = nansum(QD);

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

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

Answer 1

[Ingrid]

there are many problems with this code

for example:

your t vector has 5 elements so the length of Q_sun will also be 5, however you are trying to assign this vector to one element

Q_sun(c) = (Q_sun+Q_D+Q_d+Q_r);

there is also a problem with the two inner loops i and j as these end AFTER assigning to Q_sun(c)

after all it is not clear what you are trying to do but only that the current code is flawed at many points

[Ingrid]

in an attempt to clean up your code without knowing which processes you are modelling, it should look something like this:

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
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=10;
k = 10;
zi = k:k:360;
x=0:h:l;
Fun1 = matlabFunction(Fun);
area = 0;
Q_D = 0;
Q_d=0;
Q_r=0;
Q_earth = 0;
Q_sky = 0;
Q_IR = 0;
Q_aIR = 0;
Q_oc = 0;
t = 10:1:14;
Q_sun = zeros(1,length(t));
% not changing so should not be in loop
phi = 18.975;                                                      % lattitude
delta = 23.45.*sin((pi/180).*(360/365).*(284+n));                  %declination angle
MA = 2.*pi.*n/365;
albedo = albedo_g.*(1-CF).^2+albedo_c.*CF;
for c = 1:(length(t))
    % all these variables are only changing in the first loop so avoid
    % recalculating them in the next loops!
    w = 15.*(12-t(c)); % hour 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
    E = [cos(omega)*cos(psi), -sin(psi)*cos(omega), sin(omega)];       % solar vector -> is this supposed to give imaginary numbers????
    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;
    I_d = (1/2).*ID.*sin(omega).*(1.-Tau_atm.^M)./(1.-1.4.*log(Tau_atm));
    Ir = albedo.*(ID.*sin(omega)+I_d);
    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;

              thetai(i) = (pi/2)+asin(vec(2)/(sqrt(vec(1)^2+vec(3)^2)));
              beta(i) = acos((dot(E,vec)/(norm(E).*norm(vec))));                    % angle between two vector
          end
      end
          QD = alpha.*dA11.*ID.*sin(omega).*cos(pi-beta);                   % direct radiation in Watt
          Qd = alpha.*I_d.*dA11.*(1+cos(thetai))/2;                         % diffuse radiation
      Q_D = nansum(QD); % why take sum when only one value??? -> therefore I assume this should be after loops
      Q_d = nansum(Qd);

      Qr = alpha.*Ir.*dA11.*(1-cos(thetai))/2;                          % reflected radiation
      Q_r = nansum(Qr);
      if beta <= pi/2                                                    % this loop is for consideration of direct radiation,
          Q_sun(c) = (Q_sun+Q_D+Q_d+Q_r);
      else
          Q_sun = (Q_sun+Q_d+Q_r);
      end
  end
  disp(area)
  disp(Q_sun)