/[radiance]/ray/src/gen/gendaylit.c
ViewVC logotype

Contents of /ray/src/gen/gendaylit.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2.19 - (show annotations) (download)
Wed Dec 4 18:11:57 2019 UTC (10 days, 10 hours ago) by greg
Branch: MAIN
CVS Tags: HEAD
Changes since 2.18: +5 -5 lines
File MIME type: text/plain
Added -y option to usage

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

webmaster-btp@lbl.gov
ViewVC Help
Powered by ViewVC 1.1.27