ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gendaylit.c
Revision: 2.21
Committed: Thu Jan 28 19:03:15 2021 UTC (3 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4, HEAD
Changes since 2.20: +10 -2 lines
Log Message:
feat: Jan added -d option for setting dewpoint as requested by Zack Rogers

File Contents

# Content
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
7 * print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018
8 * version 2.6 (2021/01/29): dew point dependency added according to Perez publication 1990 (W -> atm_preci_water=exp(0.07*Td-0.075) ). by J. Wienold, EPFL
9 */
10
11 #define _USE_MATH_DEFINES
12 #include <stdio.h>
13 #include <string.h>
14 #include <math.h>
15 #include <stdlib.h>
16
17 #include "color.h"
18 #include "sun.h"
19 #include "paths.h"
20
21 #define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
22
23 double normsc();
24
25 /*static char *rcsid="$Header: /usr/local/cvs/radiance/ray/src/gen/gendaylit.c,v 2.20 2020/11/18 17:18:41 greg Exp $";*/
26
27 float coeff_perez[] = {
28 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,
29 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,
30 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,
31 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,
32 -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,
33 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,
34 -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,
35 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,
36 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,
37 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,
38 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,
39 -14.5000,-46.1148,55.3750,-7.2312,0.4050,13.3500,0.6234,1.5000,-0.6426,1.8564,0.5636};
40
41
42 float defangle_theta[] = {
43 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,
44 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,
45 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,
46 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,
47 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,
48 24, 24, 24, 24, 24, 24, 24, 24, 12, 12, 12, 12, 12, 12, 0};
49
50 float defangle_phi[] = {
51 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264,
52 276, 288, 300, 312, 324, 336, 348, 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180,
53 192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 0, 15, 30, 45, 60, 75, 90, 105,
54 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 15, 30, 45, 60, 75,
55 90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 20, 40, 60,
56 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 0, 30, 60, 90, 120, 150, 180, 210,
57 240, 270, 300, 330, 0, 60, 120, 180, 240, 300, 0};
58 /* default values for Berlin */
59 float locus[] = {
60 -4.843e9,2.5568e6,0.24282e3,0.23258,-4.843e9,2.5568e6,0.24282e3,0.23258,-1.2848,1.7519,-0.093786};
61
62
63
64 /* Perez sky parametrization: epsilon and delta calculations from the direct and diffuse irradiances */
65 double sky_brightness();
66 double sky_clearness();
67
68 /* calculation of the direct and diffuse components from the Perez parametrization */
69 double diffuse_irradiance_from_sky_brightness();
70 double direct_irradiance_from_sky_clearness();
71
72 /* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : */
73 /* input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */
74
75 double glob_h_effi_PEREZ();
76 double glob_h_diffuse_effi_PEREZ();
77 double direct_n_effi_PEREZ();
78
79 /*likelihood check of the epsilon, delta, direct and diffuse components*/
80 void check_parametrization();
81 void check_irradiances();
82 void check_illuminances();
83 void illu_to_irra_index();
84 void print_error_sky();
85
86 double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]);
87 void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]);
88 double radians(double degres);
89 double degres(double radians);
90 void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z);
91 double integ_lv(float *lv,float *theta);
92
93 void printdefaults();
94 void check_sun_position();
95 void computesky();
96 void printhead(int ac, char** av);
97 void usage_error(char* msg);
98 void printsky();
99
100 FILE * frlibopen(char* fname);
101
102 /* astronomy and geometry*/
103 double get_eccentricity();
104 double air_mass();
105
106 double solar_sunset(int month, int day);
107 double solar_sunrise(int month, int day);
108
109 const double AU = 149597890E3;
110 const double solar_constant_e = 1367; /* solar constant W/m^2 */
111 const double solar_constant_l = 127500; /* solar constant lux */
112
113 const double half_sun_angle = 0.2665;
114 const double half_direct_angle = 2.85;
115
116 const double skyclearinf = 1.0; /* limitations for the variation of the Perez parameters */
117 const double skyclearsup = 12.01;
118 const double skybriginf = 0.01;
119 const double skybrigsup = 0.6;
120
121
122
123 /* required values */
124 int year = 0; /* year (optional) */
125 int month, day; /* date */
126 double hour; /* time */
127 int tsolar; /* 0=standard, 1=solar */
128 double altitude, azimuth; /* or solar angles */
129
130
131
132 /* definition of the sky conditions through the Perez parametrization */
133 double skyclearness = 0;
134 double skybrightness = 0;
135 double solarradiance;
136 double diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance;
137 double sunzenith, daynumber, atm_preci_water, Td=10.97353115;
138
139 /*double sunaltitude_border = 0;*/
140 double diffnormalization = 0;
141 double dirnormalization = 0;
142 double *c_perez;
143
144 int output=0; /* define the unit of the output (sky luminance or radiance): */
145 /* visible watt=0, solar watt=1, lumen=2 */
146 int input=0; /* define the input for the calulation */
147 int color_output=0;
148 int suppress_warnings=0;
149
150 /* default values */
151 int cloudy = 0; /* 1=standard, 2=uniform */
152 int dosun = 1;
153 double zenithbr = -1.0;
154 double betaturbidity = 0.1;
155 double gprefl = 0.2;
156 int S_INTER=0;
157
158
159 /* computed values */
160 double sundir[3];
161 double groundbr = 0;
162 double F2;
163 double solarbr = 0.0;
164 int u_solar = 0; /* -1=irradiance, 1=radiance */
165 float timeinterval = 0;
166
167 char *progname;
168 char errmsg[128];
169
170 double st;
171
172
173 int main(int argc, char** argv)
174 {
175 int i;
176
177 progname = argv[0];
178 if (argc == 2 && !strcmp(argv[1], "-defaults")) {
179 printdefaults();
180 return 0;
181 }
182 if (argc < 4)
183 usage_error("arg count");
184 if (!strcmp(argv[1], "-ang")) {
185 altitude = atof(argv[2]) * (M_PI/180);
186 azimuth = atof(argv[3]) * (M_PI/180);
187 month = 0;
188 } else {
189 month = atoi(argv[1]);
190 if (month < 1 || month > 12)
191 usage_error("bad month");
192 day = atoi(argv[2]);
193 if (day < 1 || day > 31)
194 usage_error("bad day");
195 hour = atof(argv[3]);
196 if (hour < 0 || hour >= 24)
197 usage_error("bad hour");
198 tsolar = argv[3][0] == '+';
199 }
200 for (i = 4; i < argc; i++)
201 if (argv[i][0] == '-' || argv[i][0] == '+')
202 switch (argv[i][1]) {
203 case 'd':
204 Td = atof(argv[++i]);
205 if (Td < -40 || Td > 40) {
206 Td=10.97353115; }
207 break;
208 case 's':
209 cloudy = 0;
210 dosun = argv[i][0] == '+';
211 break;
212 case 'y':
213 year = atoi(argv[++i]);
214 break;
215 case 'R':
216 u_solar = argv[i][1] == 'R' ? -1 : 1;
217 solarbr = atof(argv[++i]);
218 break;
219 case 'c':
220 cloudy = argv[i][0] == '+' ? 2 : 1;
221 dosun = 0;
222 break;
223 case 'C':
224 if (argv[i][2] == 'I' && argv[i][3] == 'E' ) {
225 locus[0] = -4.607e9;
226 locus[1] = 2.9678e6;
227 locus[2] = 0.09911e3;
228 locus[3] = 0.244063;
229 locus[4] = -2.0064e9;
230 locus[5] = 1.9018e6;
231 locus[6] = 0.24748e3;
232 locus[7] = 0.23704;
233 locus[8] = -3.0;
234 locus[9] = 2.87;
235 locus[10] = -0.275;
236 }else{ color_output = 1;
237 }
238 break;
239 case 'l':
240 locus[0] = atof(argv[++i]);
241 locus[1] = atof(argv[++i]);
242 locus[2] = atof(argv[++i]);
243 locus[3] = atof(argv[++i]);
244 locus[4] = locus[0];
245 locus[5] = locus[1];
246 locus[6] = locus[2];
247 locus[7] = locus[3];
248 locus[8] = atof(argv[++i]);
249 locus[9] = atof(argv[++i]);
250 locus[10] = atof(argv[++i]);
251 break;
252
253 case 't':
254 betaturbidity = atof(argv[++i]);
255 break;
256 case 'w':
257 suppress_warnings = 1;
258 break;
259 case 'b':
260 zenithbr = atof(argv[++i]);
261 break;
262 case 'g':
263 gprefl = atof(argv[++i]);
264 break;
265 case 'a':
266 s_latitude = atof(argv[++i]) * (M_PI/180);
267 break;
268 case 'o':
269 s_longitude = atof(argv[++i]) * (M_PI/180);
270 break;
271 case 'm':
272 s_meridian = atof(argv[++i]) * (M_PI/180);
273 break;
274
275 case 'O':
276 output = atof(argv[++i]); /*define the unit of the output of the program:
277 sky and sun luminance/radiance
278 (0==W visible, 1==W solar radiation, 2==lm) */
279 break;
280
281 case 'P':
282 input = 0; /* Perez parameters: epsilon, delta */
283 skyclearness = atof(argv[++i]);
284 skybrightness = atof(argv[++i]);
285 break;
286
287 case 'W': /* direct normal Irradiance [W/m^2] */
288 input = 1; /* diffuse horizontal Irrad. [W/m^2] */
289 directirradiance = atof(argv[++i]);
290 diffuseirradiance = atof(argv[++i]);
291 break;
292
293 case 'L': /* direct normal Illuminance [Lux] */
294 input = 2; /* diffuse horizontal Ill. [Lux] */
295 directilluminance = atof(argv[++i]);
296 diffuseilluminance = atof(argv[++i]);
297 break;
298
299 case 'G': /* direct horizontal Irradiance [W/m^2] */
300 input = 3; /* diffuse horizontal Irrad. [W/m^2] */
301 directirradiance = atof(argv[++i]);
302 diffuseirradiance = atof(argv[++i]);
303 break;
304
305 case 'E': /* Erbs model based on the */
306 input = 4; /* global-horizontal irradiance [W/m^2] */
307 globalirradiance = atof(argv[++i]);
308 break;
309
310 case 'i':
311 timeinterval = atof(argv[++i]);
312 break;
313
314
315 default:
316 sprintf(errmsg, "unknown option: %s", argv[i]);
317 usage_error(errmsg);
318 }
319 else
320 usage_error("bad option");
321
322 if (month && !tsolar && fabs(s_meridian-s_longitude) > 45*M_PI/180)
323 fprintf(stderr,"%s: warning: %.1f hours btwn. standard meridian and longitude\n",
324 progname, (s_longitude-s_meridian)*12/M_PI);
325
326
327 /* dynamic memory allocation for the pointers */
328 if ( (c_perez = calloc(5, sizeof(double))) == NULL )
329 { fprintf(stderr,"Out of memory error in function main"); return 1; }
330
331
332 atm_preci_water=exp(0.07*Td-0.075);
333 printhead(argc, argv);
334 computesky();
335 printsky();
336 return 0;
337
338 }
339
340
341
342
343
344 void computesky()
345 {
346
347 int j;
348
349 float *lv_mod; /* 145 luminance values */
350 float *theta_o, *phi_o;
351 double dzeta, gamma;
352 double normfactor;
353 double erbs_s0, erbs_kt;
354
355
356 /* compute solar direction */
357
358 if (month) { /* from date and time */
359 double sd;
360
361 st = hour;
362 if (year) { /* Michalsky algorithm? */
363 double mjd = mjdate(year, month, day, hour);
364 if (tsolar)
365 sd = msdec(mjd, NULL);
366 else
367 sd = msdec(mjd, &st);
368 } else {
369 int jd = jdate(month, day); /* Julian date */
370 sd = sdec(jd); /* solar declination */
371 if (!tsolar) /* get solar time? */
372 st = hour + stadj(jd);
373 }
374
375 if(timeinterval) {
376
377 if(timeinterval<0) {
378 fprintf(stderr, "time interval negative\n");
379 exit(1);
380 }
381
382 if(fabs(solar_sunrise(month,day)-st)<=timeinterval/120) {
383 st= (st+timeinterval/120+solar_sunrise(month,day))/2;
384 if(suppress_warnings==0)
385 { fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); }
386 }
387
388 if(fabs(solar_sunset(month,day)-st)<timeinterval/120) {
389 st= (st-timeinterval/120+solar_sunset(month,day))/2;
390 if(suppress_warnings==0)
391 { fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); }
392 }
393
394 if((st<solar_sunrise(month,day)-timeinterval/120) || (st>solar_sunset(month,day)+timeinterval/120)) {
395 if(suppress_warnings==0)
396 { fprintf(stderr, "Warning: sun position too low, printing error sky at %d %d %.3f\n",month,day,hour); }
397 altitude = salt(sd, st);
398 azimuth = sazi(sd, st);
399 print_error_sky();
400 exit(0);
401 }
402 }
403 else
404
405 if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) {
406 if(suppress_warnings==0)
407 { fprintf(stderr, "Warning: sun altitude below zero at time step %i %i %.2f, printing error sky\n",month,day,hour); }
408 altitude = salt(sd, st);
409 azimuth = sazi(sd, st);
410 print_error_sky();
411 exit(0);
412 }
413
414 altitude = salt(sd, st);
415 azimuth = sazi(sd, st);
416
417 daynumber = (double)jdate(month, day);
418
419 }
420
421
422
423
424
425 if (!cloudy && altitude > 87.*M_PI/180.) {
426
427 if (suppress_warnings==0) {
428 fprintf(stderr,
429 "%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n",
430 progname);
431 }
432 altitude = 87.*M_PI/180.;
433 }
434
435
436
437 sundir[0] = -sin(azimuth)*cos(altitude);
438 sundir[1] = -cos(azimuth)*cos(altitude);
439 sundir[2] = sin(altitude);
440
441
442 /* calculation for the new functions */
443 sunzenith = 90 - altitude*180/M_PI;
444
445
446 /* compute the inputs for the calculation of the light distribution over the sky*/
447 if (input==0) /* P */
448 {
449 check_parametrization();
450 diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
451 directirradiance = direct_irradiance_from_sky_clearness();
452 check_irradiances();
453
454 if (output==0 || output==2)
455 {
456 diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
457 directilluminance = directirradiance*direct_n_effi_PEREZ();
458 check_illuminances();
459 }
460 }
461
462
463 else if (input==1) /* W */
464 {
465 check_irradiances();
466 skybrightness = sky_brightness();
467 skyclearness = sky_clearness();
468
469 check_parametrization();
470
471 if (output==0 || output==2)
472 {
473 diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
474 directilluminance = directirradiance*direct_n_effi_PEREZ();
475 check_illuminances();
476 }
477
478 }
479
480
481 else if (input==2) /* L */
482 {
483 check_illuminances();
484 illu_to_irra_index();
485 check_parametrization();
486 }
487
488
489 else if (input==3) /* G */
490 {
491 if (altitude<=0)
492 {
493 if (suppress_warnings==0)
494 fprintf(stderr, "Warning: sun altitude < 0, proceed with irradiance values of zero\n");
495 directirradiance = 0;
496 diffuseirradiance = 0;
497 } else {
498
499 directirradiance=directirradiance/sin(altitude);
500 }
501
502 check_irradiances();
503 skybrightness = sky_brightness();
504 skyclearness = sky_clearness();
505 check_parametrization();
506
507 if (output==0 || output==2)
508 {
509 diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
510 directilluminance = directirradiance*direct_n_effi_PEREZ();
511 check_illuminances();
512 }
513
514 }
515
516
517 else if (input==4) /* E */ /* Implementation of the Erbs model. W.Sprenger (04/13) */
518 {
519
520 if (altitude<=0)
521 {
522 if (suppress_warnings==0 && globalirradiance > 50)
523 fprintf(stderr, "Warning: global irradiance higher than 50 W/m^2 while the sun altitude is lower than zero\n");
524 globalirradiance = 0; diffuseirradiance = 0; directirradiance = 0;
525
526 } else {
527
528 erbs_s0 = solar_constant_e*get_eccentricity()*sin(altitude);
529
530 if (globalirradiance>erbs_s0)
531 {
532 if (suppress_warnings==0)
533 fprintf(stderr, "Warning: global irradiance is higher than the time-dependent solar constant s0\n");
534 globalirradiance=erbs_s0*0.999;
535 }
536
537 erbs_kt=globalirradiance/erbs_s0;
538
539 if (erbs_kt<=0.22) diffuseirradiance=globalirradiance*(1-0.09*erbs_kt);
540 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));
541 else if (erbs_kt<1) diffuseirradiance=globalirradiance*(0.165);
542
543 directirradiance=globalirradiance-diffuseirradiance;
544
545 printf("# erbs_s0, erbs_kt, irr_dir_h, irr_diff: %.3f %.3f %.3f %.3f\n", erbs_s0, erbs_kt, directirradiance, diffuseirradiance);
546 printf("# WARNING: the -E option is only recommended for a rough estimation!\n");
547
548 directirradiance=directirradiance/sin(altitude);
549
550 }
551
552 check_irradiances();
553 skybrightness = sky_brightness();
554 skyclearness = sky_clearness();
555 check_parametrization();
556
557 if (output==0 || output==2)
558 {
559 diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
560 directilluminance = directirradiance*direct_n_effi_PEREZ();
561 check_illuminances();
562 }
563
564 }
565
566
567
568
569 else { fprintf(stderr,"error at the input arguments"); exit(1); }
570
571
572
573 /* normalization factor for the relative sky luminance distribution, diffuse part*/
574
575 if ( (lv_mod = malloc(145*sizeof(float))) == NULL)
576 {
577 fprintf(stderr,"Out of memory in function main");
578 exit(1);
579 }
580
581 /* read the angles */
582 theta_o = defangle_theta;
583 phi_o = defangle_phi;
584
585
586 /* parameters for the perez model */
587 coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez);
588
589
590
591 /*calculation of the modelled luminance */
592 for (j=0;j<145;j++)
593 {
594 theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith));
595
596 *(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez);
597
598 /* fprintf(stderr,"theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); */
599 }
600
601 /* integration of luminance for the normalization factor, diffuse part of the sky*/
602
603 diffnormalization = integ_lv(lv_mod, theta_o);
604
605
606
607 /*normalization coefficient in lumen or in watt*/
608 if (output==0)
609 {
610 diffnormalization = diffuseilluminance/diffnormalization/WHTEFFICACY;
611 }
612 else if (output==1)
613 {
614 diffnormalization = diffuseirradiance/diffnormalization;
615 }
616 else if (output==2)
617 {
618 diffnormalization = diffuseilluminance/diffnormalization;
619 }
620
621 else {fprintf(stderr,"Wrong output specification.\n"); exit(1);}
622
623
624
625
626 /* calculation for the solar source */
627 if (output==0)
628 solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY;
629
630 else if (output==1)
631 solarradiance = directirradiance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
632
633 else
634 solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
635
636
637
638 /* Compute the ground radiance */
639 zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
640 zenithbr*=diffnormalization;
641
642 if (skyclearness==1)
643 normfactor = 0.777778;
644
645 if (skyclearness>=6)
646 {
647 F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]);
648 normfactor = normsc()/F2/M_PI;
649 }
650
651 if ( (skyclearness>1) && (skyclearness<6) )
652 {
653 S_INTER=1;
654 F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude));
655 normfactor = normsc()/F2/M_PI;
656 }
657
658 groundbr = zenithbr*normfactor;
659
660 if (dosun&&(skyclearness>1))
661 groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];
662
663 groundbr *= gprefl;
664
665
666
667 if(*(c_perez+1)>0)
668 {
669 if(suppress_warnings==0)
670 { fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1));}
671 print_error_sky();
672 exit(0);
673 }
674
675
676 return;
677 }
678
679
680
681
682
683 double solar_sunset(int month,int day)
684 {
685 float W;
686 extern double s_latitude;
687 W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day))));
688 return(12+(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15));
689 }
690
691
692
693
694 double solar_sunrise(int month,int day)
695 {
696 float W;
697 extern double s_latitude;
698 W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day))));
699 return(12-(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15));
700 }
701
702
703
704
705 void printsky()
706 {
707
708 printf("# Local solar time: %.2f\n", st);
709 printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI);
710 printf("# epsilon, delta, atmospheric precipitable water content : %.4f %.4f %.4f \n", skyclearness, skybrightness,atm_preci_water );
711
712
713 if (dosun&&(skyclearness>1))
714 {
715 printf("\nvoid light solar\n");
716 printf("0\n0\n");
717 printf("3 %.3e %.3e %.3e\n", solarradiance, solarradiance, solarradiance);
718 printf("\nsolar source sun\n");
719 printf("0\n0\n");
720 printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
721 } else if (dosun) {
722 printf("\nvoid light solar\n");
723 printf("0\n0\n");
724 printf("3 0.0 0.0 0.0\n");
725 printf("\nsolar source sun\n");
726 printf("0\n0\n");
727 printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
728 }
729 /* print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018 */
730 if (color_output==1 && skyclearness < 4.5 && skyclearness >1.065 )
731 {
732 fprintf(stderr, " warning: sky clearness(epsilon)= %f \n",skyclearness);
733 fprintf(stderr, " warning: intermediate sky!! \n");
734 fprintf(stderr, " warning: color model for intermediate sky pending \n");
735 fprintf(stderr, " warning: no color output ! \n");
736 color_output=0;
737 }
738 if (color_output==1)
739 {
740 printf("\nvoid colorfunc skyfunc\n");
741 printf("4 skybright_r skybright_g skybright_b perezlum_c.cal\n");
742 printf("0\n");
743 printf("22 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n", diffnormalization, groundbr,
744 *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
745 sundir[0], sundir[1], sundir[2],skyclearness,locus[0],locus[1],locus[2],locus[3],locus[4],locus[5],locus[6],locus[7],locus[8],locus[9],locus[10]);
746 }else{
747 printf("\nvoid brightfunc skyfunc\n");
748 printf("2 skybright perezlum.cal\n");
749 printf("0\n");
750 printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr,
751 *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
752 sundir[0], sundir[1], sundir[2]);
753 }
754
755 }
756
757
758
759 void print_error_sky()
760 {
761
762
763 sundir[0] = -sin(azimuth)*cos(altitude);
764 sundir[1] = -cos(azimuth)*cos(altitude);
765 sundir[2] = sin(altitude);
766
767 printf("# Local solar time: %.2f\n", st);
768 printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI);
769
770 printf("\nvoid brightfunc skyfunc\n");
771 printf("2 skybright perezlum.cal\n");
772 printf("0\n");
773 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]);
774 }
775
776
777
778
779
780 void printdefaults() /* print default values */
781 {
782 printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
783 if (zenithbr > 0.0)
784 printf("-b %f\t\t\t# Zenith radiance (watts/ster/m^2\n", zenithbr);
785 else
786 printf("-t %f\t\t\t# Atmospheric betaturbidity\n", betaturbidity);
787 printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/M_PI));
788 printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/M_PI));
789 printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/M_PI));
790 }
791
792
793
794
795 void usage_error(char* msg) /* print usage error and quit */
796 {
797 if (msg != NULL)
798 fprintf(stderr, "%s: Use error - %s\n\n", progname, msg);
799 fprintf(stderr, "Usage: %s month day hour [-y year] [...]\n", progname);
800 fprintf(stderr, " or: %s -ang altitude azimuth [...]\n", progname);
801 fprintf(stderr, " followed by: -P epsilon delta [options]\n");
802 fprintf(stderr, " or: [-W|-L|-G] direct_value diffuse_value [options]\n");
803 fprintf(stderr, " or: -E global_irradiance [options]\n\n");
804 fprintf(stderr, " Description:\n");
805 fprintf(stderr, " -P epsilon delta (these are the Perez parameters) \n");
806 fprintf(stderr, " -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
807 fprintf(stderr, " -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n");
808 fprintf(stderr, " -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
809 fprintf(stderr, " -E global-horizontal-irradiance (W/m^2)\n\n");
810 fprintf(stderr, " Output specification with option:\n");
811 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");
812 fprintf(stderr, " gendaylit version 2.5 (2018/04/18) \n\n");
813 exit(1);
814 }
815
816
817
818
819 double normsc() /* compute normalization factor (E0*F2/L0) */
820 {
821 static double nfc[2][5] = {
822 /* clear sky approx. */
823 {2.766521, 0.547665, -0.369832, 0.009237, 0.059229},
824 /* intermediate sky approx. */
825 {3.5556, -2.7152, -1.3081, 1.0660, 0.60227},
826 };
827 register double *nf;
828 double x, nsc;
829 register int i;
830 /* polynomial approximation */
831 nf = nfc[S_INTER];
832 x = (altitude - M_PI/4.0)/(M_PI/4.0);
833 nsc = nf[i=4];
834 while (i--)
835 nsc = nsc*x + nf[i];
836
837 return(nsc);
838 }
839
840
841
842
843
844 void printhead(int ac, char** av) /* print command header */
845 {
846 putchar('#');
847 while (ac--) {
848 putchar(' ');
849 fputs(*av++, stdout);
850 }
851 putchar('\n');
852 }
853
854
855
856
857
858
859 /* Perez models */
860
861 /* Perez global horizontal luminous efficacy model */
862 double glob_h_effi_PEREZ()
863 {
864
865 double value;
866 double category_bounds[10], a[10], b[10], c[10], d[10];
867 int category_total_number, category_number, i;
868
869 check_parametrization();
870
871
872 /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
873 fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n"); */
874
875
876 /* initialize category bounds (clearness index bounds) */
877
878 category_total_number = 8;
879
880 category_bounds[1] = 1;
881 category_bounds[2] = 1.065;
882 category_bounds[3] = 1.230;
883 category_bounds[4] = 1.500;
884 category_bounds[5] = 1.950;
885 category_bounds[6] = 2.800;
886 category_bounds[7] = 4.500;
887 category_bounds[8] = 6.200;
888 category_bounds[9] = 12.01;
889
890
891 /* initialize model coefficients */
892 a[1] = 96.63;
893 a[2] = 107.54;
894 a[3] = 98.73;
895 a[4] = 92.72;
896 a[5] = 86.73;
897 a[6] = 88.34;
898 a[7] = 78.63;
899 a[8] = 99.65;
900
901 b[1] = -0.47;
902 b[2] = 0.79;
903 b[3] = 0.70;
904 b[4] = 0.56;
905 b[5] = 0.98;
906 b[6] = 1.39;
907 b[7] = 1.47;
908 b[8] = 1.86;
909
910 c[1] = 11.50;
911 c[2] = 1.79;
912 c[3] = 4.40;
913 c[4] = 8.36;
914 c[5] = 7.10;
915 c[6] = 6.06;
916 c[7] = 4.93;
917 c[8] = -4.46;
918
919 d[1] = -9.16;
920 d[2] = -1.19;
921 d[3] = -6.95;
922 d[4] = -8.31;
923 d[5] = -10.94;
924 d[6] = -7.60;
925 d[7] = -11.37;
926 d[8] = -3.15;
927
928
929
930 for (i=1; i<=category_total_number; i++)
931 {
932 if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
933 category_number = i;
934 }
935
936 value = a[category_number] + b[category_number]*atm_preci_water +
937 c[category_number]*cos(sunzenith*M_PI/180) + d[category_number]*log(skybrightness);
938
939 return(value);
940 }
941
942
943
944
945 /* global horizontal diffuse efficacy model, according to PEREZ */
946 double glob_h_diffuse_effi_PEREZ()
947 {
948 double value;
949 double category_bounds[10], a[10], b[10], c[10], d[10];
950 int category_total_number, category_number, i;
951
952 check_parametrization();
953
954
955 /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
956 fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */
957
958 /* initialize category bounds (clearness index bounds) */
959
960 category_total_number = 8;
961
962 //XXX: category_bounds > 0.1
963 category_bounds[1] = 1;
964 category_bounds[2] = 1.065;
965 category_bounds[3] = 1.230;
966 category_bounds[4] = 1.500;
967 category_bounds[5] = 1.950;
968 category_bounds[6] = 2.800;
969 category_bounds[7] = 4.500;
970 category_bounds[8] = 6.200;
971 category_bounds[9] = 12.01;
972
973
974 /* initialize model coefficients */
975 a[1] = 97.24;
976 a[2] = 107.22;
977 a[3] = 104.97;
978 a[4] = 102.39;
979 a[5] = 100.71;
980 a[6] = 106.42;
981 a[7] = 141.88;
982 a[8] = 152.23;
983
984 b[1] = -0.46;
985 b[2] = 1.15;
986 b[3] = 2.96;
987 b[4] = 5.59;
988 b[5] = 5.94;
989 b[6] = 3.83;
990 b[7] = 1.90;
991 b[8] = 0.35;
992
993 c[1] = 12.00;
994 c[2] = 0.59;
995 c[3] = -5.53;
996 c[4] = -13.95;
997 c[5] = -22.75;
998 c[6] = -36.15;
999 c[7] = -53.24;
1000 c[8] = -45.27;
1001
1002 d[1] = -8.91;
1003 d[2] = -3.95;
1004 d[3] = -8.77;
1005 d[4] = -13.90;
1006 d[5] = -23.74;
1007 d[6] = -28.83;
1008 d[7] = -14.03;
1009 d[8] = -7.98;
1010
1011
1012
1013 category_number = -1;
1014 for (i=1; i<=category_total_number; i++)
1015 {
1016 if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
1017 category_number = i;
1018 }
1019
1020 if (category_number == -1) {
1021 if (suppress_warnings==0)
1022 fprintf(stderr, "Warning: sky clearness (= %.3f) too high, printing error sky\n", skyclearness);
1023 print_error_sky();
1024 exit(0);
1025 }
1026
1027
1028 value = a[category_number] + b[category_number]*atm_preci_water + c[category_number]*cos(sunzenith*M_PI/180) +
1029 d[category_number]*log(skybrightness);
1030
1031 return(value);
1032
1033 }
1034
1035
1036
1037
1038
1039
1040 /* direct normal efficacy model, according to PEREZ */
1041
1042 double direct_n_effi_PEREZ()
1043
1044 {
1045 double value;
1046 double category_bounds[10], a[10], b[10], c[10], d[10];
1047 int category_total_number, category_number, i;
1048
1049
1050 /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
1051 fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");*/
1052
1053
1054 /* initialize category bounds (clearness index bounds) */
1055
1056 category_total_number = 8;
1057
1058 category_bounds[1] = 1;
1059 category_bounds[2] = 1.065;
1060 category_bounds[3] = 1.230;
1061 category_bounds[4] = 1.500;
1062 category_bounds[5] = 1.950;
1063 category_bounds[6] = 2.800;
1064 category_bounds[7] = 4.500;
1065 category_bounds[8] = 6.200;
1066 category_bounds[9] = 12.1;
1067
1068
1069 /* initialize model coefficients */
1070 a[1] = 57.20;
1071 a[2] = 98.99;
1072 a[3] = 109.83;
1073 a[4] = 110.34;
1074 a[5] = 106.36;
1075 a[6] = 107.19;
1076 a[7] = 105.75;
1077 a[8] = 101.18;
1078
1079 b[1] = -4.55;
1080 b[2] = -3.46;
1081 b[3] = -4.90;
1082 b[4] = -5.84;
1083 b[5] = -3.97;
1084 b[6] = -1.25;
1085 b[7] = 0.77;
1086 b[8] = 1.58;
1087
1088 c[1] = -2.98;
1089 c[2] = -1.21;
1090 c[3] = -1.71;
1091 c[4] = -1.99;
1092 c[5] = -1.75;
1093 c[6] = -1.51;
1094 c[7] = -1.26;
1095 c[8] = -1.10;
1096
1097 d[1] = 117.12;
1098 d[2] = 12.38;
1099 d[3] = -8.81;
1100 d[4] = -4.56;
1101 d[5] = -6.16;
1102 d[6] = -26.73;
1103 d[7] = -34.44;
1104 d[8] = -8.29;
1105
1106
1107
1108 for (i=1; i<=category_total_number; i++)
1109 {
1110 if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
1111 category_number = i;
1112 }
1113
1114 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;
1115
1116 if (value < 0) value = 0;
1117
1118 return(value);
1119 }
1120
1121
1122 /*check the range of epsilon and delta indexes of the perez parametrization*/
1123 void check_parametrization()
1124 {
1125
1126 if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup)
1127 {
1128
1129 /* limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */
1130
1131 if (skyclearness<skyclearinf){
1132 /* if (suppress_warnings==0)
1133 fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */
1134 skyclearness=skyclearinf;
1135 }
1136 if (skyclearness>skyclearsup){
1137 /* if (suppress_warnings==0)
1138 fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */
1139 skyclearness=skyclearsup-0.001;
1140 }
1141 if (skybrightness<skybriginf){
1142 /* if (suppress_warnings==0)
1143 fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */
1144 skybrightness=skybriginf;
1145 }
1146 if (skybrightness>skybrigsup){
1147 /* if (suppress_warnings==0)
1148 fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */
1149 skybrightness=skybrigsup;
1150 }
1151
1152 return; }
1153 else return;
1154 }
1155
1156
1157
1158
1159
1160 /* validity of the direct and diffuse components */
1161 void check_illuminances()
1162 {
1163 if (directilluminance < 0) {
1164 if(suppress_warnings==0)
1165 { fprintf(stderr,"Warning: direct illuminance < 0. Using 0.0\n"); }
1166 directilluminance = 0.0;
1167 }
1168 if (diffuseilluminance < 0) {
1169 if(suppress_warnings==0)
1170 { fprintf(stderr,"Warning: diffuse illuminance < 0. Using 0.0\n"); }
1171 diffuseilluminance = 0.0;
1172 }
1173
1174 if (directilluminance+diffuseilluminance==0 && altitude > 0) {
1175 if(suppress_warnings==0)
1176 { fprintf(stderr,"Warning: zero illuminance at sun altitude > 0, printing error sky\n"); }
1177 print_error_sky();
1178 exit(0);
1179 }
1180
1181 if (directilluminance > solar_constant_l) {
1182 if(suppress_warnings==0)
1183 { fprintf(stderr,"Warning: direct illuminance exceeds solar constant\n"); }
1184 print_error_sky();
1185 exit(0);
1186 }
1187 }
1188
1189
1190 void check_irradiances()
1191 {
1192 if (directirradiance < 0) {
1193 if(suppress_warnings==0)
1194 { fprintf(stderr,"Warning: direct irradiance < 0. Using 0.0\n"); }
1195 directirradiance = 0.0;
1196 }
1197 if (diffuseirradiance < 0) {
1198 if(suppress_warnings==0)
1199 { fprintf(stderr,"Warning: diffuse irradiance < 0. Using 0.0\n"); }
1200 diffuseirradiance = 0.0;
1201 }
1202
1203 if (directirradiance+diffuseirradiance==0 && altitude > 0) {
1204 if(suppress_warnings==0)
1205 { fprintf(stderr,"Warning: zero irradiance at sun altitude > 0, printing error sky\n"); }
1206 print_error_sky();
1207 exit(0);
1208 }
1209
1210 if (directirradiance > solar_constant_e) {
1211 if(suppress_warnings==0)
1212 { fprintf(stderr,"Warning: direct irradiance exceeds solar constant\n"); }
1213 print_error_sky();
1214 exit(0);
1215 }
1216 }
1217
1218
1219
1220 /* Perez sky's brightness */
1221 double sky_brightness()
1222 {
1223 double value;
1224
1225 value = diffuseirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
1226
1227 return(value);
1228 }
1229
1230
1231 /* Perez sky's clearness */
1232 double sky_clearness()
1233 {
1234 double value;
1235
1236 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) ;
1237
1238 return(value);
1239 }
1240
1241
1242
1243 /* diffus horizontal irradiance from Perez sky's brightness */
1244 double diffuse_irradiance_from_sky_brightness()
1245 {
1246 double value;
1247
1248 value = skybrightness / air_mass() * ( solar_constant_e*get_eccentricity());
1249
1250 return(value);
1251 }
1252
1253
1254 /* direct normal irradiance from Perez sky's clearness */
1255 double direct_irradiance_from_sky_clearness()
1256 {
1257 double value;
1258
1259 value = diffuse_irradiance_from_sky_brightness();
1260 value = value * ( (skyclearness-1) * (1+1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) );
1261
1262 return(value);
1263 }
1264
1265
1266
1267
1268 void illu_to_irra_index()
1269 {
1270 double test1=0.1, test2=0.1, d_eff;
1271 int counter=0;
1272
1273 diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l);
1274 directirradiance = directilluminance*solar_constant_e/(solar_constant_l);
1275 skyclearness = sky_clearness();
1276 skybrightness = sky_brightness();
1277 check_parametrization();
1278
1279
1280 while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
1281 || (!(skyclearness<skyclearinf || skyclearness>skyclearsup))
1282 || (!(skybrightness<skybriginf || skybrightness>skybrigsup)) )
1283 && !(counter==9) )
1284 {
1285
1286 test1=diffuseirradiance;
1287 test2=directirradiance;
1288 counter++;
1289
1290 diffuseirradiance = diffuseilluminance/glob_h_diffuse_effi_PEREZ();
1291 d_eff = direct_n_effi_PEREZ();
1292
1293
1294 if (d_eff < 0.1)
1295 directirradiance = 0;
1296 else
1297 directirradiance = directilluminance/d_eff;
1298
1299 skybrightness = sky_brightness();
1300 skyclearness = sky_clearness();
1301 check_parametrization();
1302
1303 }
1304
1305
1306 return;
1307 }
1308
1309 static int get_numlin(float epsilon)
1310 {
1311 if (epsilon < 1.065)
1312 return 0;
1313 else if (epsilon < 1.230)
1314 return 1;
1315 else if (epsilon < 1.500)
1316 return 2;
1317 else if (epsilon < 1.950)
1318 return 3;
1319 else if (epsilon < 2.800)
1320 return 4;
1321 else if (epsilon < 4.500)
1322 return 5;
1323 else if (epsilon < 6.200)
1324 return 6;
1325 return 7;
1326 }
1327
1328 /* sky luminance perez model */
1329 double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[])
1330 {
1331
1332 float x[5][4];
1333 int i,j,num_lin;
1334 double c_perez[5];
1335
1336 if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) )
1337 {
1338 fprintf(stderr,"Error: epsilon out of range in function calc_rel_lum_perez!\n");
1339 exit(1);
1340 }
1341
1342 /* correction de modele de Perez solar energy ...*/
1343 if ( (epsilon > 1.065) && (epsilon < 2.8) )
1344 {
1345 if ( Delta < 0.2 ) Delta = 0.2;
1346 }
1347
1348
1349 num_lin = get_numlin(epsilon);
1350
1351 for (i=0;i<5;i++)
1352 for (j=0;j<4;j++)
1353 {
1354 x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j);
1355 /* fprintf(stderr,"x %d %d vaut %f\n",i,j,x[i][j]); */
1356 }
1357
1358
1359 if (num_lin)
1360 {
1361 for (i=0;i<5;i++)
1362 c_perez[i] = x[i][0] + x[i][1]*Z + Delta * (x[i][2] + x[i][3]*Z);
1363 }
1364 else
1365 {
1366 c_perez[0] = x[0][0] + x[0][1]*Z + Delta * (x[0][2] + x[0][3]*Z);
1367 c_perez[1] = x[1][0] + x[1][1]*Z + Delta * (x[1][2] + x[1][3]*Z);
1368 c_perez[4] = x[4][0] + x[4][1]*Z + Delta * (x[4][2] + x[4][3]*Z);
1369 c_perez[2] = exp( pow(Delta*(x[2][0]+x[2][1]*Z),x[2][2])) - x[2][3];
1370 c_perez[3] = -exp( Delta*(x[3][0]+x[3][1]*Z) )+x[3][2]+Delta*x[3][3];
1371 }
1372
1373
1374 return (1 + c_perez[0]*exp(c_perez[1]/cos(dzeta)) ) *
1375 (1 + c_perez[2]*exp(c_perez[3]*gamma) +
1376 c_perez[4]*cos(gamma)*cos(gamma) );
1377 }
1378
1379
1380
1381 /* coefficients for the sky luminance perez model */
1382 void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[])
1383 {
1384 float x[5][4];
1385 int i,j,num_lin;
1386
1387 if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) )
1388 {
1389 fprintf(stderr,"Error: epsilon out of range in function coeff_lum_perez!\n");
1390 exit(1);
1391 }
1392
1393 /* correction du modele de Perez solar energy ...*/
1394 if ( (epsilon > 1.065) && (epsilon < 2.8) )
1395 {
1396 if ( Delta < 0.2 ) Delta = 0.2;
1397 }
1398
1399
1400 num_lin = get_numlin(epsilon);
1401
1402 /*fprintf(stderr,"numlin %d\n", num_lin);*/
1403
1404 for (i=0;i<5;i++)
1405 for (j=0;j<4;j++)
1406 {
1407 x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j);
1408 /* printf("x %d %d vaut %f\n",i,j,x[i][j]); */
1409 }
1410
1411
1412 if (num_lin)
1413 {
1414 for (i=0;i<5;i++)
1415 *(c_perez+i) = x[i][0] + x[i][1]*Z + Delta * (x[i][2] + x[i][3]*Z);
1416
1417 }
1418 else
1419 {
1420 *(c_perez+0) = x[0][0] + x[0][1]*Z + Delta * (x[0][2] + x[0][3]*Z);
1421 *(c_perez+1) = x[1][0] + x[1][1]*Z + Delta * (x[1][2] + x[1][3]*Z);
1422 *(c_perez+4) = x[4][0] + x[4][1]*Z + Delta * (x[4][2] + x[4][3]*Z);
1423 *(c_perez+2) = exp( pow(Delta*(x[2][0]+x[2][1]*Z),x[2][2])) - x[2][3];
1424 *(c_perez+3) = -exp( Delta*(x[3][0]+x[3][1]*Z) )+x[3][2]+Delta*x[3][3];
1425
1426
1427 }
1428
1429
1430 return;
1431 }
1432
1433
1434
1435 /* degrees into radians */
1436 double radians(double degres)
1437 {
1438 return degres*(M_PI/180.);
1439 }
1440
1441
1442 /* radian into degrees */
1443 double degres(double radians)
1444 {
1445 return radians*(180./M_PI);
1446 }
1447
1448
1449 /* calculation of the angles dzeta and gamma */
1450 void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z)
1451 {
1452 *dzeta = theta; /* dzeta = phi */
1453 if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1 && (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi) < 1.1 ) )
1454 *gamma = 0;
1455 else if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1.1 )
1456 {
1457 printf("error in calculation of gamma (angle between point and sun");
1458 exit(1);
1459 }
1460 else
1461 *gamma = acos(cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi));
1462 }
1463
1464
1465
1466 double integ_lv(float *lv,float *theta)
1467 {
1468 int i;
1469 double buffer=0.0;
1470
1471 for (i=0;i<145;i++)
1472 {
1473 buffer += (*(lv+i))*cos(radians(*(theta+i)));
1474 }
1475
1476 return buffer*(2.*M_PI/145.);
1477 }
1478
1479
1480
1481 /* enter day number(double), return E0 = square(R0/R): eccentricity correction factor */
1482
1483 double get_eccentricity()
1484 {
1485 double day_angle;
1486 double E0;
1487
1488 day_angle = 2*M_PI*(daynumber - 1)/365;
1489 E0 = 1.00011+0.034221*cos(day_angle)+0.00128*sin(day_angle)+
1490 0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle);
1491
1492 return (E0);
1493 }
1494
1495
1496 /* enter sunzenith angle (degrees) return relative air mass (double) */
1497 double air_mass()
1498 {
1499 double m;
1500 if (sunzenith>90)
1501 {
1502 if(suppress_warnings==0)
1503 { fprintf(stderr, "Warning: air mass has reached the maximal value\n"); }
1504 sunzenith=90;
1505 }
1506 m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1507 return(m);
1508 }
1509
1510
1511