nested for loops and solar radiation

Hi I've got a function that works for part of my model. It works with the latitude 0 and 45 but not for 90 and gives me an error message of Attempted to access Q_sun(:,2); index out of bounds because numel(Q_sun)=1.

Also my code is giving me the wrong numerical results, if anyone has any thing that would help it is much appreciated. below is the function that I am running. Thank you.

function [ abs_sun_rad ] = working_sun_radiation(nc,lat,sec,day )
%calculates the absorbed solar radiation as a function of time of day and
%year, cloud  cover, and latitude
%   Inputs: 
    %day= day number
    %lat=latitude
    %sec=seconds passed from midnight
    %nc=cloud cover 0 or 100

%outputs: abs_rad= absorbtion solar radiation units of (W/m^2)
day=365;
hr=24;
abs_sun_rad= nan(day,hr); 

for day = 1:365;

       for hr=1:24; 
           sec= hr*3600;

           del = -23.45*(pi/180)*(cos((2*pi/365)*(day+9)));
              S=1357+45*(cos(((2*pi)/86400)*day));

          sin_h =[0,(sind(lat)*sind(del)+ cosd(lat)*cosd(del)*cosd((2*pi/86400)*(sec+43200)))];
          h=asin(sin_h);
          n=[1,1];
          m=(n/sin_h);
          Cext=0.128-(0.0235*(log(m)));

          i=((pi/2)-h);
          j=asin(0.75*(sin(i)));
          r=(0.5)*(abs((sin(i-j).^2)/(sin(i+j).^2)+(tan(i-j).^2/tan(i+j).^2)));

           for h = h

              if h>0 ;         
                  Insd=S*(exp(-Cext*m))*(sin_h*(1-0.71*nc));
              end
              if h<=0;
                  Insd=0;
              end

           end 

          Insg = 0.52*nc*Insd;
          Q_sun=(1-r)*Insd+(0.97*(Insg));
         abs_sun_rad(day,hr)=Q_sun(:,2);     

       end 
  end
  end

[Star Strider]

What does ‘the wrong numerical results’ mean? What are you expecting and what is your function returning?

You define ‘Q_sun’ as a scalar. What do you want it to be?

What happens at 90°? What variable does that calculation?

[Alyssa]

there should be zeros for the first 6 hours, because there is no sun light at that point. Im looking to have a matrix of 365x24 for. this is what happens at 90 degrees (lat), Attempted to access Q_sun(:,2); index out of bounds because numel(Q_sun)=1. thanks

Answer 1

[Star Strider]

Since ‘sin_h’ is a (1x2) vector, so are ‘Insd’, ‘Insg’ and Q_sun, if h>0.

If h<=0, all of these are scalars equal to zero.

See if changing:

Insd=0;

to:

Insd=[0 0];

improves things.