1 |
– |
#ifndef lint |
2 |
– |
static const char RCSid[] = "$Id$"; |
3 |
– |
#endif |
1 |
|
/* Copyright (c) 1994,2006 *Fraunhofer Institut for Solar Energy Systems |
2 |
|
* Heidenhofstr. 2, D-79110 Freiburg, Germany |
3 |
|
* *Agence de l'Environnement et de la Maitrise de l'Energie |
4 |
|
* Centre de Valbonne, 500 route des Lucioles, 06565 Sophia Antipolis Cedex, France |
5 |
|
* *BOUYGUES |
6 |
|
* 1 Avenue Eugene Freyssinet, Saint-Quentin-Yvelines, France |
10 |
– |
* |
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 |
7 |
|
*/ |
8 |
|
|
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 |
– |
|
9 |
|
#include <stdio.h> |
10 |
|
#include <string.h> |
11 |
|
#include <math.h> |
12 |
|
#include <stdlib.h> |
13 |
|
|
14 |
|
#include "color.h" |
15 |
+ |
#include "sun.h" |
16 |
|
#include "paths.h" |
17 |
|
|
18 |
|
#define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]) |
22 |
|
/*static char *rcsid="$Header$";*/ |
23 |
|
|
24 |
|
float coeff_perez[] = { |
25 |
< |
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, |
26 |
< |
-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, |
27 |
< |
-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, |
28 |
< |
-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, |
29 |
< |
-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, |
30 |
< |
-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, |
31 |
< |
-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, |
32 |
< |
-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}; |
25 |
> |
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, |
26 |
> |
0.9738,0.2809,0.0356,-0.1246,-0.5718,0.9938,-1.2219,-0.7730,1.4148,1.1016,-0.2054,0.0367,-3.9128,0.9156, |
27 |
> |
6.9750,0.1774,6.4477,-0.1239,-1.5798,-0.5081,-1.7812,0.1080,0.2624,0.0672,-0.2190,-0.4285,-1.1000,-0.2515, |
28 |
> |
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, |
29 |
> |
-2.6204,-0.0156,0.1597,0.4199,-0.5562,-0.5484,-0.6654,-0.2672,0.7117,0.7234,-0.6219,-5.6812,2.6297, |
30 |
> |
33.3389,-18.3000,-62.2500,52.0781,-3.5000,0.0016,1.1477,0.1062,0.4659,-0.3296,-0.0876,-0.0329,-0.6000, |
31 |
> |
-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, |
32 |
> |
1.4062,0.3988,0.0032,0.0766,-0.0656,-0.1294,-1.0156,-0.3670,1.0078,1.4051,0.2875,-0.5328,-3.8500,3.3750, |
33 |
> |
14.0000,-0.9999,-7.1406,7.5469,-3.4000,-0.1078,-1.0750,1.5702,-0.0672,0.4016,0.3017,-0.4844,-1.0000, |
34 |
> |
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, |
35 |
> |
0.2656,1.0468,-0.3788,-2.4517,1.4656,-1.0500,0.0289,0.4260,0.3590,-0.3250,0.1156,0.7781,0.0025,31.0625, |
36 |
> |
-14.5000,-46.1148,55.3750,-7.2312,0.4050,13.3500,0.6234,1.5000,-0.6426,1.8564,0.5636}; |
37 |
|
|
38 |
|
|
39 |
< |
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}; |
39 |
> |
float defangle_theta[] = { |
40 |
> |
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, |
41 |
> |
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, |
42 |
> |
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, |
43 |
> |
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, |
44 |
> |
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, |
45 |
> |
24, 24, 24, 24, 24, 24, 24, 24, 12, 12, 12, 12, 12, 12, 0}; |
46 |
|
|
47 |
< |
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}; |
47 |
> |
float defangle_phi[] = { |
48 |
> |
0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264, |
49 |
> |
276, 288, 300, 312, 324, 336, 348, 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, |
50 |
> |
192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 0, 15, 30, 45, 60, 75, 90, 105, |
51 |
> |
120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 15, 30, 45, 60, 75, |
52 |
> |
90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 20, 40, 60, |
53 |
> |
80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 0, 30, 60, 90, 120, 150, 180, 210, |
54 |
> |
240, 270, 300, 330, 0, 60, 120, 180, 240, 300, 0}; |
55 |
|
|
56 |
|
|
57 |
|
|
63 |
|
double diffuse_irradiance_from_sky_brightness(); |
64 |
|
double direct_irradiance_from_sky_clearness(); |
65 |
|
|
66 |
+ |
/* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : */ |
67 |
+ |
/* input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */ |
68 |
|
|
64 |
– |
/* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */ |
69 |
|
double glob_h_effi_PEREZ(); |
70 |
|
double glob_h_diffuse_effi_PEREZ(); |
71 |
|
double direct_n_effi_PEREZ(); |
72 |
+ |
|
73 |
|
/*likelihood check of the epsilon, delta, direct and diffuse components*/ |
74 |
|
void check_parametrization(); |
75 |
|
void check_irradiances(); |
77 |
|
void illu_to_irra_index(); |
78 |
|
void print_error_sky(); |
79 |
|
|
80 |
< |
|
76 |
< |
/* Perez sky luminance model */ |
77 |
< |
double calc_rel_lum_perez(double dzeta,double gamma,double Z, |
78 |
< |
double epsilon,double Delta,float coeff_perez[]); |
79 |
< |
/* coefficients for the sky luminance perez model */ |
80 |
> |
double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]); |
81 |
|
void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]); |
82 |
|
double radians(double degres); |
83 |
|
double degres(double radians); |
85 |
|
double integ_lv(float *lv,float *theta); |
86 |
|
|
87 |
|
void printdefaults(); |
88 |
+ |
void check_sun_position(); |
89 |
|
void computesky(); |
90 |
|
void printhead(int ac, char** av); |
91 |
|
void userror(char* msg); |
97 |
|
double get_eccentricity(); |
98 |
|
double air_mass(); |
99 |
|
|
100 |
< |
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); |
102 |
< |
extern double sazi(double sd, double st); |
100 |
> |
double solar_sunset(); |
101 |
> |
double solar_sunrise(); |
102 |
> |
double stadj(); |
103 |
> |
int jdate(int month, int day); |
104 |
|
|
105 |
|
|
106 |
+ |
|
107 |
|
/* sun calculation constants */ |
108 |
< |
extern double s_latitude; |
109 |
< |
extern double s_longitude; |
110 |
< |
extern double s_meridian; |
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 */ |
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 */ |
119 |
> |
const double skyclearinf = 1.0; /* 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; |
134 |
|
/* definition of the sky conditions through the Perez parametrization */ |
135 |
|
double skyclearness = 0; |
136 |
|
double skybrightness = 0; |
137 |
< |
double solarradiance; /*radiance of the sun disk and of the circumsolar area*/ |
138 |
< |
double diffuseilluminance, directilluminance, diffuseirradiance, directirradiance; |
139 |
< |
double sunzenith, daynumber=150, atm_preci_water=2; |
137 |
> |
double solarradiance; |
138 |
> |
double diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance; |
139 |
> |
double sunzenith, daynumber, atm_preci_water=2; |
140 |
|
|
141 |
< |
double sunaltitude_border = 0; |
141 |
> |
/*double sunaltitude_border = 0;*/ |
142 |
|
double diffnormalization = 0; |
143 |
< |
double dirnormalization = 0; |
143 |
> |
double dirnormalization = 0; |
144 |
|
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*/ |
146 |
> |
int output=0; /* define the unit of the output (sky luminance or radiance): */ |
147 |
> |
/* visible watt=0, solar watt=1, lumen=2 */ |
148 |
> |
int input=0; /* define the input for the calulation */ |
149 |
|
|
150 |
|
int suppress_warnings=0; |
151 |
|
|
152 |
|
/* default values */ |
153 |
< |
int cloudy = 0; /* 1=standard, 2=uniform */ |
154 |
< |
int dosun = 1; |
153 |
> |
int cloudy = 0; /* 1=standard, 2=uniform */ |
154 |
> |
int dosun = 1; |
155 |
|
double zenithbr = -1.0; |
156 |
|
double betaturbidity = 0.1; |
157 |
|
double gprefl = 0.2; |
163 |
|
double F2; |
164 |
|
double solarbr = 0.0; |
165 |
|
int u_solar = 0; /* -1=irradiance, 1=radiance */ |
166 |
+ |
float timeinterval = 0; |
167 |
|
|
168 |
< |
char *progname; |
169 |
< |
char errmsg[128]; |
168 |
> |
char *progname; |
169 |
> |
char errmsg[128]; |
170 |
|
|
171 |
|
|
172 |
+ |
|
173 |
+ |
|
174 |
|
int main(int argc, char** argv) |
175 |
|
{ |
176 |
|
int i; |
205 |
|
cloudy = 0; |
206 |
|
dosun = argv[i][0] == '+'; |
207 |
|
break; |
202 |
– |
case 'r': |
208 |
|
case 'R': |
209 |
|
u_solar = argv[i][1] == 'R' ? -1 : 1; |
210 |
|
solarbr = atof(argv[++i]); |
264 |
|
diffuseirradiance = atof(argv[++i]); |
265 |
|
break; |
266 |
|
|
267 |
< |
case 'l': |
267 |
> |
case 'E': /* Erbs model based on the */ |
268 |
> |
input = 4; /* global-horizontal irradiance [W/m^2] */ |
269 |
> |
globalirradiance = atof(argv[++i]); |
270 |
> |
break; |
271 |
> |
|
272 |
> |
/* |
273 |
> |
case 'l': |
274 |
|
sunaltitude_border = atof(argv[++i]); |
275 |
|
break; |
276 |
+ |
*/ |
277 |
+ |
|
278 |
+ |
case 'i': |
279 |
+ |
timeinterval = atof(argv[++i]); |
280 |
+ |
break; |
281 |
|
|
282 |
|
|
283 |
|
default: |
293 |
|
progname, (s_longitude-s_meridian)*12/M_PI); |
294 |
|
|
295 |
|
|
296 |
< |
/* allocation dynamique de memoire pour les pointeurs */ |
296 |
> |
/* dynamic memory allocation for the pointers */ |
297 |
> |
|
298 |
|
if ( (c_perez = calloc(5, sizeof(double))) == NULL ) |
299 |
|
{ |
300 |
< |
fprintf(stderr,"Out of memory error in function main !"); |
300 |
> |
fprintf(stderr,"Out of memory error in function main"); |
301 |
|
return 1; |
302 |
|
} |
303 |
|
|
304 |
|
printhead(argc, argv); |
288 |
– |
|
305 |
|
computesky(); |
306 |
+ |
|
307 |
+ |
if(*(c_perez+1)>0) |
308 |
+ |
{ |
309 |
+ |
fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1)); |
310 |
+ |
print_error_sky(); |
311 |
+ |
exit(1); |
312 |
+ |
} |
313 |
+ |
|
314 |
|
printsky(); |
291 |
– |
|
315 |
|
return 0; |
316 |
|
} |
317 |
|
|
318 |
|
|
319 |
< |
void computesky() /* compute sky parameters */ |
319 |
> |
|
320 |
> |
|
321 |
> |
|
322 |
> |
|
323 |
> |
|
324 |
> |
|
325 |
> |
void computesky() |
326 |
|
{ |
327 |
|
|
299 |
– |
/* new variables */ |
328 |
|
int j; |
329 |
< |
float *lv_mod; /* 145 luminance values*/ |
330 |
< |
/* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */ |
329 |
> |
|
330 |
> |
float *lv_mod; /* 145 luminance values */ |
331 |
|
float *theta_o, *phi_o; |
332 |
|
double dzeta, gamma; |
333 |
|
double normfactor; |
334 |
+ |
double erbs_s0, erbs_kt; |
335 |
|
|
336 |
|
|
308 |
– |
|
337 |
|
/* compute solar direction */ |
338 |
< |
|
338 |
> |
|
339 |
|
if (month) { /* from date and time */ |
340 |
|
int jd; |
341 |
|
double sd, st; |
346 |
|
st = hour; |
347 |
|
else |
348 |
|
st = hour + stadj(jd); |
349 |
+ |
|
350 |
+ |
|
351 |
+ |
if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) { |
352 |
+ |
print_error_sky(); |
353 |
+ |
exit(1); |
354 |
+ |
} |
355 |
+ |
|
356 |
+ |
|
357 |
+ |
if(timeinterval) { |
358 |
+ |
|
359 |
+ |
if(timeinterval<0) { |
360 |
+ |
fprintf(stderr, "time interval negative\n"); |
361 |
+ |
exit(1); |
362 |
+ |
} |
363 |
+ |
|
364 |
+ |
if(fabs(solar_sunrise(month,day)-st)<timeinterval/60) { |
365 |
+ |
|
366 |
+ |
fprintf(stderr, "Solar position corrected at %d %d %.3f\n",month,day,hour); |
367 |
+ |
st= (st+timeinterval/120+solar_sunrise(month,day))/2; |
368 |
+ |
} |
369 |
+ |
|
370 |
+ |
if(fabs(solar_sunset(month,day)-st)<timeinterval/60) { |
371 |
+ |
fprintf(stderr, "Solar position corrected at %d %d %.3f\n",month,day,hour); |
372 |
+ |
st= (st-timeinterval/120+solar_sunset(month,day))/2; |
373 |
+ |
} |
374 |
+ |
} |
375 |
+ |
|
376 |
+ |
|
377 |
|
altitude = salt(sd, st); |
378 |
|
azimuth = sazi(sd, st); |
379 |
|
|
380 |
|
daynumber = (double)jdate(month, day); |
381 |
< |
|
381 |
> |
|
382 |
|
} |
383 |
|
|
384 |
|
|
385 |
< |
|
386 |
< |
|
387 |
< |
|
388 |
< |
/* if loop for the -l option. 01/2013 Sprenger */ |
333 |
< |
|
385 |
> |
|
386 |
> |
|
387 |
> |
/* if loop for the -l option. W.Sprenger (01/2013) */ |
388 |
> |
/* |
389 |
|
if (altitude*180/M_PI < sunaltitude_border) { |
390 |
|
|
391 |
< |
if (suppress_warnings==0) |
392 |
< |
fprintf(stderr, "Warning: sun altitude (%.3f degrees) below the border (%.3f degrees)\n",altitude*180/M_PI,sunaltitude_border); |
391 |
> |
if (suppress_warnings==0) { |
392 |
> |
fprintf(stderr, "Warning: sun altitude (%.3f degrees) below the border (%.3f degrees)\n",altitude*180/M_PI,sunaltitude_border); |
393 |
> |
} |
394 |
|
print_error_sky(); |
395 |
< |
exit(0); |
395 |
> |
exit(1); |
396 |
|
} |
397 |
< |
|
398 |
< |
|
397 |
> |
*/ |
398 |
> |
|
399 |
|
|
344 |
– |
|
345 |
– |
|
400 |
|
if (!cloudy && altitude > 87.*M_PI/180.) { |
401 |
|
|
402 |
|
if (suppress_warnings==0) { |
418 |
|
|
419 |
|
|
420 |
|
/* compute the inputs for the calculation of the light distribution over the sky*/ |
421 |
< |
if (input==0) |
421 |
> |
if (input==0) /* P */ |
422 |
|
{ |
423 |
|
check_parametrization(); |
424 |
|
diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/ |
434 |
|
} |
435 |
|
|
436 |
|
|
437 |
< |
else if (input==1) |
437 |
> |
else if (input==1) /* W */ |
438 |
|
{ |
439 |
|
check_irradiances(); |
440 |
|
skybrightness = sky_brightness(); |
451 |
|
} |
452 |
|
|
453 |
|
|
454 |
< |
else if (input==2) |
454 |
> |
else if (input==2) /* L */ |
455 |
|
{ |
456 |
|
check_illuminances(); |
457 |
|
illu_to_irra_index(); |
459 |
|
} |
460 |
|
|
461 |
|
|
462 |
< |
else if (input==3) |
462 |
> |
else if (input==3) /* G */ |
463 |
|
{ |
464 |
|
if (altitude<=0) |
465 |
|
{ |
468 |
|
directirradiance = 0; |
469 |
|
diffuseirradiance = 0; |
470 |
|
} else { |
471 |
< |
directirradiance=directirradiance/sin(altitude); |
471 |
> |
|
472 |
> |
directirradiance=directirradiance/sin(altitude); |
473 |
|
} |
474 |
+ |
|
475 |
|
check_irradiances(); |
476 |
|
skybrightness = sky_brightness(); |
477 |
|
skyclearness = sky_clearness(); |
486 |
|
|
487 |
|
} |
488 |
|
|
489 |
+ |
|
490 |
+ |
else if (input==4) /* E */ /* Implementation of the Erbs model. W.Sprenger (04/13) */ |
491 |
+ |
{ |
492 |
+ |
|
493 |
+ |
if (altitude<=0) |
494 |
+ |
{ |
495 |
+ |
if (suppress_warnings==0 && globalirradiance > 50) |
496 |
+ |
fprintf(stderr, "Warning: global irradiance higher than 50 W/m^2 while the sun altitude is lower than zero\n"); |
497 |
+ |
globalirradiance = 0; diffuseirradiance = 0; directirradiance = 0; |
498 |
+ |
|
499 |
+ |
} else { |
500 |
+ |
|
501 |
+ |
erbs_s0 = solar_constant_e*get_eccentricity()*sin(altitude); |
502 |
+ |
|
503 |
+ |
if (globalirradiance>erbs_s0) |
504 |
+ |
{ |
505 |
+ |
if (suppress_warnings==0) |
506 |
+ |
fprintf(stderr, "Warning: global irradiance is higher than the time-dependent solar constant s0\n"); |
507 |
+ |
|
508 |
+ |
globalirradiance=erbs_s0*0.999; |
509 |
+ |
} |
510 |
+ |
|
511 |
+ |
erbs_kt=globalirradiance/erbs_s0; |
512 |
+ |
|
513 |
+ |
if (erbs_kt<=0.22) diffuseirradiance=globalirradiance*(1-0.09*erbs_kt); |
514 |
+ |
else if (erbs_kt<=0.8) diffuseirradiance=globalirradiance*(0.9511-0.1604*erbs_kt+4.388*pow(erbs_kt,2)-16.638*pow(erbs_kt,3)+12.336*pow(erbs_kt,4)); |
515 |
+ |
else if (erbs_kt<1) diffuseirradiance=globalirradiance*(0.165); |
516 |
+ |
|
517 |
+ |
directirradiance=globalirradiance-diffuseirradiance; |
518 |
+ |
|
519 |
+ |
printf("# erbs_s0, erbs_kt, irr_dir_h, irr_diff: %.3f %.3f %.3f %.3f\n", erbs_s0, erbs_kt, directirradiance, diffuseirradiance); |
520 |
+ |
printf("# WARNING: the -E option is only recommended for a rough estimation!"); |
521 |
+ |
|
522 |
+ |
directirradiance=directirradiance/sin(altitude); |
523 |
+ |
|
524 |
+ |
} |
525 |
+ |
|
526 |
+ |
check_irradiances(); |
527 |
+ |
skybrightness = sky_brightness(); |
528 |
+ |
skyclearness = sky_clearness(); |
529 |
+ |
check_parametrization(); |
530 |
+ |
|
531 |
+ |
if (output==0 || output==2) |
532 |
+ |
{ |
533 |
+ |
diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
534 |
+ |
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
535 |
+ |
check_illuminances(); |
536 |
+ |
} |
537 |
+ |
|
538 |
+ |
} |
539 |
+ |
|
540 |
+ |
|
541 |
+ |
|
542 |
|
|
543 |
< |
else {fprintf(stderr,"error in giving the input arguments"); exit(1);} |
543 |
> |
else {fprintf(stderr,"error at the input arguments"); exit(1);} |
544 |
|
|
545 |
|
|
546 |
|
|
547 |
|
/* normalization factor for the relative sky luminance distribution, diffuse part*/ |
548 |
< |
|
548 |
> |
|
549 |
|
if ( (lv_mod = malloc(145*sizeof(float))) == NULL) |
550 |
|
{ |
551 |
|
fprintf(stderr,"Out of memory in function main"); |
555 |
|
/* read the angles */ |
556 |
|
theta_o = defangle_theta; |
557 |
|
phi_o = defangle_phi; |
558 |
+ |
|
559 |
|
|
560 |
|
/* parameters for the perez model */ |
561 |
|
coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez); |
562 |
|
|
563 |
+ |
|
564 |
+ |
|
565 |
+ |
|
566 |
+ |
|
567 |
|
/*calculation of the modelled luminance */ |
568 |
|
for (j=0;j<145;j++) |
569 |
|
{ |
570 |
|
theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith)); |
571 |
+ |
|
572 |
|
*(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez); |
573 |
< |
// printf("theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); |
573 |
> |
|
574 |
> |
/* fprintf(stderr,"theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); */ |
575 |
|
} |
576 |
< |
|
576 |
> |
|
577 |
|
/* integration of luminance for the normalization factor, diffuse part of the sky*/ |
578 |
+ |
|
579 |
|
diffnormalization = integ_lv(lv_mod, theta_o); |
463 |
– |
/*printf("perez integration %lf\n", diffnormalization);*/ |
580 |
|
|
465 |
– |
|
581 |
|
|
582 |
|
|
583 |
|
/*normalization coefficient in lumen or in watt*/ |
661 |
|
|
662 |
|
|
663 |
|
|
664 |
+ |
|
665 |
+ |
|
666 |
+ |
|
667 |
+ |
|
668 |
+ |
double solar_sunset(int month,int day) |
669 |
+ |
{ |
670 |
+ |
float W; |
671 |
+ |
extern double s_latitude; |
672 |
+ |
W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day)))); |
673 |
+ |
return(12+(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15)); |
674 |
+ |
} |
675 |
+ |
|
676 |
+ |
|
677 |
+ |
double solar_sunrise(int month,int day) |
678 |
+ |
{ |
679 |
+ |
float W; |
680 |
+ |
extern double s_latitude; |
681 |
+ |
W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day)))); |
682 |
+ |
return(12-(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15)); |
683 |
+ |
} |
684 |
+ |
|
685 |
+ |
|
686 |
+ |
|
687 |
+ |
|
688 |
+ |
|
689 |
+ |
|
690 |
+ |
|
691 |
+ |
|
692 |
+ |
|
693 |
+ |
|
694 |
+ |
|
695 |
|
void printsky() /* print out sky */ |
696 |
|
{ |
697 |
|
if (dosun&&(skyclearness>1)) |
711 |
|
printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle); |
712 |
|
} |
713 |
|
|
714 |
+ |
|
715 |
|
printf("\nvoid brightfunc skyfunc\n"); |
716 |
|
printf("2 skybright perezlum.cal\n"); |
717 |
|
printf("0\n"); |
718 |
|
printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr, |
719 |
|
*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4), |
720 |
|
sundir[0], sundir[1], sundir[2]); |
721 |
+ |
|
722 |
|
} |
723 |
|
|
724 |
|
|
738 |
|
void userror(char* msg) /* print usage error and quit */ |
739 |
|
{ |
740 |
|
if (msg != NULL) |
741 |
< |
fprintf(stderr, "%s: Use error - %s\n", progname, msg); |
742 |
< |
fprintf(stderr, "Usage: %s month day hour [-P|-W|-L|-G] direct_value diffuse_value [options]\n", progname); |
743 |
< |
fprintf(stderr, "or: %s -ang altitude azimuth [-P|-W|-L|-G] direct_value diffuse_value [options]\n", progname); |
741 |
> |
fprintf(stderr, "%s: Use error - %s\n\n", progname, msg); |
742 |
> |
fprintf(stderr, "Usage: %s month day hour [...]\n", progname); |
743 |
> |
fprintf(stderr, " or: %s -ang altitude azimuth [...]\n", progname); |
744 |
> |
fprintf(stderr, " followed by: -P epsilon delta [options]\n"); |
745 |
> |
fprintf(stderr, " or: [-W|-L|-G] direct_value diffuse_value [options]\n"); |
746 |
> |
fprintf(stderr, " or: -E global_irradiance [options]\n\n"); |
747 |
> |
fprintf(stderr, " Description:\n"); |
748 |
|
fprintf(stderr, " -P epsilon delta (these are the Perez parameters) \n"); |
749 |
|
fprintf(stderr, " -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n"); |
750 |
|
fprintf(stderr, " -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n"); |
751 |
|
fprintf(stderr, " -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n"); |
752 |
+ |
fprintf(stderr, " -E global-horizontal-irradiance (W/m^2)\n\n"); |
753 |
+ |
fprintf(stderr, " Output specification with option:\n"); |
754 |
|
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"); |
755 |
< |
fprintf(stderr, " gendaylit version 2.00 (2013/01/28) \n"); |
755 |
> |
fprintf(stderr, " gendaylit version 2.3 (2013/08/08) \n\n"); |
756 |
|
exit(1); |
757 |
|
} |
758 |
|
|
794 |
|
|
795 |
|
|
796 |
|
|
643 |
– |
|
644 |
– |
|
797 |
|
/* Perez models */ |
798 |
|
|
799 |
|
/* Perez global horizontal luminous efficacy model */ |
803 |
|
double value; |
804 |
|
double category_bounds[10], a[10], b[10], c[10], d[10]; |
805 |
|
int category_total_number, category_number, i; |
806 |
< |
|
807 |
< |
|
808 |
< |
if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
809 |
< |
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n"); |
810 |
< |
|
806 |
> |
|
807 |
> |
check_parametrization(); |
808 |
> |
|
809 |
> |
|
810 |
> |
/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
811 |
> |
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n"); */ |
812 |
> |
|
813 |
> |
|
814 |
|
/* initialize category bounds (clearness index bounds) */ |
815 |
|
|
816 |
|
category_total_number = 8; |
886 |
|
double category_bounds[10], a[10], b[10], c[10], d[10]; |
887 |
|
int category_total_number, category_number, i; |
888 |
|
|
889 |
+ |
|
890 |
+ |
|
891 |
|
|
892 |
< |
if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
893 |
< |
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_effi_PEREZ \n"); |
892 |
> |
check_parametrization(); |
893 |
> |
|
894 |
> |
|
895 |
> |
/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
896 |
> |
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */ |
897 |
> |
|
898 |
> |
|
899 |
|
|
900 |
|
/* initialize category bounds (clearness index bounds) */ |
901 |
|
|
902 |
|
category_total_number = 8; |
903 |
|
|
904 |
< |
//XXX: category_bounds > 0.1 |
904 |
> |
//XXX: category_bounds > 0.1 |
905 |
|
category_bounds[1] = 1; |
906 |
|
category_bounds[2] = 1.065; |
907 |
|
category_bounds[3] = 1.230; |
962 |
|
|
963 |
|
if (category_number == -1) { |
964 |
|
if (suppress_warnings==0) |
965 |
< |
fprintf(stderr, "ERROR: Model parameters out of range\n"); |
965 |
> |
fprintf(stderr, "ERROR: Model parameters out of range, skyclearness = %lf \n", skyclearness); |
966 |
|
print_error_sky(); |
967 |
|
exit(1); |
968 |
|
} |
972 |
|
d[category_number]*log(skybrightness); |
973 |
|
|
974 |
|
return(value); |
975 |
+ |
|
976 |
|
} |
977 |
|
|
978 |
|
|
979 |
+ |
|
980 |
|
/* direct normal efficacy model, according to PEREZ */ |
981 |
|
|
982 |
|
double direct_n_effi_PEREZ() |
1067 |
|
|
1068 |
|
/* limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */ |
1069 |
|
if (skyclearness<skyclearinf){ |
906 |
– |
skyclearness=skyclearinf; |
1070 |
|
if (suppress_warnings==0) |
1071 |
< |
fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); |
1071 |
> |
/* fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */ |
1072 |
> |
skyclearness=skyclearinf; |
1073 |
|
} |
1074 |
|
if (skyclearness>skyclearsup){ |
911 |
– |
skyclearness=skyclearsup-0.1; |
1075 |
|
if (suppress_warnings==0) |
1076 |
< |
fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); |
1076 |
> |
/* fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */ |
1077 |
> |
skyclearness=skyclearsup-0.1; |
1078 |
|
} |
1079 |
|
if (skybrightness<skybriginf){ |
916 |
– |
skybrightness=skybriginf; |
1080 |
|
if (suppress_warnings==0) |
1081 |
< |
fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); |
1081 |
> |
/* fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */ |
1082 |
> |
skybrightness=skybriginf; |
1083 |
|
} |
1084 |
|
if (skybrightness>skybrigsup){ |
921 |
– |
skybrightness=skybrigsup; |
1085 |
|
if (suppress_warnings==0) |
1086 |
< |
fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); |
1086 |
> |
/* fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */ |
1087 |
> |
skybrightness=skybrigsup; |
1088 |
|
} |
1089 |
|
|
1090 |
|
return; } |
1174 |
|
} |
1175 |
|
|
1176 |
|
|
1177 |
+ |
|
1178 |
+ |
|
1179 |
|
void illu_to_irra_index() |
1180 |
|
{ |
1181 |
|
double test1=0.1, test2=0.1, d_eff; |
1186 |
|
skyclearness = sky_clearness(); |
1187 |
|
skybrightness = sky_brightness(); |
1188 |
|
check_parametrization(); |
1189 |
< |
|
1189 |
> |
|
1190 |
> |
|
1191 |
|
while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10) |
1192 |
< |
|| skyclearness>skyclearinf || skyclearness<skyclearsup |
1193 |
< |
|| skybrightness>skybriginf || skybrightness<skybrigsup ) |
1194 |
< |
&& !(counter==5) ) |
1192 |
> |
|| (!(skyclearness<skyclearinf || skyclearness>skyclearsup)) |
1193 |
> |
|| (!(skybrightness<skybriginf || skybrightness>skybrigsup)) ) |
1194 |
> |
&& !(counter==9) ) |
1195 |
|
{ |
1196 |
< |
|
1196 |
> |
|
1197 |
|
test1=diffuseirradiance; |
1198 |
|
test2=directirradiance; |
1199 |
|
counter++; |
1200 |
|
|
1201 |
|
diffuseirradiance = diffuseilluminance/glob_h_diffuse_effi_PEREZ(); |
1202 |
|
d_eff = direct_n_effi_PEREZ(); |
1203 |
+ |
|
1204 |
+ |
|
1205 |
|
if (d_eff < 0.1) |
1206 |
|
directirradiance = 0; |
1207 |
< |
else |
1207 |
> |
else |
1208 |
|
directirradiance = directilluminance/d_eff; |
1209 |
|
|
1210 |
|
skybrightness = sky_brightness(); |
1211 |
|
skyclearness = sky_clearness(); |
1212 |
|
check_parametrization(); |
1213 |
|
|
1214 |
+ |
/*fprintf(stderr,"skyclearness = %lf, skybrightness = %lf, directirradiance = %lf, diffuseirradiance = %lf\n",skyclearness, skybrightness, directirradiance, diffuseirradiance);*/ |
1215 |
+ |
|
1216 |
|
} |
1217 |
|
|
1218 |
|
|
1241 |
|
/* sky luminance perez model */ |
1242 |
|
double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]) |
1243 |
|
{ |
1244 |
+ |
|
1245 |
|
float x[5][4]; |
1246 |
|
int i,j,num_lin; |
1247 |
|
double c_perez[5]; |
1257 |
|
{ |
1258 |
|
if ( Delta < 0.2 ) Delta = 0.2; |
1259 |
|
} |
1260 |
< |
|
1260 |
> |
|
1261 |
> |
|
1262 |
|
num_lin = get_numlin(epsilon); |
1263 |
< |
|
1263 |
> |
|
1264 |
|
for (i=0;i<5;i++) |
1265 |
|
for (j=0;j<4;j++) |
1266 |
|
{ |
1267 |
|
x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j); |
1268 |
< |
/* printf("x %d %d vaut %f\n",i,j,x[i][j]); */ |
1268 |
> |
/* fprintf(stderr,"x %d %d vaut %f\n",i,j,x[i][j]); */ |
1269 |
|
} |
1270 |
|
|
1271 |
|
|
1308 |
|
{ |
1309 |
|
if ( Delta < 0.2 ) Delta = 0.2; |
1310 |
|
} |
1311 |
< |
|
1311 |
> |
|
1312 |
> |
|
1313 |
|
num_lin = get_numlin(epsilon); |
1314 |
|
|
1315 |
< |
//fprintf(stderr,"numlin %d\n", num_lin); |
1315 |
> |
/*fprintf(stderr,"numlin %d\n", num_lin);*/ |
1316 |
|
|
1317 |
|
for (i=0;i<5;i++) |
1318 |
|
for (j=0;j<4;j++) |
1344 |
|
} |
1345 |
|
|
1346 |
|
|
1347 |
+ |
|
1348 |
|
/* degrees into radians */ |
1349 |
|
double radians(double degres) |
1350 |
|
{ |
1351 |
|
return degres*M_PI/180.0; |
1352 |
|
} |
1353 |
|
|
1354 |
+ |
|
1355 |
|
/* radian into degrees */ |
1356 |
|
double degres(double radians) |
1357 |
|
{ |
1358 |
|
return radians/M_PI*180.0; |
1359 |
|
} |
1360 |
|
|
1361 |
+ |
|
1362 |
|
/* calculation of the angles dzeta and gamma */ |
1363 |
|
void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z) |
1364 |
|
{ |
1376 |
|
|
1377 |
|
|
1378 |
|
|
1202 |
– |
/********************************************************************************/ |
1203 |
– |
/* Fonction: integ_lv */ |
1204 |
– |
/* */ |
1205 |
– |
/* In: float *lv,*theta */ |
1206 |
– |
/* int sun_pos */ |
1207 |
– |
/* */ |
1208 |
– |
/* Out: double */ |
1209 |
– |
/* */ |
1210 |
– |
/* Update: 29/08/93 */ |
1211 |
– |
/* */ |
1212 |
– |
/* Rem: */ |
1213 |
– |
/* */ |
1214 |
– |
/* But: calcul l'integrale de luminance relative sans la dir. du soleil */ |
1215 |
– |
/* */ |
1216 |
– |
/********************************************************************************/ |
1379 |
|
double integ_lv(float *lv,float *theta) |
1380 |
|
{ |
1381 |
|
int i; |
1382 |
|
double buffer=0.0; |
1383 |
< |
|
1383 |
> |
|
1384 |
|
for (i=0;i<145;i++) |
1385 |
+ |
{ |
1386 |
|
buffer += (*(lv+i))*cos(radians(*(theta+i))); |
1387 |
< |
|
1387 |
> |
} |
1388 |
> |
|
1389 |
|
return buffer*2*M_PI/144; |
1226 |
– |
|
1390 |
|
} |
1391 |
|
|
1392 |
|
|
1393 |
|
|
1231 |
– |
|
1232 |
– |
|
1233 |
– |
|
1394 |
|
/* enter day number(double), return E0 = square(R0/R): eccentricity correction factor */ |
1395 |
|
|
1396 |
|
double get_eccentricity() |
1403 |
|
0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle); |
1404 |
|
|
1405 |
|
return (E0); |
1246 |
– |
|
1406 |
|
} |
1407 |
|
|
1408 |
|
|
1410 |
|
double air_mass() |
1411 |
|
{ |
1412 |
|
double m; |
1254 |
– |
|
1413 |
|
if (sunzenith>90) |
1414 |
|
{ |
1415 |
|
fprintf(stderr, "Solar zenith angle larger than 90 degrees in function air_mass()\n"); |
1416 |
|
exit(1); |
1417 |
|
} |
1260 |
– |
|
1418 |
|
m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) ); |
1419 |
|
return(m); |
1420 |
|
} |