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 次查看（过去 30 天）
omkar bichkar
omkar bichkar 2015-6-16
评论： Ingrid 2015-6-16
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)

请先登录，再回答此问题。

回答（1 个）

Ingrid
Ingrid 2015-6-16
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 个评论
Ingrid
Ingrid 2015-6-16
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)

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Propagation and Channel Models 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by