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.

1 view (last 30 days)
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)

Answers (1)

Ingrid
Ingrid on 16 Jun 2015
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
  1 Comment
Ingrid
Ingrid on 16 Jun 2015
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)

Sign in to comment.

Categories

Find more on Get Started with MATLAB in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by