Unsure on how to achieve goal - Function output isn't right

Hello,

I'm trying to work out what I'm doing wrong with my approach to a question. I am trying to create a function to calculate solar radiation absorption as a function of time of day&year, cloud cover, and latitude.

My code is as follows:

function [abs_sol_rad]=solarrad(nc,lat,t,I)
% Function to calculate absorbed solar radiation given cloud cover,
% latitude, seconds from midnight and day in the year. 
% Inputs:
    % I= day number (365 days)
    %lat = latitude
    %t = seconds from mightnight
    %nc = cloud cover from 0 (cloudless) to 100 (full cloud cover)
% Outputs are abs_sol_rad=absorption solar radiation with units W/m^2
d0=23.45; %angle between equatorial plane and earth rotational axis
wd=(2*pi)/86400;wy=(2*pi)/365; %wd= day rotational angle frequency, wy=year rotation angle frequency
I=365;
hr=24;
abs_sol_rad=nan(hr,I);
for I=1:365 %
    for hr=1:24
        t=hr*3600;
    end
    delta=d0*(pi/180)*(cos(wy*(I+9))); %declination of the sun from winter solstice (day 356)
    S=1357+45*(cos(wy*I)); %solar constant
    sin_h=max(0,(sind(lat)*sind(delta))+(cosd(lat)*cos(delta)*cos(wd*(t+43200)))); %sun height (angle) equation
    m=asin(sin_h); %optical length
    Cext=0.128-(0.0235.*log(m)); %integrated extinction coefficient
    for h=m
    i=((pi/2)-h); %angle of incidence
    j=asin(0.75*(sin(i))); %angle of reflection
    r=0.5*(abs((sin(i-j).^2)/(sin(i+j).^2)+(tan(i-j).^2/tan(i+j).^2))); %reflectance
    end
    if h>0
        Insd=S*(exp(-Cext*m))*(sin_h*(1-0.71*nc));
    end
    if h<=0
        Insd=0;
    end
 Insg=0.52*nc*Insd; %Global radiation
 Qsun=(Insd*(1-r))+(0.97*Insg); %absorbed radiation in water
 abs_sol_rad(hr,I)=Qsun;
end
end

The output I'm getting is just for the last time step I'm wanting (want a time step of 1 hour each day, 365 days). I need the output to be a matrix 24x365. Which I believe shows absorption when given latitude and cloud cover I guess.

Do I need to create a secondary script to use this function to create an array?

Should I just create a meshgrid to produce an 24x365 array and then just apply the function to each cell?

[Walter Roberson]

What is the point of multiplying hr*3600, 24 times in a row?

[CheeseToastie]

To get seconds for t later on.

Need it to run 24 times to get the hourly data, need t to be in seconds for further in the solar equation

Answer 1

[Miriam]

Have the second for loop end after assigning the value of abs_sol_rad, ie:

function [abs_sol_rad]=solarrad(nc,lat,t,I)

%% ...

for I=1:365 %

for hr=1:24

t=hr*3600;

delta=d0*(pi/180)*(cos(wy*(I+9)));

%% ...

abs_sol_rad(hr,I)=Qsun;

end

end

end

[CheeseToastie]

Did that and got an error, addedd a +1 to the abs_sol_rad(I,hr+1) and it's finally worked, I hope this makes sense, even if I don't undetstand why.

[Miriam]

Are you still using hr=1:24? I'm assuming you've adjusted your code slightly since your indexing has reversed. If you also, say, adjusted hr=0:23 that would explain why you need to add one, as MATLAB indexing begins from 1.