solarradiation

 PUPROSE: Calculate solar radiation for a digital elevation model (DEM)
          over one year for clear sky conditions in W/m2
 -------------------------------------------------------------------
 USAGE: srad = solarradiation(dem,lat,cs,r)
 where: dem is the DEM to calculate hillshade for
        lat is the latitude vector for the DEM - same size as size(dem,2)
        cs is the cellsize in meters
        r is the ground reflectance (global value or map, default is 0.2)

       srad is the solar radiation in W/m2 over one year per grid cell

 EXAMPLE:
       srad = solarradiation(peaks(50)*100,54.9:-0.1:50,1000,0.2);
       - calculates the solar radiation for an example 50x50 peak surface.

 See also: GRADIENT, CART2POL

 Note: Follows the approach of Kumar et al 1997. Calculates clear sky
       radiation corrected for the incident angle (selfshading) plus
       diffuse and reflected radiation. Insolation is depending on time of year (and day),
       latitude, elevation, slope and aspect.
       Relief shading is not considered.
       Script uses simple unweighed gradient of 4 nearest neighbours for slope
       calculation.

 Reference: Kumar, L, Skidmore AK and Knowles E 1997: Modelling topographic variation in solar radiation in
            a GIS environment. Int.J.Geogr.Info.Sys. 11(5), 475-497


 Felix Hebeler, Dept. of Geography, University Zurich, May 2008.

0001 function srad = solarradiation(dem,lat,cs,r)
0002 % PUPROSE: Calculate solar radiation for a digital elevation model (DEM)
0003 %          over one year for clear sky conditions in W/m2
0004 % -------------------------------------------------------------------
0005 % USAGE: srad = solarradiation(dem,lat,cs,r)
0006 % where: dem is the DEM to calculate hillshade for
0007 %        lat is the latitude vector for the DEM - same size as size(dem,2)
0008 %        cs is the cellsize in meters
0009 %        r is the ground reflectance (global value or map, default is 0.2)
0010 %
0011 %       srad is the solar radiation in W/m2 over one year per grid cell
0012 %
0013 % EXAMPLE:
0014 %       srad = solarradiation(peaks(50)*100,54.9:-0.1:50,1000,0.2);
0015 %       - calculates the solar radiation for an example 50x50 peak surface.
0016 %
0017 % See also: GRADIENT, CART2POL
0018 %
0019 % Note: Follows the approach of Kumar et al 1997. Calculates clear sky
0020 %       radiation corrected for the incident angle (selfshading) plus
0021 %       diffuse and reflected radiation. Insolation is depending on time of year (and day),
0022 %       latitude, elevation, slope and aspect.
0023 %       Relief shading is not considered.
0024 %       Script uses simple unweighed gradient of 4 nearest neighbours for slope
0025 %       calculation.
0026 %
0027 % Reference: Kumar, L, Skidmore AK and Knowles E 1997: Modelling topographic variation in solar radiation in
0028 %            a GIS environment. Int.J.Geogr.Info.Sys. 11(5), 475-497
0029 %
0030 %
0031 % Felix Hebeler, Dept. of Geography, University Zurich, May 2008.
0032 
0033 
0034 %% parameters
0035 %It ;               % total hours of daily sunshine (calculated inline)
0036 %M ;                % air mass ratio parameter (calculated inline)
0037 %r = 0.20;          % ground reflectance coefficient (more sensible to give as input)
0038 %L=lat;             %latitude
0039 n = 1;              % timestep of calculation over sunshine hours: 1=hourly, 0.5=30min, 2=2hours etc
0040 tau_a    = 365;     %length of the year in days
0041 S0 = 1367;          % solar constant W m^-2   default 1367
0042 
0043 dr= 0.0174532925;   % degree to radians conversion factor
0044 
0045 %%  convert factors
0046 [slop,asp]=get_ders(dem,cs);   % calculate slope and aspect in radians using given cellsize cs
0047 [dummy,L]=meshgrid(1:size(dem,2),lat);   % grid latitude
0048 clear dummy;
0049 L=L*dr;                     % convert to radians
0050 fcirc = 360*dr; % 360 degrees in radians
0051 
0052 %% some setup calculations
0053 srad=0;
0054 sinL=sin(L);
0055 cosL=cos(L);
0056 tanL=tan(L);
0057 sinSlop=sin(slop);
0058 cosSlop=cos(slop);
0059 cosSlop2=cosSlop.*cosSlop;
0060 sinSlop2=sinSlop.*sinSlop;
0061 sinAsp=sin(asp);
0062 cosAsp=cos(asp);
0063 term1 = ( sinL.*cosSlop - cosL.*sinSlop.*cosAsp);
0064 term2 = ( cosL.*cosSlop + sinL.*sinSlop.*cosAsp);
0065 term3 = sinSlop.*sinAsp;
0066 %% loop over year
0067 for d = 1:tau_a; 
0068     %display(['Calculating melt for day ',num2str(d)])
0069     % clear sky solar radiation
0070     I0 = S0 * (1 + 0.0344*cos(fcirc*d/tau_a)); % extraterr rad per day
0071     % sun declination dS
0072     dS = 23.45 * dr* sin(fcirc * ( (284+d)/tau_a ) ); %in radians, correct/verified
0073     % angle at sunrise/sunset
0074     % t = 1:It; % sun hour
0075     hsr = real(acos(-tanL*tan(dS)));  % angle at sunrise
0076     % this only works for latitudes up to 66.5 deg N! Workaround:
0077     % hsr(hsr<-1)=acos(-1);
0078     % hsr(hsr>1)=acos(1);
0079     It=round(12*(1+mean(hsr(:))/pi)-12*(1-mean(hsr(:))/pi)); % calc daylength
0080 %%  daily loop
0081     I=0;
0082     for t=1:n:It % loop over sunshine hours
0083         % if accounting for shading should be included, calc hillshade here
0084         % hourangle of sun hs
0085         hs=hsr-(pi*t/It);               % hs(t)
0086         %solar angle and azimuth
0087         %alpha = asin(sinL*sin(dS)+cosL*cos(dS)*cos(hs));% solar altitude angle
0088         sinAlpha = sinL.*sin(dS)+cosL.*cos(dS).*cos(hs);
0089         %alpha_s = asin(cos(dS)*sin(hs)/cos(alpha)); % solar azimuth angle
0090         % correction  using atmospheric transmissivity taub_b
0091         M=sqrt(1229+((614.*sinAlpha)).^2)-614.*sinAlpha; % Air mass ratio
0092         tau_b = 0.56 * (exp(-0.65*M) + exp(-0.095*M));
0093         tau_d = 0.271-0.294*tau_b; % radiation diffusion coefficient for diffuse insolation
0094         tau_r = 0.271+0.706*tau_b; % reflectance transmitivity
0095         % correct for local incident angle
0096         cos_i = (sin(dS).*term1) + (cos(dS).*cos(hs).*term2) + (cos(dS).*term3.*sin(hs));
0097         Is = I0 * tau_b; % potential incoming shortwave radiation at surface normal (equator)
0098         % R = potential clear sky solar radiation W m2
0099         R = Is .* cos_i;
0100         R(R<0)=0;  % kick out negative values
0101         Id = I0 .* tau_d .* cosSlop2./ 2.*sinAlpha; %diffuse radiation;
0102         Ir = I0 .* r .* tau_r .* sinSlop2./ 2.* sinAlpha; % reflectance
0103         R= R + Id + Ir;
0104         R(R<0)=0; 
0105         I=I+R;% solar radiation per day (sunshine hours)
0106      end % end of sun hours in day loop
0107 %%  add up radiation part melt for every day
0108     srad = srad + I;
0109 end   % end of days in year loop
0110 
0111 
0112 %%
0113 function [grad,asp] = get_ders(dem,cs)
0114 % calculate slope and aspect (deg) using GRADIENT function
0115 [fx,fy] = gradient(dem,cs,cs); % uses simple, unweighted gradient of immediate neighbours
0116 [asp,grad]=cart2pol(fy,fx); % convert to carthesian coordinates
0117 grad=atan(grad); %steepest slope
0118 asp=asp.*-1+pi; % convert asp 0 facing south