Commit e9d2307e authored by Jana Fischereit's avatar Jana Fischereit

Commit source code including license, readme and test

parents
This article is distributed under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/.
SURM
===============
Simple Urban Radiation Model is a model to calculate mean radiant temperature (Tmrt) in an idealised symmetric streat canyon.
The model can be used either stand-alone with standard meteorological parameters as inputs or for radiation modifications in a built-up area with radiation fluxes from a mesoscale model or measurements.
The model is written in Matlab but it can also be used in the freely available GNU Octave environment (www.gnu.org/software/octave/).
## Getting started
The program `model/get_started_example.m` provides an exampe, which variables have to be set for SURM and how to call the main routine.
`tests/plausibility_checks.m` provides another example for the plausibility tests described in Section 3.2 in Fischereit (2019).
Detailed description of variables and parameterisations can be found in Fischereit (2019) and in `model/calcTmrt`.
Example Linke turbidity factors of air mass 2 for Hamburg are located in `model/tlam2.txt`. For other locations those need to be adjusted as described in Fischereit (2019).
## Reference
Fischereit, J. (2019): The Simple Urban Radiation Model for estimating Mean Radiant
Temperature in idealized street canyons, Environmental Modelling and Software,
submitted
function[tmrt, varargout]=calcTmrt(h, w, ori,...
epsilon_w, epsilon_g, alpha_w, alpha_g, ...
T_a2m, e_2mPa,...
phi, jd, hh, ...
alpha_k, epsilon_p, pos_rel, ...
varargin)
% CALCTMRT Calculates Tmrt for an idealized symmetrical street
% canyon with radiation parametersiations for cloudless
% conditions
%
% REQUIRED INPUTS (in this order):
%
% h - Building height [m]
% w - Building width [m]
% ori - Orientation angle of the street with respect to north [°]
% epsilon_w - Emissivity of the wall
% epsilon_g - Emissivity of the ground
% alpha_w - Mean albedo of walls
% alpha_g - Mean albedo of ground
%
% T_a2m - Temperature in 2m [K]
% e_2mPa - Water vapour pressure in 2m [Pa]
%
% phi - Geographical latitude [°]
% jd - Julian day number
% hh - hour of interest
%
% alpha_k - short wave absorption coefficient of the human body [default: 0.7]
% epsilon_p - long wave absorption coefficient (=Emissivity) of human body [default: 0.97]
% pos_rel - Position of person (between 0 and 1 as seen
% from right wall looking into the direction of ori)
%
% OPTIONAL INPUTS (Needed only for specific parameterizations, input as
% parameter-value-pairs)
%
% T_w - Temperature of sunny wall (mean) [K]
% T_g - Temperature of sunny ground (mean) [K]
% v_10 - wind speed in 10 m, only used if
% basis_ts, to estimate wall and ground
% temperature [m/s]
% TLam2 - Linke turbidity factor of relative optical
% air mass 2 can be taken from
% http://www.soda-pro.com/web-services/atmosphere/linke-turbidity-factor-ozone-water-vapor-and-angstroembeta
% bo - Bowen ratio: 0.5 grass, 1.3 urban
% z_nn - Station height above sea level [m]
%
%
% Optional arguments in this order:
%
% basis_declin_calc- Type of declination calculation (valid are:
% 'MITRAS','vdi3789')
% basis_dir_rad - Define basis of short wave absorption
% applied here (valid are: 'MITRAS',
% 'vdi1994','esra')
% basis_diff_rad - Define basis of diffuse SW parameterization
% applied here (valid are: 'MITRAS',
% 'vdi1994','esra')
% basis_lw_sky - Define basis of long wave from sky
% (valid are: 'MITRAS', 'Idso1981', 'vdi3789')
% basis_ts - Handling of surface temperatures of ground and
% walls ('prescribed' (prescribed values for
% sun-lit wall and ground during the day and
% entire surfaces during the night) and
% 'prognostic' (prognostic calculation of ts)).
% basis_pvf - Define basis of view factor calculation for person
% (valid are: 'Fanger',
% 'Fangerfast','Cannistraro')
% l_dispValues - Logical: Assign calculated values of different
% terms to argout (true,false)
%
%
% OUTPUT: Optional outputs are returned, if l_dispValues=true
%
% tmrt - Mean radiant temperature at chosen point
% within street canyon in [K]
% is_pos_lit [optional] - Is position of person lit
% lw [optional] - longwave radiation component [W/m2]
% dif [optional] - diffusive radiation component [W/m2]
% dis [optional] - direct radiation component [W/m2]
% ges [optional] - total radiation [W/m2]
% T_g [optional] - Ground Temperature [K]
% T_w [optional] - Wall Temperature [K]
% Glob [optional] - Global radiation [W/m2]
% D_s [optional] - Diffuse radiation from sky [W/m2]
% E_s [optional] - Longwave radiation from sky [W/m2]
% R_sw [optional] - Direct radiation [W/m2]
% D_reflg [optional] - Reflected SW radiation at ground [W/m2]
% E_g [optional] - Longwave radiation from ground [W/m2]
% ill_wall [optional] - Illumunated wall ['l','r','n']
% f_p [optional] - Projection factor
% cos_eta [optional] - Cos(eta) angle for inclined surfaces
% alpha_sun [optional] - Elevation angle of the sun [rad]
% vfp2 [optional] - View factors of person to surrounding surfaces
% vf_shg [optional] - View factor of shaded ground to other surfaces
% vf_ilg [optional] - View factor of lit ground to other surfaces
% vf_ilw [optional] - View factor of lit wall to other surfaces
% vf_shwe [optional] - View factor of entirely shaded wall to other surfaces
%
% BACKGROUND:
%
% Tmrt is defined as the uniform temperature of a fictive black-body
% radiation enclosure in which a subject would experience the same
% net-radiation energy exchange as in the actual more complex
% radiation environment (according to Kantor, N., J. Unger, 2011)
%
% REFERENCES:
%
% Fischereit, J., 2019: The Simple Urban Radiation Model for
% estimating Mean Radiant Temperature in idealized
% street canyons, Environmental Software and
% Modelling, submitted
%
%% DEFAULT OPTIONS AND OPTIONAL INPUT ARGUMENTS
[T_w, T_g, v_10, TLam2, bo, z_nn, ...
basis_declin_calc, basis_dir_rad, basis_diff_rad, ...
basis_lw_sky, basis_ts, basis_pvf, l_dispValues] = se_parse_Input2(varargin);
%% CONSTANTS
p0 = 1013.25; % pressure at NN [hPa]
sigma = 5.67*10^(-8); % Boltzmann
pb = p0 * exp(-z_nn/8434.5); % derived pressure in station height [hPa]
I_infty = 1367; % solar constant
l = 'infinite'; % Length of street fixed length [m];
% currently not not fully implemented
h_man = 1.1; % height of the center point of a person
%
% set fixed width in case of open area for view factor calculation
%if h == 0
% w = 100;
%end
%
%
%% DIRECT RADIATION + STREET DESIGN PARAMETERS
% astronomic parameters related to sw radiation
[alpha_sun, psi_sun] = se_calc_astronomic(jd,hh,phi,basis_declin_calc);
% sw radiation reaching ground after absorption in atm
[R_sw,I_0] = se_calc_sw_ground(I_infty, alpha_sun,jd, pb, p0, TLam2,basis_dir_rad);
% calculation of fraction of lit and shaded part of street canyon
[frac_wall, frac_floor, is_pos_lit, cos_eta,ill_wall] ...
= se_calc_illuminated(h, w, ori, pos_rel, psi_sun, alpha_sun);
% Estimate surface-to-surface view factors
[vf_shg, vf_ilg, vf_ilw, vf_shw, vf_shwe] ...
= se_calc_vf_a2a(h,w,frac_floor,frac_wall);
% person to surface view factors
[vfp2] = se_calc_vf_p2a(h,h_man,w,l,pos_rel,frac_wall,frac_floor,ill_wall,basis_pvf);
% direct solar radiation reaching the human body
[I_t, f_p] = se_calc_dir_man(R_sw,alpha_sun, is_pos_lit);
%
%% DIFFUSE RADIATION
% diffuse solar radiation from sky
D_s = se_calc_diff_sky(R_sw,alpha_sun,pb,p0,TLam2,I_0,basis_diff_rad);
% total diffuse solar radiation reaching the human body:
% 1) sky
% 2) reflected (by walls and ground based on part that is lit)
[D_t] =...
se_calc_diff_man(D_s,R_sw,alpha_w,alpha_g,vfp2,...
vf_shg, vf_ilg, vf_ilw, vf_shw, vf_shwe,cos_eta);
%% LONGWAVE RADIATION
% longwave radiation from sky
E_s = se_calc_lon_sky(e_2mPa,T_a2m,sigma, basis_lw_sky);
% longwave radiation from ground and walls
[E_g, E_w,T_g,T_w]=se_calc_lon_morph(basis_ts,T_g,T_w,...
T_a2m,epsilon_g,epsilon_w,sigma,frac_wall,frac_floor,E_s,...
D_s,alpha_g,R_sw,cos_eta,ill_wall,...
alpha_w,v_10,bo, h,...
vf_shg, vf_ilg, vf_ilw, vf_shw, vf_shwe);
% total longwave radiation reaching the human body:
% 1) sky
% 2) ground
% 3) walls
% 4) (reflected at walls + ground) depends on used version
E_t = se_calc_lon_man(E_s,E_g,E_w,epsilon_w,epsilon_g,vfp2,...
vf_shg, vf_ilg, vf_ilw, vf_shw, vf_shwe);
%
%% T_MRT [K]
tmrt = (1/sigma*(E_t + alpha_k/epsilon_p*(D_t + I_t))).^0.25;
%% IO
varargout{1}=is_pos_lit;
if l_dispValues
lw = E_t;
dfs = alpha_k/epsilon_p*(D_t);
dis = alpha_k/epsilon_p*(I_t);
ges = E_t+alpha_k/epsilon_p*(D_t+I_t);
varargout{2}=lw;
varargout{3}=dfs;
varargout{4}=dis;
varargout{5}=ges;
varargout{6}=T_g;
varargout{7}=T_w;
%glob
varargout{8}=is_pos_lit*R_sw*cos_eta(2)+D_s;
% diffuse radiation from sky
varargout{9}=D_s;
% longwave downward radiation
varargout{10}=E_s;
% direct radiation
varargout{11}=R_sw;
% SW from ground
varargout{12}= varargout{8}*alpha_g;
% LW from ground
varargout{13}=E_g;
varargout{14}=ill_wall;
varargout{15}=f_p;
varargout{16}=cos_eta;
varargout{17}=alpha_sun;
varargout{18}=vfp2;
varargout{19}=vf_shg;
varargout{20}=vf_ilg;
varargout{21}=vf_ilw;
varargout{22}=vf_shw;
varargout{23}=vf_shwe;
end
%
end
%% Example program to set-up a Tmrt calculation with SURM
%
% This program shows how to set different input variables for SURM
% and how to call the main SURM calculation
%
% For details on the calculation of Tmrt and the input parameter see
% the main of SURM: calcTmrt
%
%% USER INPUT SECTION
clear all
%
%------- street design ------
%
h = 10; % building height [m]
w = 10; % street width [m]
ori = 57; % steet angle towards north [°]
%
%------- Built-up environment properties ------
%
T_w_init = 293.15; % wall temperature [K]
T_g_init = 293.15; % ground temperature [K]
epsilon_w = 0.9; % emissivity of wall
epsilon_g = 0.95; % emissivity of ground
alpha_w = 0.15; % mean wall albedo
alpha_g = 0.15; % mean ground albedo
%
%------- Metorological input --------
%
T_a2m = 20 + 273.15; % air temperature in 2m [K]
RH_a2m = 50; % Relative humidity in 2m [%]
z_nn = 11; % Station height [m]; e.g. Hamburg, Fuhlsbuettel = 11 m
%
%------- Radiation input --------
%
phi = 53; % latitude [°]
year = 2015; % ATTENTION: No leap years!
month = 6; % Month to simulate
day = 21; % Day to simulate
hour = 12; % Hour to simulate
fileID = fopen('tlam2.txt'); % Linke turbidity factor taken from http://www.soda-pro.com/web-services/atmosphere/linke-turbidity-factor-ozone-water-vapor-and-angstroembeta
tlam2_cell= textscan(fileID,'%4s %f','HeaderLines',7,'Delimiter',' ');
%
%------ human characteristics input ------
%
alpha_k = 0.7; % absorption coefficient of the irradiated body surface
% for short-wave radiation (default: 0.7)
epsilon_p = 0.97; % absorption coefficient for long-wave radiation (= emissivity) of the human body (default: 0.97)
pos_rel = 0.5; % position within the street canyon as seen from right to
% left in a N-S street canyon, valid are values
% between 0 and 1
%
%------ Optional inputs ------
%
% For a description of optional inputs see calcTmrt
%
%% ############# END OF USER INPUT SECTION ##############################
%
% conversions
% 1) julian day in radian
jd = datenum(year,month,day) - datenum(year, 1, 1) + 1;
%
% 2) RH in e using magnus equation
e_2mPa = ...
(RH_a2m*6.112*exp(17.62*(T_a2m-273.15)/(243.12+(T_a2m-273.15))))...
/100*100;
%
%
% calculation of tmrt
tic % to measure calculation time
tmrt = calcTmrt(h,w,ori,...
epsilon_w,epsilon_g,alpha_w,alpha_g,...
T_a2m,e_2mPa,...
phi,jd,hour,...
alpha_k,epsilon_p,pos_rel,...
'T_w',T_w_init,'T_g',T_g_init,...
'TLam2',tlam2_cell{1,2}(month),'z_nn',z_nn)
toc
function [T_s] = se_calc_Ts(E_swnet, epsilon, ...
E_lwin, sigma, T_a2m, bo, v_10, surface_flag)
% SE_CALC_TW Calculates wall/ground Temperature with Newton
% approximation assuming immidiate equilibrium
%
% INPUT:
%
% E_swnet - SW net radiation balance at surface [Wm^-2]
% epsilon - emissivity
% E_lwnet - Lw incoming at surface [Wm^-2]
% sigma - Stefan-Boltzmann-constant [Wm^-2K^-4]
% T_a2m - Air Temperature [K]
% bo - Bowen ratio
% v_10 - wind speed in 10 m [m/s]
% surface_flag - Flag: 'w' (wall), 'g' (ground)
%
%
% OUTPUT:
%
% T_s - Surface temperature [K]
%
% REFERENCES
%
% Gierisch, A., 2011: Mikroskalige Modellierung meteorologischer und
% anthropogener Einflüsse auf die Wä̈rmeabgabe eines Gebäudes
% Master’s thesis, Geowissenschaften, University of Hamburg,
% Hamburg, Germany, p. 12
%
% Jendritzky, G., G. Menz, W. Schmidt-Kessen, H. Schirmer, 1990:
% Methodik zur raumbezogenen Bewertung der thermischen Komponente
% des Bioklima des Menschen (Fortgeschriebenes Klima-Michel-Modell).
% Technical Report 114, Akademie für Raumforschung und
% Landesplanung, Hannover.
%
% Staiger, H., 2014: Die Strahlungskomponente im thermischen
% Wirkungskomplex für operationelle Anwendungen in der Human-
% Biometeorologie Ph.D. thesis, Geowissenschaften,
% Albert-Ludwigs-Universit ̈at Freiburg im Breisgau, Germany.
%
%% initalize
T_s = T_a2m;
%
%% Newton approximation
if strcmp(surface_flag, 'w')
% PARAMETERS
% heat transfer between interior and outside
U = 1;
C = 1/(1/U-0.04);
% indoor room temperature
T_in = 20+273.15;
% heat transfer coefficient
h_aw = 4 + 4 * v_10; % wall
% CALCULATION
for i = 1:3 % three iteration steps enough according to Staiger (2014)
T_s = T_s - ...
((E_swnet + ...
epsilon * (E_lwin - sigma * T_s^4)) ...
- C * (T_in - T_s) ...
+ h_aw * (T_a2m - T_s))...
/...
(-4 * epsilon * sigma * T_s^3 ...
+ C - h_aw);
end
else
% PARAMETERS
% heat transfer coefficient
h_ag = 6.2 + 4.26 * v_10; % ground
% parameter of soil characteristics
if E_swnet > 0 %during the day
b_terr = -0.19;
else
b_terr = -0.32;
end
% CALCULATION
for i = 1:3 % three iteration steps enough according to Staiger (2014)
T_s = T_s - ...
((E_swnet ...
+ epsilon * (E_lwin - sigma * T_s^4))...
* (1 + b_terr) + ...
+ (1 + 1/bo) * h_ag * (T_a2m - T_s))...
/...
(-4 * epsilon * sigma * (1 + b_terr) * T_s^3 ...
- (1 + 1/bo) * h_ag);
end
end
end
function [alpha_sun, psi_sun, varargout] = se_calc_astronomic(...
jd,hh,phi,basis_declin_calc)
% SE_CALC_ASTRONOMIC Calculates astronomic parameters of direct
% radiation
%
% INPUT:
%
% jd - Julian day number
% hh - hour of interest
% phi - Geographical latitude [°]
% basis_declin_calc - Calculation Method for declination, valid are
% {MITRAS,vdi3789}
%
%
% OUTPUT:
%
% alpha_sun - Elevation angle of the sun [rad]
% psi_sun - Azimuth angle of the sun [rad]
% cos_zenit_sun - Cosine of zenith angle of the sun []
%
% REFERENCES
%
% Schluenzen, K. H., D.D. Flagg, B.H. Fock, A. Gierisch, C. Lüpkes,
% V. Reinhardt, C. Spensberger, 2012: Scientific Documentation of
% the Multiscale Model System M-SYS, Meteorological Insitute,
% University of Hamburg, MEMI Technical Report 4, P. 94
%
% VDI 3789, Part 3 (2001): Environmental meteorology - Interactions
% between atmosphere and surfaces: Calculation of spectral
% irradiances in the solar wavelength range
%
%
%% DEFINITION OF D2R
% because: matlab intrinsic function deg2rad is slow.
d2r = (pi()/180);
r2d = (180/pi());
%% CALCULATION OF HOUR ANGLE FOR specific hour
omega = ((hh - 12) * 15); % from VDI3789, Part 3 p 54
%% CALCULATION OF EARTH'S DECLINATION
switch basis_declin_calc
case lower('MITRAS')
jd_rad = 2*pi*(jd-1)/365;
delta = 0.006918 ...
- 0.399912 * cos(jd_rad) + 0.070257 * sin(jd_rad)...
- 0.006758 * cos(2*jd_rad) + 0.000907 * sin(2*jd_rad) ...
- 0.002697 * cos(3*jd_rad) + 0.001480 * sin(3*jd_rad);
case lower('vdi3789')
delta = (asin(0.3978 * sin(d2r*((0.9856 * jd - 2.72) - 77.51...
+ 1.92 * sin(d2r*(0.9856 * jd - 2.72))))));
otherwise
error('You used an unknown value for declin_calc')
end
%% CALCULATION OF ELEVATION ANGLE
sinalpha = sin(delta).*sin(d2r*phi) + cos(delta).*cos(d2r*phi).*cos(d2r*omega);
%% CALCULATION OF AZIMUTH ANGLE
%
% here pi/2-zenit is used for alpha_sun, because alpha_sun is set to 0
% during the night already
cos_psi = (sinalpha.*sin(d2r*phi)-sin(delta)) ...
./ (cos(asin(sinalpha)).*cos(d2r*phi));
% for the correct sign of psi
psi_sun = acos(cos_psi).*sign(omega);
%% ADJUSTMENT OF ELEVATION ANGLE DURING THE NIGHT
alpha_sun = asin(sinalpha);
%% OPTIONAL OUTPUT: TIME OF SUNRISE
hh_sr = (((acos(sin(1*pi)-sin(delta).*sin(d2r*phi)/(cos(delta).*cos(d2r*phi))))*r2d)/15+12-24)*-1;
hh_ss = ((acos(-sin(delta).*sin(d2r*phi)/(cos(delta).*cos(d2r*phi))))*r2d)/15+12;
varargout{1} = hh_sr;
varargout{2} = hh_ss;
end
function[D_t]=se_calc_diff_man(D_s,R_sw,alpha_w,alpha_g,vfp2,...
vf_shag, vf_litg, vf_litw, vf_shaw, vf_shawe,cos_eta)
% SE_CALC_DIFF_MAN Calculates total diffuse radiation
%
% INPUT:
%
% D_s - Diffuse SW radiation from sky (actually only
% isotropic component)
% R_sw - SW dirct radiation after absorption in atm +
% corrected for sun height
% alpha_w - Mean albedo of walls
% alpha_g - Mean albedo of ground
% vfp2 - view factor of person to surfaces
% vf_shag - View Factors for shaded ground to other
% surfaces
% vf_litg - View Factors for lit ground to other surfaces
% vf_litw - View Factors for lit wall to other surfaces
% vf_shaw - View Factors for shaded wall to other surfaces
% vf_shawe - View Factors for entirley shaded wall to other
% surfaces
% cos_eta - cos(eta) angle for inclined surfaces
%
%
% OUTPUT:
%
% D_t - total diffuse Radiation reaching the human body
%
%
%% CALCUALTION
%
if isnan(cos_eta(1))
sw_dir_refl_w = 0; % sun parallel to canyon
else
sw_dir_refl_w = alpha_w * R_sw * cos_eta(1);
end
D_t = vfp2.w.lit * sw_dir_refl_w ...
+ vfp2.w.sha * vf_litw.s * alpha_w * D_s ...
+ vfp2.w.sha * vf_shaw.s * alpha_w * D_s ...
+ vfp2.w.shae* vf_shawe.s * alpha_w * D_s ...
+ vfp2.g.lit * alpha_g * R_sw * cos_eta(2) ...
+ vfp2.g.sha * vf_shag.s * alpha_g * D_s ......
+ vfp2.g.lit * vf_litg.s * alpha_g * D_s ...
+ vfp2.s * D_s;
%
%
end
function[D_s]=se_calc_diff_sky(R_sw,alpha_sun,p,p0,TLam2,I_0,basis_diff_rad)
% SE_CALC_DIFF_SKY Calculates diffuse radiaton from sky
%
% INPUT:
%
% R_sw - SW from sky reaching ground without
% consideration of city's morphology for
% different times
% alpha_sun - Elevation angle of the sun [rad]
% p - Pressure
% p0 - Mean sea level pressure
% TLam2 - Linke turbidity factor of relative optical
% air mass 2
% I_0 - solar constant with accounted annual cycle
% basis_diff_rad - Define basis of diffuse SW parameterization
% applied here (valid are: 'MITRAS',
% 'Bruse1999', 'vdi1994', 'esra')
%
%
% OUTPUT:
%
% D_s - Diffuse radiation from sky
%
% REFERENCES
%
% Bruse, M., 1999: Die Auswirkungen kleinskaliger Umweltgestaltung
% auf das Mikroklima – Entwicklung des prognostischen numerischen
% Modells ENVI-met zur Simulation der Wind-, Temperatur- und
% Feuchteverteilung in städtischen Strukturen Ph.D. thesis,
% Geowissenschaften, Ruhr-Universität Bochum, Bochum, Germany.
% p. 53
%
% Gierisch, A., 2011: Mikroskalige Modellierung meteorologischer und
% anthropogener Einflüsse auf die Wä̈rmeabgabe eines Gebäudes
% Master’s thesis, Geowissenschaften, University of Hamburg,
% Hamburg, Germany, p. 12
%
% Rigollier C, Bauer O, Wald L (2000): On the clear sky model of the
% ESRA - European Solar Radiation Atlas - with respect to the
% HELIOSAT method. Solar Energy, 68 (1): 33 - 48
%
% VDI 3789, 2 (1994): Interactions between Atmosphere and Surfaces
% Calculation of Short-wave and Long-wave Radiation
%
%%
switch basis_diff_rad
case lower('MITRAS')
% Derive diffuse part from direct shortwave radiation radiation
gamma = 1./(1+8*(sin(alpha_sun)).^0.7);
%
% warning is issued for high latitudes, when division by 0
% Comment: sin(alpha_sun) is here part of the parameterization to derive
% the amount of diffuse radiation from direct radiation and does not
% depend on the elevation angle of the surface
D_s = R_sw.*sin(alpha_sun).*(gamma/(1-gamma));
case lower('vdi1994')
% Derive diffuse part as difference from global radiation
TL = TLam2*p/p0;
G = 0.84 * I_0 * sin(alpha_sun) * exp(-0.027*TL/sin(alpha_sun));
D_s = G - R_sw * sin(alpha_sun);
case lower('esra')
% Estimate clear sky independent of direct radiation
TL = TLam2*p/p0;
%
Trd = -0.015843 + 0.030543 * TL + 0.0003797 * TL^2;
%
a0 = 0.26463-0.061581*TL+0.0031408*TL^2;
if a0*Trd < 0.002
a0 = 0.002/Trd;
end
a1 = 2.04020 + 0.018945 * TL - 0.0111610 * TL^2;
a2 = -1.33025 + 0.0392311 * TL + 0.0085079 * TL^2;
F_d = a0 + a1*sin(alpha_sun)+a2*(sin(alpha_sun))^2;
%
D_s = I_0*Trd*F_d;
end
if D_s < 0 || isnan(D_s)
D_s = 0;
end
end
function[I_t,f_p]=se_calc_dir_man(R_sw,alpha_sun,is_pos_lit)
% SE_CALC_DIR_MAN Calculates SW-radiation, which reaches
% the human body using a projection factor
% (= weighting function)
%
% INPUT: