1 |
greg |
2.1 |
#ifndef lint |
2 |
greg |
2.11 |
static const char RCSid[] = "$Id: gendaylit.c,v 2.10 2013/02/13 18:30:22 greg Exp $"; |
3 |
greg |
2.1 |
#endif |
4 |
greg |
2.9 |
/* Copyright (c) 1994,2006 *Fraunhofer Institut for Solar Energy Systems |
5 |
|
|
* Heidenhofstr. 2, D-79110 Freiburg, Germany |
6 |
greg |
2.1 |
* *Agence de l'Environnement et de la Maitrise de l'Energie |
7 |
|
|
* Centre de Valbonne, 500 route des Lucioles, 06565 Sophia Antipolis Cedex, France |
8 |
|
|
* *BOUYGUES |
9 |
|
|
* 1 Avenue Eugene Freyssinet, Saint-Quentin-Yvelines, France |
10 |
greg |
2.9 |
* |
11 |
|
|
* 24.1.2006 some adjustments for cygwin compilation, inclusion of RADIANCE3.7 libraries, by J. Wienold |
12 |
|
|
* 2011/10/08 [email protected]: |
13 |
|
|
* - integrated coeff_perez.dat and defangles.dat |
14 |
|
|
* - avoid some segfaults caused by out of range parameters and |
15 |
|
|
* - numerically dangerous range checks |
16 |
greg |
2.1 |
*/ |
17 |
|
|
|
18 |
|
|
/* |
19 |
|
|
* gendaylit.c program to generate the angular distribution of the daylight. |
20 |
|
|
* Our zenith is along the Z-axis, the X-axis |
21 |
|
|
* points east, and the Y-axis points north. |
22 |
|
|
*/ |
23 |
|
|
|
24 |
greg |
2.10 |
#define _USE_MATH_DEFINES |
25 |
|
|
|
26 |
greg |
2.1 |
#include <stdio.h> |
27 |
|
|
#include <string.h> |
28 |
|
|
#include <math.h> |
29 |
|
|
#include <stdlib.h> |
30 |
|
|
|
31 |
|
|
#include "color.h" |
32 |
|
|
#include "paths.h" |
33 |
|
|
|
34 |
greg |
2.9 |
#define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]) |
35 |
greg |
2.1 |
|
36 |
|
|
double normsc(); |
37 |
|
|
|
38 |
greg |
2.11 |
/*static char *rcsid="$Header: /cvs/radiance/ray/src/gen/gendaylit.c,v 2.10 2013/02/13 18:30:22 greg Exp $";*/ |
39 |
greg |
2.9 |
|
40 |
|
|
float coeff_perez[] = { |
41 |
|
|
1.3525,-0.2576,-0.2690,-1.4366,-0.7670,0.0007,1.2734,-0.1233,2.8000,0.6004,1.2375,1.000,1.8734,0.6297,0.9738,0.2809,0.0356,-0.1246,-0.5718,0.9938, |
42 |
|
|
-1.2219,-0.7730,1.4148,1.1016,-0.2054,0.0367,-3.9128,0.9156,6.9750,0.1774,6.4477,-0.1239,-1.5798,-0.5081,-1.7812,0.1080,0.2624,0.0672,-0.2190,-0.4285, |
43 |
|
|
-1.1000,-0.2515,0.8952,0.0156,0.2782,-0.1812,-4.5000,1.1766,24.7219,-13.0812,-37.7000,34.8438,-5.0000,1.5218,3.9229,-2.6204,-0.0156,0.1597,0.4199,-0.5562, |
44 |
|
|
-0.5484,-0.6654,-0.2672,0.7117,0.7234,-0.6219,-5.6812,2.6297,33.3389,-18.3000,-62.2500,52.0781,-3.5000,0.0016,1.1477,0.1062,0.4659,-0.3296,-0.0876,-0.0329, |
45 |
|
|
-0.6000,-0.3566,-2.5000,2.3250,0.2937,0.0496,-5.6812,1.8415,21.0000,-4.7656,-21.5906,7.2492,-3.5000,-0.1554,1.4062,0.3988,0.0032,0.0766,-0.0656,-0.1294, |
46 |
|
|
-1.0156,-0.3670,1.0078,1.4051,0.2875,-0.5328,-3.8500,3.3750,14.0000,-0.9999,-7.1406,7.5469,-3.4000,-0.1078,-1.0750,1.5702,-0.0672,0.4016,0.3017,-0.4844, |
47 |
|
|
-1.0000,0.0211,0.5025,-0.5119,-0.3000,0.1922,0.7023,-1.6317,19.0000,-5.0000,1.2438,-1.9094,-4.0000,0.0250,0.3844,0.2656,1.0468,-0.3788,-2.4517,1.4656, |
48 |
|
|
-1.0500,0.0289,0.4260,0.3590,-0.3250,0.1156,0.7781,0.0025,31.0625,-14.5000,-46.1148,55.3750,-7.2312,0.4050,13.3500,0.6234,1.5000,-0.6426,1.8564,0.5636}; |
49 |
|
|
|
50 |
|
|
|
51 |
|
|
float defangle_theta[] = {84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 12, 12, 12, 12, 12, 12, 0}; |
52 |
|
|
|
53 |
|
|
float defangle_phi[] = {0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 0, 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 0, 60, 120, 180, 240, 300, 0}; |
54 |
greg |
2.1 |
|
55 |
|
|
|
56 |
|
|
|
57 |
|
|
/* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */ |
58 |
|
|
double sky_brightness(); |
59 |
|
|
double sky_clearness(); |
60 |
|
|
|
61 |
|
|
/* calculation of the direct and diffuse components from the Perez parametrization */ |
62 |
greg |
2.9 |
double diffuse_irradiance_from_sky_brightness(); |
63 |
greg |
2.1 |
double direct_irradiance_from_sky_clearness(); |
64 |
|
|
|
65 |
|
|
|
66 |
|
|
/* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */ |
67 |
|
|
double glob_h_effi_PEREZ(); |
68 |
|
|
double glob_h_diffuse_effi_PEREZ(); |
69 |
|
|
double direct_n_effi_PEREZ(); |
70 |
|
|
/*likelihood check of the epsilon, delta, direct and diffuse components*/ |
71 |
|
|
void check_parametrization(); |
72 |
|
|
void check_irradiances(); |
73 |
|
|
void check_illuminances(); |
74 |
|
|
void illu_to_irra_index(); |
75 |
greg |
2.9 |
void print_error_sky(); |
76 |
greg |
2.1 |
|
77 |
|
|
|
78 |
|
|
/* Perez sky luminance model */ |
79 |
|
|
double calc_rel_lum_perez(double dzeta,double gamma,double Z, |
80 |
greg |
2.9 |
double epsilon,double Delta,float coeff_perez[]); |
81 |
greg |
2.1 |
/* coefficients for the sky luminance perez model */ |
82 |
greg |
2.9 |
void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]); |
83 |
greg |
2.1 |
double radians(double degres); |
84 |
|
|
double degres(double radians); |
85 |
|
|
void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z); |
86 |
|
|
double integ_lv(float *lv,float *theta); |
87 |
|
|
|
88 |
greg |
2.9 |
void printdefaults(); |
89 |
|
|
void computesky(); |
90 |
|
|
void printhead(int ac, char** av); |
91 |
|
|
void userror(char* msg); |
92 |
|
|
void printsky(); |
93 |
greg |
2.1 |
|
94 |
greg |
2.9 |
FILE * frlibopen(char* fname); |
95 |
greg |
2.1 |
|
96 |
|
|
/* astronomy and geometry*/ |
97 |
|
|
double get_eccentricity(); |
98 |
|
|
double air_mass(); |
99 |
|
|
|
100 |
greg |
2.9 |
extern int jdate(int month, int day); |
101 |
|
|
extern double stadj(int jd); |
102 |
|
|
extern double sdec(int jd); |
103 |
|
|
extern double salt(double sd, double st); |
104 |
|
|
extern double sazi(double sd, double st); |
105 |
greg |
2.1 |
|
106 |
|
|
|
107 |
|
|
/* sun calculation constants */ |
108 |
|
|
extern double s_latitude; |
109 |
|
|
extern double s_longitude; |
110 |
|
|
extern double s_meridian; |
111 |
|
|
|
112 |
|
|
const double AU = 149597890E3; |
113 |
|
|
const double solar_constant_e = 1367; /* solar constant W/m^2 */ |
114 |
|
|
const double solar_constant_l = 127.5; /* solar constant klux */ |
115 |
|
|
|
116 |
|
|
const double half_sun_angle = 0.2665; |
117 |
|
|
const double half_direct_angle = 2.85; |
118 |
|
|
|
119 |
|
|
const double skyclearinf = 1.000; /* limitations for the variation of the Perez parameters */ |
120 |
|
|
const double skyclearsup = 12.1; |
121 |
|
|
const double skybriginf = 0.01; |
122 |
|
|
const double skybrigsup = 0.6; |
123 |
|
|
|
124 |
|
|
|
125 |
|
|
|
126 |
|
|
/* required values */ |
127 |
|
|
int month, day; /* date */ |
128 |
|
|
double hour; /* time */ |
129 |
|
|
int tsolar; /* 0=standard, 1=solar */ |
130 |
|
|
double altitude, azimuth; /* or solar angles */ |
131 |
|
|
|
132 |
|
|
|
133 |
|
|
|
134 |
|
|
/* definition of the sky conditions through the Perez parametrization */ |
135 |
greg |
2.9 |
double skyclearness = 0; |
136 |
|
|
double skybrightness = 0; |
137 |
greg |
2.1 |
double solarradiance; /*radiance of the sun disk and of the circumsolar area*/ |
138 |
greg |
2.9 |
double diffuseilluminance, directilluminance, diffuseirradiance, directirradiance; |
139 |
greg |
2.1 |
double sunzenith, daynumber=150, atm_preci_water=2; |
140 |
|
|
|
141 |
greg |
2.9 |
double sunaltitude_border = 0; |
142 |
|
|
double diffnormalization = 0; |
143 |
|
|
double dirnormalization = 0; |
144 |
greg |
2.1 |
double *c_perez; |
145 |
|
|
|
146 |
|
|
int output=0; /*define the unit of the output (sky luminance or radiance): visible watt=0, solar watt=1, lumen=2*/ |
147 |
|
|
int input=0; /*define the input for the calulation*/ |
148 |
|
|
|
149 |
greg |
2.9 |
int suppress_warnings=0; |
150 |
|
|
|
151 |
greg |
2.1 |
/* default values */ |
152 |
|
|
int cloudy = 0; /* 1=standard, 2=uniform */ |
153 |
|
|
int dosun = 1; |
154 |
|
|
double zenithbr = -1.0; |
155 |
|
|
double betaturbidity = 0.1; |
156 |
|
|
double gprefl = 0.2; |
157 |
|
|
int S_INTER=0; |
158 |
|
|
|
159 |
|
|
/* computed values */ |
160 |
|
|
double sundir[3]; |
161 |
greg |
2.9 |
double groundbr = 0; |
162 |
greg |
2.1 |
double F2; |
163 |
|
|
double solarbr = 0.0; |
164 |
|
|
int u_solar = 0; /* -1=irradiance, 1=radiance */ |
165 |
|
|
|
166 |
|
|
char *progname; |
167 |
|
|
char errmsg[128]; |
168 |
|
|
|
169 |
|
|
|
170 |
greg |
2.9 |
int main(int argc, char** argv) |
171 |
greg |
2.1 |
{ |
172 |
|
|
int i; |
173 |
|
|
|
174 |
|
|
progname = argv[0]; |
175 |
|
|
if (argc == 2 && !strcmp(argv[1], "-defaults")) { |
176 |
|
|
printdefaults(); |
177 |
greg |
2.9 |
return 0; |
178 |
greg |
2.1 |
} |
179 |
|
|
if (argc < 4) |
180 |
|
|
userror("arg count"); |
181 |
|
|
if (!strcmp(argv[1], "-ang")) { |
182 |
greg |
2.9 |
altitude = atof(argv[2]) * (M_PI/180); |
183 |
|
|
azimuth = atof(argv[3]) * (M_PI/180); |
184 |
greg |
2.1 |
month = 0; |
185 |
|
|
} else { |
186 |
|
|
month = atoi(argv[1]); |
187 |
|
|
if (month < 1 || month > 12) |
188 |
|
|
userror("bad month"); |
189 |
|
|
day = atoi(argv[2]); |
190 |
|
|
if (day < 1 || day > 31) |
191 |
|
|
userror("bad day"); |
192 |
|
|
hour = atof(argv[3]); |
193 |
|
|
if (hour < 0 || hour >= 24) |
194 |
|
|
userror("bad hour"); |
195 |
|
|
tsolar = argv[3][0] == '+'; |
196 |
|
|
} |
197 |
|
|
for (i = 4; i < argc; i++) |
198 |
|
|
if (argv[i][0] == '-' || argv[i][0] == '+') |
199 |
|
|
switch (argv[i][1]) { |
200 |
|
|
case 's': |
201 |
|
|
cloudy = 0; |
202 |
|
|
dosun = argv[i][0] == '+'; |
203 |
|
|
break; |
204 |
|
|
case 'r': |
205 |
|
|
case 'R': |
206 |
|
|
u_solar = argv[i][1] == 'R' ? -1 : 1; |
207 |
|
|
solarbr = atof(argv[++i]); |
208 |
|
|
break; |
209 |
|
|
case 'c': |
210 |
|
|
cloudy = argv[i][0] == '+' ? 2 : 1; |
211 |
|
|
dosun = 0; |
212 |
|
|
break; |
213 |
|
|
case 't': |
214 |
|
|
betaturbidity = atof(argv[++i]); |
215 |
|
|
break; |
216 |
greg |
2.9 |
case 'w': |
217 |
|
|
suppress_warnings = 1; |
218 |
|
|
break; |
219 |
greg |
2.1 |
case 'b': |
220 |
|
|
zenithbr = atof(argv[++i]); |
221 |
|
|
break; |
222 |
|
|
case 'g': |
223 |
|
|
gprefl = atof(argv[++i]); |
224 |
|
|
break; |
225 |
|
|
case 'a': |
226 |
greg |
2.9 |
s_latitude = atof(argv[++i]) * (M_PI/180); |
227 |
greg |
2.1 |
break; |
228 |
|
|
case 'o': |
229 |
greg |
2.9 |
s_longitude = atof(argv[++i]) * (M_PI/180); |
230 |
greg |
2.1 |
break; |
231 |
|
|
case 'm': |
232 |
greg |
2.9 |
s_meridian = atof(argv[++i]) * (M_PI/180); |
233 |
greg |
2.1 |
break; |
234 |
|
|
|
235 |
|
|
case 'O': |
236 |
greg |
2.11 |
output = atoi(argv[++i]); /*define the unit of the output of the program : |
237 |
greg |
2.9 |
sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm) */ |
238 |
greg |
2.1 |
break; |
239 |
|
|
|
240 |
|
|
case 'P': |
241 |
|
|
input = 0; /* Perez parameters: epsilon, delta */ |
242 |
|
|
skyclearness = atof(argv[++i]); |
243 |
|
|
skybrightness = atof(argv[++i]); |
244 |
|
|
break; |
245 |
|
|
|
246 |
|
|
case 'W': /* direct normal Irradiance [W/m^2] */ |
247 |
|
|
input = 1; /* diffuse horizontal Irrad. [W/m^2] */ |
248 |
|
|
directirradiance = atof(argv[++i]); |
249 |
greg |
2.9 |
diffuseirradiance = atof(argv[++i]); |
250 |
greg |
2.1 |
break; |
251 |
|
|
|
252 |
|
|
case 'L': /* direct normal Illuminance [Lux] */ |
253 |
|
|
input = 2; /* diffuse horizontal Ill. [Lux] */ |
254 |
|
|
directilluminance = atof(argv[++i]); |
255 |
greg |
2.9 |
diffuseilluminance = atof(argv[++i]); |
256 |
greg |
2.1 |
break; |
257 |
|
|
|
258 |
|
|
case 'G': /* direct horizontal Irradiance [W/m^2] */ |
259 |
|
|
input = 3; /* diffuse horizontal Irrad. [W/m^2] */ |
260 |
|
|
directirradiance = atof(argv[++i]); |
261 |
greg |
2.9 |
diffuseirradiance = atof(argv[++i]); |
262 |
|
|
break; |
263 |
|
|
|
264 |
|
|
case 'l': |
265 |
|
|
sunaltitude_border = atof(argv[++i]); |
266 |
greg |
2.1 |
break; |
267 |
greg |
2.9 |
|
268 |
greg |
2.1 |
|
269 |
|
|
default: |
270 |
|
|
sprintf(errmsg, "unknown option: %s", argv[i]); |
271 |
|
|
userror(errmsg); |
272 |
|
|
} |
273 |
|
|
else |
274 |
|
|
userror("bad option"); |
275 |
|
|
|
276 |
greg |
2.9 |
if (fabs(s_meridian-s_longitude) > 30*M_PI/180) |
277 |
greg |
2.1 |
fprintf(stderr, |
278 |
|
|
"%s: warning: %.1f hours btwn. standard meridian and longitude\n", |
279 |
greg |
2.9 |
progname, (s_longitude-s_meridian)*12/M_PI); |
280 |
greg |
2.1 |
|
281 |
|
|
|
282 |
|
|
/* allocation dynamique de memoire pour les pointeurs */ |
283 |
greg |
2.9 |
if ( (c_perez = calloc(5, sizeof(double))) == NULL ) |
284 |
greg |
2.1 |
{ |
285 |
|
|
fprintf(stderr,"Out of memory error in function main !"); |
286 |
greg |
2.9 |
return 1; |
287 |
greg |
2.1 |
} |
288 |
|
|
|
289 |
|
|
printhead(argc, argv); |
290 |
|
|
|
291 |
|
|
computesky(); |
292 |
|
|
printsky(); |
293 |
|
|
|
294 |
greg |
2.9 |
return 0; |
295 |
greg |
2.1 |
} |
296 |
|
|
|
297 |
|
|
|
298 |
greg |
2.9 |
void computesky() /* compute sky parameters */ |
299 |
greg |
2.1 |
{ |
300 |
|
|
|
301 |
|
|
/* new variables */ |
302 |
greg |
2.6 |
int j; |
303 |
greg |
2.1 |
float *lv_mod; /* 145 luminance values*/ |
304 |
|
|
/* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */ |
305 |
greg |
2.9 |
float *theta_o, *phi_o; |
306 |
greg |
2.1 |
double dzeta, gamma; |
307 |
|
|
double normfactor; |
308 |
|
|
|
309 |
|
|
|
310 |
|
|
|
311 |
|
|
/* compute solar direction */ |
312 |
|
|
|
313 |
|
|
if (month) { /* from date and time */ |
314 |
|
|
int jd; |
315 |
|
|
double sd, st; |
316 |
|
|
|
317 |
|
|
jd = jdate(month, day); /* Julian date */ |
318 |
|
|
sd = sdec(jd); /* solar declination */ |
319 |
|
|
if (tsolar) /* solar time */ |
320 |
|
|
st = hour; |
321 |
|
|
else |
322 |
|
|
st = hour + stadj(jd); |
323 |
|
|
altitude = salt(sd, st); |
324 |
|
|
azimuth = sazi(sd, st); |
325 |
|
|
|
326 |
|
|
daynumber = (double)jdate(month, day); |
327 |
|
|
|
328 |
|
|
} |
329 |
greg |
2.9 |
|
330 |
|
|
|
331 |
|
|
|
332 |
|
|
|
333 |
|
|
|
334 |
|
|
/* if loop for the -l option. 01/2013 Sprenger */ |
335 |
|
|
|
336 |
|
|
if (altitude*180/M_PI < sunaltitude_border) { |
337 |
|
|
|
338 |
|
|
if (suppress_warnings==0) |
339 |
|
|
fprintf(stderr, "Warning: sun altitude (%.3f degrees) below the border (%.3f degrees)\n",altitude*180/M_PI,sunaltitude_border); |
340 |
|
|
print_error_sky(); |
341 |
|
|
exit(0); |
342 |
|
|
} |
343 |
|
|
|
344 |
|
|
|
345 |
|
|
|
346 |
|
|
|
347 |
|
|
|
348 |
|
|
if (!cloudy && altitude > 87.*M_PI/180.) { |
349 |
|
|
|
350 |
|
|
if (suppress_warnings==0) { |
351 |
|
|
fprintf(stderr, |
352 |
greg |
2.1 |
"%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n", |
353 |
|
|
progname); |
354 |
greg |
2.9 |
} |
355 |
|
|
altitude = 87.*M_PI/180.; |
356 |
greg |
2.1 |
} |
357 |
greg |
2.9 |
|
358 |
greg |
2.1 |
sundir[0] = -sin(azimuth)*cos(altitude); |
359 |
|
|
sundir[1] = -cos(azimuth)*cos(altitude); |
360 |
|
|
sundir[2] = sin(altitude); |
361 |
|
|
|
362 |
|
|
|
363 |
|
|
/* calculation for the new functions */ |
364 |
greg |
2.9 |
sunzenith = 90 - altitude*180/M_PI; |
365 |
greg |
2.1 |
|
366 |
|
|
|
367 |
|
|
|
368 |
greg |
2.9 |
/* compute the inputs for the calculation of the light distribution over the sky*/ |
369 |
greg |
2.1 |
if (input==0) |
370 |
|
|
{ |
371 |
|
|
check_parametrization(); |
372 |
greg |
2.9 |
diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/ |
373 |
greg |
2.1 |
directirradiance = direct_irradiance_from_sky_clearness(); |
374 |
|
|
check_irradiances(); |
375 |
|
|
|
376 |
|
|
if (output==0 || output==2) |
377 |
|
|
{ |
378 |
greg |
2.9 |
diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
379 |
greg |
2.1 |
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
380 |
|
|
check_illuminances(); |
381 |
|
|
} |
382 |
|
|
} |
383 |
|
|
|
384 |
|
|
|
385 |
|
|
else if (input==1) |
386 |
|
|
{ |
387 |
|
|
check_irradiances(); |
388 |
|
|
skybrightness = sky_brightness(); |
389 |
|
|
skyclearness = sky_clearness(); |
390 |
|
|
check_parametrization(); |
391 |
|
|
|
392 |
|
|
if (output==0 || output==2) |
393 |
|
|
{ |
394 |
greg |
2.9 |
diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
395 |
greg |
2.1 |
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
396 |
|
|
check_illuminances(); |
397 |
|
|
} |
398 |
|
|
|
399 |
|
|
} |
400 |
|
|
|
401 |
|
|
|
402 |
|
|
else if (input==2) |
403 |
|
|
{ |
404 |
|
|
check_illuminances(); |
405 |
|
|
illu_to_irra_index(); |
406 |
|
|
check_parametrization(); |
407 |
|
|
} |
408 |
|
|
|
409 |
|
|
|
410 |
|
|
else if (input==3) |
411 |
|
|
{ |
412 |
|
|
if (altitude<=0) |
413 |
|
|
{ |
414 |
greg |
2.9 |
if (suppress_warnings==0) |
415 |
|
|
fprintf(stderr, "Warning: solar zenith angle larger than 90 degrees; using zero irradiance to proceed\n"); |
416 |
|
|
directirradiance = 0; |
417 |
|
|
diffuseirradiance = 0; |
418 |
|
|
} else { |
419 |
|
|
directirradiance=directirradiance/sin(altitude); |
420 |
greg |
2.1 |
} |
421 |
|
|
check_irradiances(); |
422 |
|
|
skybrightness = sky_brightness(); |
423 |
|
|
skyclearness = sky_clearness(); |
424 |
|
|
check_parametrization(); |
425 |
|
|
|
426 |
|
|
if (output==0 || output==2) |
427 |
|
|
{ |
428 |
greg |
2.9 |
diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
429 |
greg |
2.1 |
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
430 |
|
|
check_illuminances(); |
431 |
|
|
} |
432 |
|
|
|
433 |
|
|
} |
434 |
|
|
|
435 |
|
|
|
436 |
|
|
else {fprintf(stderr,"error in giving the input arguments"); exit(1);} |
437 |
|
|
|
438 |
|
|
|
439 |
|
|
|
440 |
greg |
2.9 |
/* normalization factor for the relative sky luminance distribution, diffuse part*/ |
441 |
greg |
2.1 |
|
442 |
|
|
if ( (lv_mod = malloc(145*sizeof(float))) == NULL) |
443 |
|
|
{ |
444 |
|
|
fprintf(stderr,"Out of memory in function main"); |
445 |
|
|
exit(1); |
446 |
|
|
} |
447 |
|
|
|
448 |
|
|
/* read the angles */ |
449 |
greg |
2.9 |
theta_o = defangle_theta; |
450 |
|
|
phi_o = defangle_phi; |
451 |
greg |
2.1 |
|
452 |
greg |
2.9 |
/* parameters for the perez model */ |
453 |
greg |
2.1 |
coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez); |
454 |
|
|
|
455 |
greg |
2.9 |
/*calculation of the modelled luminance */ |
456 |
greg |
2.1 |
for (j=0;j<145;j++) |
457 |
|
|
{ |
458 |
|
|
theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith)); |
459 |
|
|
*(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez); |
460 |
greg |
2.9 |
// printf("theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); |
461 |
greg |
2.1 |
} |
462 |
|
|
|
463 |
|
|
/* integration of luminance for the normalization factor, diffuse part of the sky*/ |
464 |
|
|
diffnormalization = integ_lv(lv_mod, theta_o); |
465 |
|
|
/*printf("perez integration %lf\n", diffnormalization);*/ |
466 |
|
|
|
467 |
|
|
|
468 |
|
|
|
469 |
|
|
|
470 |
greg |
2.9 |
/*normalization coefficient in lumen or in watt*/ |
471 |
greg |
2.1 |
if (output==0) |
472 |
|
|
{ |
473 |
greg |
2.9 |
diffnormalization = diffuseilluminance/diffnormalization/WHTEFFICACY; |
474 |
greg |
2.1 |
} |
475 |
|
|
else if (output==1) |
476 |
|
|
{ |
477 |
greg |
2.9 |
diffnormalization = diffuseirradiance/diffnormalization; |
478 |
greg |
2.1 |
} |
479 |
|
|
else if (output==2) |
480 |
|
|
{ |
481 |
greg |
2.9 |
diffnormalization = diffuseilluminance/diffnormalization; |
482 |
greg |
2.1 |
} |
483 |
|
|
|
484 |
greg |
2.9 |
else {fprintf(stderr,"Wrong output specification.\n"); exit(1);} |
485 |
greg |
2.1 |
|
486 |
|
|
|
487 |
|
|
|
488 |
|
|
|
489 |
greg |
2.9 |
/* calculation for the solar source */ |
490 |
greg |
2.1 |
if (output==0) |
491 |
greg |
2.9 |
solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY; |
492 |
greg |
2.1 |
|
493 |
|
|
else if (output==1) |
494 |
greg |
2.9 |
solarradiance = directirradiance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180))); |
495 |
greg |
2.1 |
|
496 |
|
|
else |
497 |
greg |
2.9 |
solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180))); |
498 |
greg |
2.1 |
|
499 |
|
|
|
500 |
|
|
|
501 |
|
|
|
502 |
|
|
/* Compute the ground radiance */ |
503 |
|
|
zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez); |
504 |
|
|
zenithbr*=diffnormalization; |
505 |
greg |
2.9 |
|
506 |
greg |
2.1 |
if (skyclearness==1) |
507 |
|
|
normfactor = 0.777778; |
508 |
|
|
|
509 |
|
|
if (skyclearness>=6) |
510 |
|
|
{ |
511 |
greg |
2.9 |
F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]); |
512 |
|
|
normfactor = normsc()/F2/M_PI; |
513 |
greg |
2.1 |
} |
514 |
|
|
|
515 |
|
|
if ( (skyclearness>1) && (skyclearness<6) ) |
516 |
|
|
{ |
517 |
|
|
S_INTER=1; |
518 |
greg |
2.9 |
F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude)); |
519 |
|
|
normfactor = normsc()/F2/M_PI; |
520 |
greg |
2.1 |
} |
521 |
|
|
|
522 |
|
|
groundbr = zenithbr*normfactor; |
523 |
|
|
|
524 |
|
|
if (dosun&&(skyclearness>1)) |
525 |
greg |
2.9 |
groundbr += 6.8e-5/M_PI*solarradiance*sundir[2]; |
526 |
greg |
2.1 |
|
527 |
|
|
groundbr *= gprefl; |
528 |
|
|
|
529 |
|
|
|
530 |
|
|
|
531 |
|
|
return; |
532 |
|
|
} |
533 |
|
|
|
534 |
|
|
|
535 |
|
|
|
536 |
|
|
|
537 |
greg |
2.9 |
void print_error_sky() |
538 |
|
|
{ |
539 |
|
|
sundir[0] = -sin(azimuth)*cos(altitude); |
540 |
|
|
sundir[1] = -cos(azimuth)*cos(altitude); |
541 |
|
|
sundir[2] = sin(altitude); |
542 |
|
|
|
543 |
|
|
printf("\nvoid brightfunc skyfunc\n"); |
544 |
|
|
printf("2 skybright perezlum.cal\n"); |
545 |
|
|
printf("0\n"); |
546 |
|
|
printf("10 0.00 0.00 0.000 0.000 0.000 0.000 0.000 %f %f %f \n", sundir[0], sundir[1], sundir[2]); |
547 |
|
|
} |
548 |
|
|
|
549 |
greg |
2.1 |
|
550 |
|
|
|
551 |
greg |
2.9 |
void printsky() /* print out sky */ |
552 |
greg |
2.1 |
{ |
553 |
|
|
if (dosun&&(skyclearness>1)) |
554 |
greg |
2.9 |
{ |
555 |
greg |
2.1 |
printf("\nvoid light solar\n"); |
556 |
|
|
printf("0\n0\n"); |
557 |
|
|
printf("3 %.3e %.3e %.3e\n", solarradiance, solarradiance, solarradiance); |
558 |
|
|
printf("\nsolar source sun\n"); |
559 |
|
|
printf("0\n0\n"); |
560 |
|
|
printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle); |
561 |
greg |
2.9 |
} else if (dosun) { |
562 |
greg |
2.1 |
printf("\nvoid light solar\n"); |
563 |
|
|
printf("0\n0\n"); |
564 |
|
|
printf("3 0.0 0.0 0.0\n"); |
565 |
|
|
printf("\nsolar source sun\n"); |
566 |
|
|
printf("0\n0\n"); |
567 |
|
|
printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle); |
568 |
greg |
2.9 |
} |
569 |
greg |
2.1 |
|
570 |
|
|
printf("\nvoid brightfunc skyfunc\n"); |
571 |
|
|
printf("2 skybright perezlum.cal\n"); |
572 |
|
|
printf("0\n"); |
573 |
|
|
printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr, |
574 |
greg |
2.9 |
*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4), |
575 |
|
|
sundir[0], sundir[1], sundir[2]); |
576 |
greg |
2.1 |
} |
577 |
|
|
|
578 |
|
|
|
579 |
greg |
2.9 |
void printdefaults() /* print default values */ |
580 |
greg |
2.1 |
{ |
581 |
|
|
printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl); |
582 |
|
|
if (zenithbr > 0.0) |
583 |
|
|
printf("-b %f\t\t\t# Zenith radiance (watts/ster/m^2\n", zenithbr); |
584 |
|
|
else |
585 |
|
|
printf("-t %f\t\t\t# Atmospheric betaturbidity\n", betaturbidity); |
586 |
greg |
2.9 |
printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/M_PI)); |
587 |
|
|
printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/M_PI)); |
588 |
|
|
printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/M_PI)); |
589 |
greg |
2.1 |
} |
590 |
|
|
|
591 |
|
|
|
592 |
greg |
2.9 |
void userror(char* msg) /* print usage error and quit */ |
593 |
greg |
2.1 |
{ |
594 |
|
|
if (msg != NULL) |
595 |
|
|
fprintf(stderr, "%s: Use error - %s\n", progname, msg); |
596 |
greg |
2.9 |
fprintf(stderr, "Usage: %s month day hour [-P|-W|-L|-G] direct_value diffuse_value [options]\n", progname); |
597 |
|
|
fprintf(stderr, "or: %s -ang altitude azimuth [-P|-W|-L|-G] direct_value diffuse_value [options]\n", progname); |
598 |
greg |
2.1 |
fprintf(stderr, " -P epsilon delta (these are the Perez parameters) \n"); |
599 |
|
|
fprintf(stderr, " -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n"); |
600 |
|
|
fprintf(stderr, " -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n"); |
601 |
|
|
fprintf(stderr, " -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n"); |
602 |
|
|
fprintf(stderr, " -O [0|1|2] (0=output in W/m^2/sr visible, 1=output in W/m^2/sr solar, 2=output in candela/m^2), default is 0 \n"); |
603 |
greg |
2.9 |
fprintf(stderr, " gendaylit version 2.00 (2013/01/28) \n"); |
604 |
greg |
2.1 |
exit(1); |
605 |
|
|
} |
606 |
|
|
|
607 |
|
|
|
608 |
|
|
|
609 |
greg |
2.9 |
double normsc() /* compute normalization factor (E0*F2/L0) */ |
610 |
greg |
2.1 |
{ |
611 |
|
|
static double nfc[2][5] = { |
612 |
|
|
/* clear sky approx. */ |
613 |
|
|
{2.766521, 0.547665, -0.369832, 0.009237, 0.059229}, |
614 |
|
|
/* intermediate sky approx. */ |
615 |
|
|
{3.5556, -2.7152, -1.3081, 1.0660, 0.60227}, |
616 |
|
|
}; |
617 |
|
|
register double *nf; |
618 |
|
|
double x, nsc; |
619 |
|
|
register int i; |
620 |
|
|
/* polynomial approximation */ |
621 |
|
|
nf = nfc[S_INTER]; |
622 |
greg |
2.9 |
x = (altitude - M_PI/4.0)/(M_PI/4.0); |
623 |
greg |
2.1 |
nsc = nf[i=4]; |
624 |
|
|
while (i--) |
625 |
|
|
nsc = nsc*x + nf[i]; |
626 |
|
|
|
627 |
|
|
return(nsc); |
628 |
|
|
} |
629 |
|
|
|
630 |
|
|
|
631 |
|
|
|
632 |
greg |
2.9 |
void printhead(int ac, char** av) /* print command header */ |
633 |
greg |
2.1 |
{ |
634 |
|
|
putchar('#'); |
635 |
|
|
while (ac--) { |
636 |
|
|
putchar(' '); |
637 |
|
|
fputs(*av++, stdout); |
638 |
|
|
} |
639 |
|
|
putchar('\n'); |
640 |
|
|
} |
641 |
|
|
|
642 |
|
|
|
643 |
|
|
|
644 |
|
|
|
645 |
|
|
|
646 |
|
|
|
647 |
|
|
/* Perez models */ |
648 |
|
|
|
649 |
|
|
/* Perez global horizontal luminous efficacy model */ |
650 |
|
|
double glob_h_effi_PEREZ() |
651 |
|
|
{ |
652 |
|
|
|
653 |
|
|
double value; |
654 |
|
|
double category_bounds[10], a[10], b[10], c[10], d[10]; |
655 |
|
|
int category_total_number, category_number, i; |
656 |
|
|
|
657 |
|
|
|
658 |
greg |
2.9 |
if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
659 |
|
|
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n"); |
660 |
greg |
2.1 |
|
661 |
|
|
/* initialize category bounds (clearness index bounds) */ |
662 |
|
|
|
663 |
|
|
category_total_number = 8; |
664 |
|
|
|
665 |
|
|
category_bounds[1] = 1; |
666 |
|
|
category_bounds[2] = 1.065; |
667 |
|
|
category_bounds[3] = 1.230; |
668 |
|
|
category_bounds[4] = 1.500; |
669 |
|
|
category_bounds[5] = 1.950; |
670 |
|
|
category_bounds[6] = 2.800; |
671 |
|
|
category_bounds[7] = 4.500; |
672 |
|
|
category_bounds[8] = 6.200; |
673 |
|
|
category_bounds[9] = 12.01; |
674 |
|
|
|
675 |
|
|
|
676 |
|
|
/* initialize model coefficients */ |
677 |
|
|
a[1] = 96.63; |
678 |
|
|
a[2] = 107.54; |
679 |
|
|
a[3] = 98.73; |
680 |
|
|
a[4] = 92.72; |
681 |
|
|
a[5] = 86.73; |
682 |
|
|
a[6] = 88.34; |
683 |
|
|
a[7] = 78.63; |
684 |
|
|
a[8] = 99.65; |
685 |
|
|
|
686 |
|
|
b[1] = -0.47; |
687 |
|
|
b[2] = 0.79; |
688 |
|
|
b[3] = 0.70; |
689 |
|
|
b[4] = 0.56; |
690 |
|
|
b[5] = 0.98; |
691 |
|
|
b[6] = 1.39; |
692 |
|
|
b[7] = 1.47; |
693 |
|
|
b[8] = 1.86; |
694 |
|
|
|
695 |
|
|
c[1] = 11.50; |
696 |
|
|
c[2] = 1.79; |
697 |
|
|
c[3] = 4.40; |
698 |
|
|
c[4] = 8.36; |
699 |
|
|
c[5] = 7.10; |
700 |
|
|
c[6] = 6.06; |
701 |
|
|
c[7] = 4.93; |
702 |
|
|
c[8] = -4.46; |
703 |
|
|
|
704 |
|
|
d[1] = -9.16; |
705 |
|
|
d[2] = -1.19; |
706 |
|
|
d[3] = -6.95; |
707 |
|
|
d[4] = -8.31; |
708 |
|
|
d[5] = -10.94; |
709 |
|
|
d[6] = -7.60; |
710 |
|
|
d[7] = -11.37; |
711 |
|
|
d[8] = -3.15; |
712 |
|
|
|
713 |
|
|
|
714 |
|
|
|
715 |
|
|
|
716 |
|
|
for (i=1; i<=category_total_number; i++) |
717 |
|
|
{ |
718 |
|
|
if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) ) |
719 |
|
|
category_number = i; |
720 |
|
|
} |
721 |
|
|
|
722 |
|
|
value = a[category_number] + b[category_number]*atm_preci_water + |
723 |
greg |
2.9 |
c[category_number]*cos(sunzenith*M_PI/180) + d[category_number]*log(skybrightness); |
724 |
greg |
2.1 |
|
725 |
|
|
return(value); |
726 |
|
|
} |
727 |
|
|
|
728 |
|
|
|
729 |
|
|
/* global horizontal diffuse efficacy model, according to PEREZ */ |
730 |
|
|
double glob_h_diffuse_effi_PEREZ() |
731 |
|
|
{ |
732 |
|
|
double value; |
733 |
|
|
double category_bounds[10], a[10], b[10], c[10], d[10]; |
734 |
|
|
int category_total_number, category_number, i; |
735 |
|
|
|
736 |
|
|
|
737 |
greg |
2.9 |
if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
738 |
|
|
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_effi_PEREZ \n"); |
739 |
greg |
2.1 |
|
740 |
|
|
/* initialize category bounds (clearness index bounds) */ |
741 |
|
|
|
742 |
|
|
category_total_number = 8; |
743 |
|
|
|
744 |
greg |
2.9 |
//XXX: category_bounds > 0.1 |
745 |
greg |
2.1 |
category_bounds[1] = 1; |
746 |
|
|
category_bounds[2] = 1.065; |
747 |
|
|
category_bounds[3] = 1.230; |
748 |
|
|
category_bounds[4] = 1.500; |
749 |
|
|
category_bounds[5] = 1.950; |
750 |
|
|
category_bounds[6] = 2.800; |
751 |
|
|
category_bounds[7] = 4.500; |
752 |
|
|
category_bounds[8] = 6.200; |
753 |
|
|
category_bounds[9] = 12.01; |
754 |
|
|
|
755 |
|
|
|
756 |
|
|
/* initialize model coefficients */ |
757 |
|
|
a[1] = 97.24; |
758 |
|
|
a[2] = 107.22; |
759 |
|
|
a[3] = 104.97; |
760 |
|
|
a[4] = 102.39; |
761 |
|
|
a[5] = 100.71; |
762 |
|
|
a[6] = 106.42; |
763 |
|
|
a[7] = 141.88; |
764 |
|
|
a[8] = 152.23; |
765 |
|
|
|
766 |
|
|
b[1] = -0.46; |
767 |
|
|
b[2] = 1.15; |
768 |
|
|
b[3] = 2.96; |
769 |
|
|
b[4] = 5.59; |
770 |
|
|
b[5] = 5.94; |
771 |
|
|
b[6] = 3.83; |
772 |
|
|
b[7] = 1.90; |
773 |
|
|
b[8] = 0.35; |
774 |
|
|
|
775 |
|
|
c[1] = 12.00; |
776 |
|
|
c[2] = 0.59; |
777 |
|
|
c[3] = -5.53; |
778 |
|
|
c[4] = -13.95; |
779 |
|
|
c[5] = -22.75; |
780 |
|
|
c[6] = -36.15; |
781 |
|
|
c[7] = -53.24; |
782 |
|
|
c[8] = -45.27; |
783 |
|
|
|
784 |
|
|
d[1] = -8.91; |
785 |
|
|
d[2] = -3.95; |
786 |
|
|
d[3] = -8.77; |
787 |
|
|
d[4] = -13.90; |
788 |
|
|
d[5] = -23.74; |
789 |
|
|
d[6] = -28.83; |
790 |
|
|
d[7] = -14.03; |
791 |
|
|
d[8] = -7.98; |
792 |
|
|
|
793 |
|
|
|
794 |
|
|
|
795 |
|
|
|
796 |
greg |
2.9 |
category_number = -1; |
797 |
greg |
2.1 |
for (i=1; i<=category_total_number; i++) |
798 |
|
|
{ |
799 |
|
|
if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) ) |
800 |
|
|
category_number = i; |
801 |
|
|
} |
802 |
|
|
|
803 |
greg |
2.9 |
if (category_number == -1) { |
804 |
|
|
if (suppress_warnings==0) |
805 |
|
|
fprintf(stderr, "ERROR: Model parameters out of range\n"); |
806 |
|
|
print_error_sky(); |
807 |
|
|
exit(1); |
808 |
|
|
} |
809 |
|
|
|
810 |
|
|
|
811 |
|
|
value = a[category_number] + b[category_number]*atm_preci_water + c[category_number]*cos(sunzenith*M_PI/180) + |
812 |
greg |
2.1 |
d[category_number]*log(skybrightness); |
813 |
|
|
|
814 |
|
|
return(value); |
815 |
|
|
} |
816 |
|
|
|
817 |
|
|
|
818 |
|
|
/* direct normal efficacy model, according to PEREZ */ |
819 |
|
|
|
820 |
|
|
double direct_n_effi_PEREZ() |
821 |
|
|
|
822 |
|
|
{ |
823 |
|
|
double value; |
824 |
|
|
double category_bounds[10], a[10], b[10], c[10], d[10]; |
825 |
|
|
int category_total_number, category_number, i; |
826 |
|
|
|
827 |
|
|
|
828 |
greg |
2.9 |
if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
829 |
|
|
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n"); |
830 |
greg |
2.1 |
|
831 |
|
|
|
832 |
|
|
/* initialize category bounds (clearness index bounds) */ |
833 |
|
|
|
834 |
|
|
category_total_number = 8; |
835 |
|
|
|
836 |
|
|
category_bounds[1] = 1; |
837 |
|
|
category_bounds[2] = 1.065; |
838 |
|
|
category_bounds[3] = 1.230; |
839 |
|
|
category_bounds[4] = 1.500; |
840 |
|
|
category_bounds[5] = 1.950; |
841 |
|
|
category_bounds[6] = 2.800; |
842 |
|
|
category_bounds[7] = 4.500; |
843 |
|
|
category_bounds[8] = 6.200; |
844 |
|
|
category_bounds[9] = 12.1; |
845 |
|
|
|
846 |
|
|
|
847 |
|
|
/* initialize model coefficients */ |
848 |
|
|
a[1] = 57.20; |
849 |
|
|
a[2] = 98.99; |
850 |
|
|
a[3] = 109.83; |
851 |
|
|
a[4] = 110.34; |
852 |
|
|
a[5] = 106.36; |
853 |
|
|
a[6] = 107.19; |
854 |
|
|
a[7] = 105.75; |
855 |
|
|
a[8] = 101.18; |
856 |
|
|
|
857 |
|
|
b[1] = -4.55; |
858 |
|
|
b[2] = -3.46; |
859 |
|
|
b[3] = -4.90; |
860 |
|
|
b[4] = -5.84; |
861 |
|
|
b[5] = -3.97; |
862 |
|
|
b[6] = -1.25; |
863 |
|
|
b[7] = 0.77; |
864 |
|
|
b[8] = 1.58; |
865 |
|
|
|
866 |
|
|
c[1] = -2.98; |
867 |
|
|
c[2] = -1.21; |
868 |
|
|
c[3] = -1.71; |
869 |
|
|
c[4] = -1.99; |
870 |
|
|
c[5] = -1.75; |
871 |
|
|
c[6] = -1.51; |
872 |
|
|
c[7] = -1.26; |
873 |
|
|
c[8] = -1.10; |
874 |
|
|
|
875 |
|
|
d[1] = 117.12; |
876 |
|
|
d[2] = 12.38; |
877 |
|
|
d[3] = -8.81; |
878 |
|
|
d[4] = -4.56; |
879 |
|
|
d[5] = -6.16; |
880 |
|
|
d[6] = -26.73; |
881 |
|
|
d[7] = -34.44; |
882 |
|
|
d[8] = -8.29; |
883 |
|
|
|
884 |
|
|
|
885 |
|
|
|
886 |
|
|
for (i=1; i<=category_total_number; i++) |
887 |
|
|
{ |
888 |
|
|
if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) ) |
889 |
|
|
category_number = i; |
890 |
|
|
} |
891 |
|
|
|
892 |
greg |
2.9 |
value = a[category_number] + b[category_number]*atm_preci_water + c[category_number]*exp(5.73*sunzenith*M_PI/180 - 5) + d[category_number]*skybrightness; |
893 |
greg |
2.1 |
|
894 |
|
|
if (value < 0) value = 0; |
895 |
|
|
|
896 |
|
|
return(value); |
897 |
|
|
} |
898 |
|
|
|
899 |
|
|
|
900 |
|
|
/*check the range of epsilon and delta indexes of the perez parametrization*/ |
901 |
|
|
void check_parametrization() |
902 |
|
|
{ |
903 |
greg |
2.9 |
if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) |
904 |
greg |
2.1 |
{ |
905 |
greg |
2.9 |
|
906 |
|
|
/* limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */ |
907 |
|
|
if (skyclearness<skyclearinf){ |
908 |
|
|
skyclearness=skyclearinf; |
909 |
|
|
if (suppress_warnings==0) |
910 |
|
|
fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); |
911 |
|
|
} |
912 |
|
|
if (skyclearness>skyclearsup){ |
913 |
|
|
skyclearness=skyclearsup-0.1; |
914 |
|
|
if (suppress_warnings==0) |
915 |
|
|
fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); |
916 |
|
|
} |
917 |
|
|
if (skybrightness<skybriginf){ |
918 |
|
|
skybrightness=skybriginf; |
919 |
|
|
if (suppress_warnings==0) |
920 |
|
|
fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); |
921 |
|
|
} |
922 |
|
|
if (skybrightness>skybrigsup){ |
923 |
|
|
skybrightness=skybrigsup; |
924 |
|
|
if (suppress_warnings==0) |
925 |
|
|
fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); |
926 |
greg |
2.1 |
} |
927 |
greg |
2.9 |
|
928 |
|
|
return; } |
929 |
greg |
2.1 |
else return; |
930 |
|
|
} |
931 |
|
|
|
932 |
|
|
|
933 |
greg |
2.9 |
/* validity of the direct and diffuse components */ |
934 |
greg |
2.1 |
void check_illuminances() |
935 |
|
|
{ |
936 |
greg |
2.9 |
if (directilluminance < 0) { |
937 |
|
|
fprintf(stderr,"WARNING: direct illuminance < 0. Using 0.0\n"); |
938 |
|
|
directilluminance = 0.0; |
939 |
|
|
} |
940 |
|
|
if (diffuseilluminance < 0) { |
941 |
|
|
fprintf(stderr,"WARNING: diffuse illuminance < 0. Using 0.0\n"); |
942 |
|
|
diffuseilluminance = 0.0; |
943 |
|
|
} |
944 |
|
|
if (directilluminance > solar_constant_l*1000.0) { |
945 |
|
|
fprintf(stderr,"ERROR: direct illuminance exceeds solar constant\n"); |
946 |
|
|
exit(1); |
947 |
greg |
2.1 |
} |
948 |
|
|
} |
949 |
|
|
|
950 |
|
|
|
951 |
|
|
void check_irradiances() |
952 |
|
|
{ |
953 |
greg |
2.9 |
if (directirradiance < 0) { |
954 |
|
|
fprintf(stderr,"WARNING: direct irradiance < 0. Using 0.0\n"); |
955 |
|
|
directirradiance = 0.0; |
956 |
|
|
} |
957 |
|
|
if (diffuseirradiance < 0) { |
958 |
|
|
fprintf(stderr,"WARNING: diffuse irradiance < 0. Using 0.0\n"); |
959 |
|
|
diffuseirradiance = 0.0; |
960 |
|
|
} |
961 |
|
|
if (directirradiance > solar_constant_e) { |
962 |
|
|
fprintf(stderr,"ERROR: direct irradiance exceeds solar constant\n"); |
963 |
|
|
exit(1); |
964 |
|
|
} |
965 |
greg |
2.1 |
} |
966 |
|
|
|
967 |
|
|
|
968 |
|
|
|
969 |
|
|
/* Perez sky's brightness */ |
970 |
|
|
double sky_brightness() |
971 |
|
|
{ |
972 |
|
|
double value; |
973 |
|
|
|
974 |
greg |
2.9 |
value = diffuseirradiance * air_mass() / ( solar_constant_e*get_eccentricity()); |
975 |
greg |
2.1 |
|
976 |
|
|
return(value); |
977 |
|
|
} |
978 |
|
|
|
979 |
|
|
|
980 |
|
|
/* Perez sky's clearness */ |
981 |
|
|
double sky_clearness() |
982 |
|
|
{ |
983 |
greg |
2.9 |
double value; |
984 |
greg |
2.1 |
|
985 |
greg |
2.9 |
value = ( (diffuseirradiance + directirradiance)/(diffuseirradiance) + 1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180 ) / (1 + 1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) ; |
986 |
greg |
2.1 |
|
987 |
greg |
2.9 |
return(value); |
988 |
greg |
2.1 |
} |
989 |
|
|
|
990 |
|
|
|
991 |
|
|
|
992 |
|
|
/* diffus horizontal irradiance from Perez sky's brightness */ |
993 |
greg |
2.9 |
double diffuse_irradiance_from_sky_brightness() |
994 |
greg |
2.1 |
{ |
995 |
|
|
double value; |
996 |
|
|
|
997 |
|
|
value = skybrightness / air_mass() * ( solar_constant_e*get_eccentricity()); |
998 |
|
|
|
999 |
|
|
return(value); |
1000 |
|
|
} |
1001 |
|
|
|
1002 |
|
|
|
1003 |
|
|
/* direct normal irradiance from Perez sky's clearness */ |
1004 |
|
|
double direct_irradiance_from_sky_clearness() |
1005 |
|
|
{ |
1006 |
|
|
double value; |
1007 |
|
|
|
1008 |
greg |
2.9 |
value = diffuse_irradiance_from_sky_brightness(); |
1009 |
|
|
value = value * ( (skyclearness-1) * (1+1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) ); |
1010 |
greg |
2.1 |
|
1011 |
|
|
return(value); |
1012 |
|
|
} |
1013 |
|
|
|
1014 |
|
|
|
1015 |
greg |
2.9 |
void illu_to_irra_index() |
1016 |
greg |
2.1 |
{ |
1017 |
greg |
2.9 |
double test1=0.1, test2=0.1, d_eff; |
1018 |
greg |
2.1 |
int counter=0; |
1019 |
|
|
|
1020 |
greg |
2.9 |
diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l*1000); |
1021 |
greg |
2.1 |
directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000); |
1022 |
|
|
skyclearness = sky_clearness(); |
1023 |
|
|
skybrightness = sky_brightness(); |
1024 |
greg |
2.9 |
check_parametrization(); |
1025 |
greg |
2.1 |
|
1026 |
greg |
2.9 |
while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10) |
1027 |
greg |
2.1 |
|| skyclearness>skyclearinf || skyclearness<skyclearsup |
1028 |
|
|
|| skybrightness>skybriginf || skybrightness<skybrigsup ) |
1029 |
|
|
&& !(counter==5) ) |
1030 |
|
|
{ |
1031 |
|
|
|
1032 |
greg |
2.9 |
test1=diffuseirradiance; |
1033 |
greg |
2.1 |
test2=directirradiance; |
1034 |
|
|
counter++; |
1035 |
|
|
|
1036 |
greg |
2.9 |
diffuseirradiance = diffuseilluminance/glob_h_diffuse_effi_PEREZ(); |
1037 |
|
|
d_eff = direct_n_effi_PEREZ(); |
1038 |
|
|
if (d_eff < 0.1) |
1039 |
|
|
directirradiance = 0; |
1040 |
|
|
else |
1041 |
|
|
directirradiance = directilluminance/d_eff; |
1042 |
greg |
2.1 |
|
1043 |
|
|
skybrightness = sky_brightness(); |
1044 |
|
|
skyclearness = sky_clearness(); |
1045 |
greg |
2.9 |
check_parametrization(); |
1046 |
|
|
|
1047 |
greg |
2.1 |
} |
1048 |
|
|
|
1049 |
|
|
|
1050 |
|
|
return; |
1051 |
|
|
} |
1052 |
|
|
|
1053 |
greg |
2.9 |
static int get_numlin(float epsilon) |
1054 |
greg |
2.1 |
{ |
1055 |
greg |
2.9 |
if (epsilon < 1.065) |
1056 |
|
|
return 0; |
1057 |
|
|
else if (epsilon < 1.230) |
1058 |
|
|
return 1; |
1059 |
|
|
else if (epsilon < 1.500) |
1060 |
|
|
return 2; |
1061 |
|
|
else if (epsilon < 1.950) |
1062 |
|
|
return 3; |
1063 |
|
|
else if (epsilon < 2.800) |
1064 |
|
|
return 4; |
1065 |
|
|
else if (epsilon < 4.500) |
1066 |
|
|
return 5; |
1067 |
|
|
else if (epsilon < 6.200) |
1068 |
|
|
return 6; |
1069 |
|
|
return 7; |
1070 |
greg |
2.1 |
} |
1071 |
|
|
|
1072 |
|
|
/* sky luminance perez model */ |
1073 |
greg |
2.9 |
double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]) |
1074 |
greg |
2.1 |
{ |
1075 |
|
|
float x[5][4]; |
1076 |
|
|
int i,j,num_lin; |
1077 |
|
|
double c_perez[5]; |
1078 |
|
|
|
1079 |
|
|
if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) ) |
1080 |
|
|
{ |
1081 |
greg |
2.9 |
fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez!\n"); |
1082 |
greg |
2.1 |
exit(1); |
1083 |
|
|
} |
1084 |
|
|
|
1085 |
|
|
/* correction de modele de Perez solar energy ...*/ |
1086 |
|
|
if ( (epsilon > 1.065) && (epsilon < 2.8) ) |
1087 |
|
|
{ |
1088 |
|
|
if ( Delta < 0.2 ) Delta = 0.2; |
1089 |
|
|
} |
1090 |
|
|
|
1091 |
greg |
2.9 |
num_lin = get_numlin(epsilon); |
1092 |
greg |
2.1 |
|
1093 |
|
|
for (i=0;i<5;i++) |
1094 |
|
|
for (j=0;j<4;j++) |
1095 |
|
|
{ |
1096 |
|
|
x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j); |
1097 |
|
|
/* printf("x %d %d vaut %f\n",i,j,x[i][j]); */ |
1098 |
|
|
} |
1099 |
|
|
|
1100 |
|
|
|
1101 |
|
|
if (num_lin) |
1102 |
|
|
{ |
1103 |
|
|
for (i=0;i<5;i++) |
1104 |
|
|
c_perez[i] = x[i][0] + x[i][1]*Z + Delta * (x[i][2] + x[i][3]*Z); |
1105 |
|
|
} |
1106 |
|
|
else |
1107 |
|
|
{ |
1108 |
|
|
c_perez[0] = x[0][0] + x[0][1]*Z + Delta * (x[0][2] + x[0][3]*Z); |
1109 |
|
|
c_perez[1] = x[1][0] + x[1][1]*Z + Delta * (x[1][2] + x[1][3]*Z); |
1110 |
|
|
c_perez[4] = x[4][0] + x[4][1]*Z + Delta * (x[4][2] + x[4][3]*Z); |
1111 |
|
|
c_perez[2] = exp( pow(Delta*(x[2][0]+x[2][1]*Z),x[2][2])) - x[2][3]; |
1112 |
|
|
c_perez[3] = -exp( Delta*(x[3][0]+x[3][1]*Z) )+x[3][2]+Delta*x[3][3]; |
1113 |
|
|
} |
1114 |
|
|
|
1115 |
|
|
|
1116 |
|
|
return (1 + c_perez[0]*exp(c_perez[1]/cos(dzeta)) ) * |
1117 |
|
|
(1 + c_perez[2]*exp(c_perez[3]*gamma) + |
1118 |
|
|
c_perez[4]*cos(gamma)*cos(gamma) ); |
1119 |
|
|
} |
1120 |
|
|
|
1121 |
|
|
|
1122 |
|
|
|
1123 |
|
|
/* coefficients for the sky luminance perez model */ |
1124 |
greg |
2.9 |
void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]) |
1125 |
greg |
2.1 |
{ |
1126 |
|
|
float x[5][4]; |
1127 |
|
|
int i,j,num_lin; |
1128 |
|
|
|
1129 |
|
|
if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) ) |
1130 |
|
|
{ |
1131 |
|
|
fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n"); |
1132 |
|
|
exit(1); |
1133 |
|
|
} |
1134 |
|
|
|
1135 |
|
|
/* correction du modele de Perez solar energy ...*/ |
1136 |
|
|
if ( (epsilon > 1.065) && (epsilon < 2.8) ) |
1137 |
|
|
{ |
1138 |
|
|
if ( Delta < 0.2 ) Delta = 0.2; |
1139 |
|
|
} |
1140 |
|
|
|
1141 |
greg |
2.9 |
num_lin = get_numlin(epsilon); |
1142 |
|
|
|
1143 |
|
|
//fprintf(stderr,"numlin %d\n", num_lin); |
1144 |
greg |
2.1 |
|
1145 |
|
|
for (i=0;i<5;i++) |
1146 |
|
|
for (j=0;j<4;j++) |
1147 |
|
|
{ |
1148 |
|
|
x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j); |
1149 |
|
|
/* printf("x %d %d vaut %f\n",i,j,x[i][j]); */ |
1150 |
|
|
} |
1151 |
|
|
|
1152 |
|
|
|
1153 |
|
|
if (num_lin) |
1154 |
|
|
{ |
1155 |
|
|
for (i=0;i<5;i++) |
1156 |
|
|
*(c_perez+i) = x[i][0] + x[i][1]*Z + Delta * (x[i][2] + x[i][3]*Z); |
1157 |
|
|
|
1158 |
|
|
} |
1159 |
|
|
else |
1160 |
|
|
{ |
1161 |
|
|
*(c_perez+0) = x[0][0] + x[0][1]*Z + Delta * (x[0][2] + x[0][3]*Z); |
1162 |
|
|
*(c_perez+1) = x[1][0] + x[1][1]*Z + Delta * (x[1][2] + x[1][3]*Z); |
1163 |
|
|
*(c_perez+4) = x[4][0] + x[4][1]*Z + Delta * (x[4][2] + x[4][3]*Z); |
1164 |
|
|
*(c_perez+2) = exp( pow(Delta*(x[2][0]+x[2][1]*Z),x[2][2])) - x[2][3]; |
1165 |
|
|
*(c_perez+3) = -exp( Delta*(x[3][0]+x[3][1]*Z) )+x[3][2]+Delta*x[3][3]; |
1166 |
|
|
|
1167 |
|
|
|
1168 |
|
|
} |
1169 |
|
|
|
1170 |
|
|
|
1171 |
|
|
return; |
1172 |
|
|
} |
1173 |
|
|
|
1174 |
|
|
|
1175 |
|
|
/* degrees into radians */ |
1176 |
|
|
double radians(double degres) |
1177 |
|
|
{ |
1178 |
greg |
2.9 |
return degres*M_PI/180.0; |
1179 |
greg |
2.1 |
} |
1180 |
|
|
|
1181 |
|
|
/* radian into degrees */ |
1182 |
|
|
double degres(double radians) |
1183 |
|
|
{ |
1184 |
greg |
2.9 |
return radians/M_PI*180.0; |
1185 |
greg |
2.1 |
} |
1186 |
|
|
|
1187 |
|
|
/* calculation of the angles dzeta and gamma */ |
1188 |
|
|
void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z) |
1189 |
|
|
{ |
1190 |
|
|
*dzeta = theta; /* dzeta = phi */ |
1191 |
|
|
if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1 && (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi) < 1.1 ) ) |
1192 |
|
|
*gamma = 0; |
1193 |
|
|
else if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1.1 ) |
1194 |
|
|
{ |
1195 |
|
|
printf("error in calculation of gamma (angle between point and sun"); |
1196 |
|
|
exit(3); |
1197 |
|
|
} |
1198 |
|
|
else |
1199 |
|
|
*gamma = acos(cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)); |
1200 |
|
|
} |
1201 |
|
|
|
1202 |
|
|
|
1203 |
|
|
|
1204 |
|
|
/********************************************************************************/ |
1205 |
|
|
/* Fonction: integ_lv */ |
1206 |
|
|
/* */ |
1207 |
|
|
/* In: float *lv,*theta */ |
1208 |
|
|
/* int sun_pos */ |
1209 |
|
|
/* */ |
1210 |
|
|
/* Out: double */ |
1211 |
|
|
/* */ |
1212 |
|
|
/* Update: 29/08/93 */ |
1213 |
|
|
/* */ |
1214 |
|
|
/* Rem: */ |
1215 |
|
|
/* */ |
1216 |
|
|
/* But: calcul l'integrale de luminance relative sans la dir. du soleil */ |
1217 |
|
|
/* */ |
1218 |
|
|
/********************************************************************************/ |
1219 |
|
|
double integ_lv(float *lv,float *theta) |
1220 |
|
|
{ |
1221 |
|
|
int i; |
1222 |
|
|
double buffer=0.0; |
1223 |
|
|
|
1224 |
|
|
for (i=0;i<145;i++) |
1225 |
|
|
buffer += (*(lv+i))*cos(radians(*(theta+i))); |
1226 |
|
|
|
1227 |
greg |
2.9 |
return buffer*2*M_PI/144; |
1228 |
greg |
2.1 |
|
1229 |
|
|
} |
1230 |
|
|
|
1231 |
|
|
|
1232 |
|
|
|
1233 |
|
|
|
1234 |
|
|
|
1235 |
|
|
|
1236 |
|
|
/* enter day number(double), return E0 = square(R0/R): eccentricity correction factor */ |
1237 |
|
|
|
1238 |
|
|
double get_eccentricity() |
1239 |
|
|
{ |
1240 |
|
|
double day_angle; |
1241 |
|
|
double E0; |
1242 |
|
|
|
1243 |
greg |
2.9 |
day_angle = 2*M_PI*(daynumber - 1)/365; |
1244 |
greg |
2.1 |
E0 = 1.00011+0.034221*cos(day_angle)+0.00128*sin(day_angle)+ |
1245 |
|
|
0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle); |
1246 |
|
|
|
1247 |
|
|
return (E0); |
1248 |
|
|
|
1249 |
|
|
} |
1250 |
|
|
|
1251 |
|
|
|
1252 |
|
|
/* enter sunzenith angle (degrees) return relative air mass (double) */ |
1253 |
|
|
double air_mass() |
1254 |
|
|
{ |
1255 |
|
|
double m; |
1256 |
|
|
|
1257 |
|
|
if (sunzenith>90) |
1258 |
|
|
{ |
1259 |
greg |
2.9 |
fprintf(stderr, "Solar zenith angle larger than 90 degrees in function air_mass()\n"); |
1260 |
greg |
2.1 |
exit(1); |
1261 |
|
|
} |
1262 |
|
|
|
1263 |
greg |
2.9 |
m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) ); |
1264 |
greg |
2.1 |
return(m); |
1265 |
|
|
} |
1266 |
|
|
|
1267 |
|
|
|
1268 |
|
|
|