ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gendaylit.c
Revision: 2.19
Committed: Wed Dec 4 18:11:57 2019 UTC (4 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R3
Changes since 2.18: +5 -5 lines
Log Message: 
Added -y option to usage

File Contents

# User Rev Content
1 greg 2.9 /* Copyright (c) 1994,2006 *Fraunhofer Institut for Solar Energy Systems
2     * Heidenhofstr. 2, D-79110 Freiburg, Germany
3 greg 2.1 * *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 greg 2.17 * 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 greg 2.1 */
9    
10 greg 2.15 #define _USE_MATH_DEFINES
11 greg 2.1 #include <stdio.h>
12     #include <string.h>
13     #include <math.h>
14     #include <stdlib.h>
15    
16     #include "color.h"
17 greg 2.12 #include "sun.h"
18 greg 2.1 #include "paths.h"
19    
20 greg 2.9 #define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
21 greg 2.1
22     double normsc();
23    
24 greg 2.19 /*static char *rcsid="$Header: /cvs/radiance/ray/src/gen/gendaylit.c,v 2.18 2019/11/07 23:15:06 greg Exp $";*/
25 greg 2.9
26     float coeff_perez[] = {
27 greg 2.12 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 greg 2.17 /* 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 greg 2.1
61    
62    
63 greg 2.14 /* Perez sky parametrization: epsilon and delta calculations from the direct and diffuse irradiances */
64 greg 2.1 double sky_brightness();
65     double sky_clearness();
66    
67     /* calculation of the direct and diffuse components from the Perez parametrization */
68 greg 2.9 double diffuse_irradiance_from_sky_brightness();
69 greg 2.1 double direct_irradiance_from_sky_clearness();
70    
71 greg 2.12 /* 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 greg 2.1
74     double glob_h_effi_PEREZ();
75     double glob_h_diffuse_effi_PEREZ();
76     double direct_n_effi_PEREZ();
77 greg 2.12
78 greg 2.1 /*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 greg 2.9 void print_error_sky();
84 greg 2.1
85 greg 2.12 double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]);
86 greg 2.9 void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]);
87 greg 2.1 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 greg 2.9 void printdefaults();
93 greg 2.12 void check_sun_position();
94 greg 2.9 void computesky();
95     void printhead(int ac, char** av);
96 greg 2.14 void usage_error(char* msg);
97 greg 2.9 void printsky();
98 greg 2.1
99 greg 2.9 FILE * frlibopen(char* fname);
100 greg 2.1
101     /* astronomy and geometry*/
102     double get_eccentricity();
103     double air_mass();
104    
105 greg 2.14 double solar_sunset(int month, int day);
106     double solar_sunrise(int month, int day);
107 greg 2.12
108 greg 2.1 const double AU = 149597890E3;
109     const double solar_constant_e = 1367; /* solar constant W/m^2 */
110 greg 2.14 const double solar_constant_l = 127500; /* solar constant lux */
111 greg 2.1
112     const double half_sun_angle = 0.2665;
113     const double half_direct_angle = 2.85;
114    
115 greg 2.14 const double skyclearinf = 1.0; /* limitations for the variation of the Perez parameters */
116     const double skyclearsup = 12.01;
117 greg 2.1 const double skybriginf = 0.01;
118     const double skybrigsup = 0.6;
119    
120    
121    
122     /* required values */
123 greg 2.18 int year = 0; /* year (optional) */
124 greg 2.1 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 greg 2.9 double skyclearness = 0;
133     double skybrightness = 0;
134 greg 2.12 double solarradiance;
135     double diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance;
136     double sunzenith, daynumber, atm_preci_water=2;
137 greg 2.1
138 greg 2.12 /*double sunaltitude_border = 0;*/
139 greg 2.9 double diffnormalization = 0;
140 greg 2.12 double dirnormalization = 0;
141 greg 2.1 double *c_perez;
142    
143 greg 2.12 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 greg 2.17 int color_output=0;
147 greg 2.9 int suppress_warnings=0;
148    
149 greg 2.1 /* default values */
150 greg 2.12 int cloudy = 0; /* 1=standard, 2=uniform */
151     int dosun = 1;
152 greg 2.1 double zenithbr = -1.0;
153     double betaturbidity = 0.1;
154     double gprefl = 0.2;
155     int S_INTER=0;
156    
157 greg 2.17
158 greg 2.1 /* computed values */
159     double sundir[3];
160 greg 2.9 double groundbr = 0;
161 greg 2.1 double F2;
162     double solarbr = 0.0;
163     int u_solar = 0; /* -1=irradiance, 1=radiance */
164 greg 2.12 float timeinterval = 0;
165    
166     char *progname;
167     char errmsg[128];
168    
169 greg 2.14 double st;
170 greg 2.1
171    
172 greg 2.9 int main(int argc, char** argv)
173 greg 2.1 {
174     int i;
175    
176     progname = argv[0];
177     if (argc == 2 && !strcmp(argv[1], "-defaults")) {
178     printdefaults();
179 greg 2.9 return 0;
180 greg 2.1 }
181     if (argc < 4)
182 greg 2.14 usage_error("arg count");
183 greg 2.1 if (!strcmp(argv[1], "-ang")) {
184 greg 2.9 altitude = atof(argv[2]) * (M_PI/180);
185     azimuth = atof(argv[3]) * (M_PI/180);
186 greg 2.1 month = 0;
187     } else {
188     month = atoi(argv[1]);
189     if (month < 1 || month > 12)
190 greg 2.14 usage_error("bad month");
191 greg 2.1 day = atoi(argv[2]);
192     if (day < 1 || day > 31)
193 greg 2.14 usage_error("bad day");
194 greg 2.1 hour = atof(argv[3]);
195     if (hour < 0 || hour >= 24)
196 greg 2.14 usage_error("bad hour");
197 greg 2.1 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 greg 2.18 case 'y':
207     year = atoi(argv[++i]);
208     break;
209 greg 2.1 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 greg 2.17 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 greg 2.1 case 't':
248     betaturbidity = atof(argv[++i]);
249     break;
250 greg 2.9 case 'w':
251     suppress_warnings = 1;
252     break;
253 greg 2.1 case 'b':
254     zenithbr = atof(argv[++i]);
255     break;
256     case 'g':
257     gprefl = atof(argv[++i]);
258     break;
259     case 'a':
260 greg 2.9 s_latitude = atof(argv[++i]) * (M_PI/180);
261 greg 2.1 break;
262     case 'o':
263 greg 2.9 s_longitude = atof(argv[++i]) * (M_PI/180);
264 greg 2.1 break;
265     case 'm':
266 greg 2.9 s_meridian = atof(argv[++i]) * (M_PI/180);
267 greg 2.1 break;
268    
269     case 'O':
270 greg 2.14 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 greg 2.1 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 greg 2.9 diffuseirradiance = atof(argv[++i]);
285 greg 2.1 break;
286    
287     case 'L': /* direct normal Illuminance [Lux] */
288     input = 2; /* diffuse horizontal Ill. [Lux] */
289     directilluminance = atof(argv[++i]);
290 greg 2.9 diffuseilluminance = atof(argv[++i]);
291 greg 2.1 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 greg 2.9 diffuseirradiance = atof(argv[++i]);
297     break;
298    
299 greg 2.12 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 greg 2.9
308 greg 2.1
309     default:
310     sprintf(errmsg, "unknown option: %s", argv[i]);
311 greg 2.14 usage_error(errmsg);
312 greg 2.1 }
313     else
314 greg 2.14 usage_error("bad option");
315 greg 2.1
316 greg 2.14 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 greg 2.9 progname, (s_longitude-s_meridian)*12/M_PI);
319 greg 2.1
320    
321 greg 2.12 /* dynamic memory allocation for the pointers */
322 greg 2.9 if ( (c_perez = calloc(5, sizeof(double))) == NULL )
323 greg 2.14 { fprintf(stderr,"Out of memory error in function main"); return 1; }
324 greg 2.1
325 greg 2.14
326 greg 2.1 printhead(argc, argv);
327     computesky();
328     printsky();
329 greg 2.9 return 0;
330 greg 2.14
331 greg 2.1 }
332    
333    
334 greg 2.12
335    
336    
337     void computesky()
338 greg 2.1 {
339    
340 greg 2.6 int j;
341 greg 2.12
342     float *lv_mod; /* 145 luminance values */
343 greg 2.9 float *theta_o, *phi_o;
344 greg 2.1 double dzeta, gamma;
345     double normfactor;
346 greg 2.12 double erbs_s0, erbs_kt;
347 greg 2.1
348    
349     /* compute solar direction */
350 greg 2.12
351 greg 2.1 if (month) { /* from date and time */
352 greg 2.14 double sd;
353 greg 2.1
354 greg 2.18 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 greg 2.12 if(timeinterval) {
369    
370     if(timeinterval<0) {
371     fprintf(stderr, "time interval negative\n");
372     exit(1);
373     }
374 greg 2.14
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 greg 2.12 }
380    
381 greg 2.14 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 greg 2.12 }
386 greg 2.14
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 greg 2.12 }
396 greg 2.14 else
397 greg 2.12
398 greg 2.14 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 greg 2.12
407 greg 2.1 altitude = salt(sd, st);
408     azimuth = sazi(sd, st);
409    
410     daynumber = (double)jdate(month, day);
411 greg 2.12
412 greg 2.1 }
413 greg 2.9
414    
415 greg 2.12
416    
417 greg 2.9
418     if (!cloudy && altitude > 87.*M_PI/180.) {
419    
420     if (suppress_warnings==0) {
421     fprintf(stderr,
422 greg 2.1 "%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n",
423     progname);
424 greg 2.9 }
425     altitude = 87.*M_PI/180.;
426 greg 2.1 }
427 greg 2.9
428 greg 2.14
429    
430 greg 2.1 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 greg 2.9 sunzenith = 90 - altitude*180/M_PI;
437 greg 2.14
438 greg 2.1
439 greg 2.9 /* compute the inputs for the calculation of the light distribution over the sky*/
440 greg 2.12 if (input==0) /* P */
441 greg 2.1 {
442     check_parametrization();
443 greg 2.9 diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
444 greg 2.1 directirradiance = direct_irradiance_from_sky_clearness();
445     check_irradiances();
446    
447     if (output==0 || output==2)
448     {
449 greg 2.9 diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
450 greg 2.1 directilluminance = directirradiance*direct_n_effi_PEREZ();
451     check_illuminances();
452     }
453     }
454    
455    
456 greg 2.12 else if (input==1) /* W */
457 greg 2.1 {
458     check_irradiances();
459     skybrightness = sky_brightness();
460     skyclearness = sky_clearness();
461 greg 2.14
462 greg 2.1 check_parametrization();
463 greg 2.14
464 greg 2.1 if (output==0 || output==2)
465     {
466 greg 2.9 diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
467 greg 2.1 directilluminance = directirradiance*direct_n_effi_PEREZ();
468     check_illuminances();
469     }
470    
471     }
472    
473    
474 greg 2.12 else if (input==2) /* L */
475 greg 2.1 {
476     check_illuminances();
477     illu_to_irra_index();
478     check_parametrization();
479     }
480    
481    
482 greg 2.12 else if (input==3) /* G */
483 greg 2.1 {
484     if (altitude<=0)
485     {
486 greg 2.9 if (suppress_warnings==0)
487 greg 2.14 fprintf(stderr, "Warning: sun altitude < 0, proceed with irradiance values of zero\n");
488 greg 2.9 directirradiance = 0;
489     diffuseirradiance = 0;
490     } else {
491 greg 2.12
492     directirradiance=directirradiance/sin(altitude);
493 greg 2.1 }
494 greg 2.12
495 greg 2.1 check_irradiances();
496     skybrightness = sky_brightness();
497     skyclearness = sky_clearness();
498     check_parametrization();
499    
500     if (output==0 || output==2)
501     {
502 greg 2.9 diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
503 greg 2.1 directilluminance = directirradiance*direct_n_effi_PEREZ();
504     check_illuminances();
505     }
506    
507     }
508    
509 greg 2.12
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 greg 2.14 printf("# WARNING: the -E option is only recommended for a rough estimation!\n");
540 greg 2.12
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 greg 2.1
562 greg 2.14 else { fprintf(stderr,"error at the input arguments"); exit(1); }
563 greg 2.1
564    
565    
566 greg 2.9 /* normalization factor for the relative sky luminance distribution, diffuse part*/
567 greg 2.12
568 greg 2.1 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 greg 2.9 theta_o = defangle_theta;
576     phi_o = defangle_phi;
577 greg 2.12
578 greg 2.1
579 greg 2.9 /* parameters for the perez model */
580 greg 2.1 coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez);
581    
582 greg 2.12
583    
584 greg 2.9 /*calculation of the modelled luminance */
585 greg 2.1 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 greg 2.12
589 greg 2.1 *(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez);
590 greg 2.12
591     /* fprintf(stderr,"theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); */
592 greg 2.1 }
593 greg 2.12
594 greg 2.1 /* integration of luminance for the normalization factor, diffuse part of the sky*/
595 greg 2.12
596 greg 2.1 diffnormalization = integ_lv(lv_mod, theta_o);
597    
598    
599    
600 greg 2.9 /*normalization coefficient in lumen or in watt*/
601 greg 2.1 if (output==0)
602     {
603 greg 2.9 diffnormalization = diffuseilluminance/diffnormalization/WHTEFFICACY;
604 greg 2.1 }
605     else if (output==1)
606     {
607 greg 2.9 diffnormalization = diffuseirradiance/diffnormalization;
608 greg 2.1 }
609     else if (output==2)
610     {
611 greg 2.9 diffnormalization = diffuseilluminance/diffnormalization;
612 greg 2.1 }
613    
614 greg 2.9 else {fprintf(stderr,"Wrong output specification.\n"); exit(1);}
615 greg 2.1
616    
617    
618    
619 greg 2.9 /* calculation for the solar source */
620 greg 2.1 if (output==0)
621 greg 2.9 solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY;
622 greg 2.1
623     else if (output==1)
624 greg 2.9 solarradiance = directirradiance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
625 greg 2.1
626     else
627 greg 2.9 solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
628 greg 2.1
629    
630    
631 greg 2.14 /* Compute the ground radiance */
632     zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
633     zenithbr*=diffnormalization;
634 greg 2.9
635 greg 2.14 if (skyclearness==1)
636 greg 2.1 normfactor = 0.777778;
637    
638 greg 2.14 if (skyclearness>=6)
639 greg 2.1 {
640 greg 2.9 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 greg 2.1 }
643    
644 greg 2.14 if ( (skyclearness>1) && (skyclearness<6) )
645 greg 2.1 {
646     S_INTER=1;
647 greg 2.9 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 greg 2.1 }
650    
651 greg 2.14 groundbr = zenithbr*normfactor;
652 greg 2.1
653 greg 2.14 if (dosun&&(skyclearness>1))
654 greg 2.9 groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];
655 greg 2.1
656 greg 2.14 groundbr *= gprefl;
657 greg 2.1
658    
659 greg 2.14
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 greg 2.1
669     return;
670     }
671    
672    
673    
674    
675 greg 2.12
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 greg 2.14
686    
687 greg 2.12 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 greg 2.14 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 greg 2.12
704    
705 greg 2.1 if (dosun&&(skyclearness>1))
706 greg 2.9 {
707 greg 2.1 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 greg 2.9 } else if (dosun) {
714 greg 2.1 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 greg 2.9 }
721 greg 2.17 /* 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 greg 2.1 printf("0\n");
735 greg 2.17 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 greg 2.9 *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
737 greg 2.17 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 greg 2.12
747 greg 2.1 }
748    
749    
750 greg 2.14
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 greg 2.9 void printdefaults() /* print default values */
773 greg 2.1 {
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 greg 2.9 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 greg 2.1 }
783    
784    
785 greg 2.14
786    
787     void usage_error(char* msg) /* print usage error and quit */
788 greg 2.1 {
789     if (msg != NULL)
790 greg 2.12 fprintf(stderr, "%s: Use error - %s\n\n", progname, msg);
791 greg 2.19 fprintf(stderr, "Usage: %s month day hour [-y year] [...]\n", progname);
792     fprintf(stderr, " or: %s -ang altitude azimuth [...]\n", progname);
793 greg 2.12 fprintf(stderr, " followed by: -P epsilon delta [options]\n");
794     fprintf(stderr, " or: [-W|-L|-G] direct_value diffuse_value [options]\n");
795 greg 2.19 fprintf(stderr, " or: -E global_irradiance [options]\n\n");
796     fprintf(stderr, " Description:\n");
797 greg 2.1 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 greg 2.12 fprintf(stderr, " -E global-horizontal-irradiance (W/m^2)\n\n");
802     fprintf(stderr, " Output specification with option:\n");
803 greg 2.1 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 greg 2.17 fprintf(stderr, " gendaylit version 2.5 (2018/04/18) \n\n");
805 greg 2.1 exit(1);
806     }
807    
808    
809    
810 greg 2.14
811 greg 2.9 double normsc() /* compute normalization factor (E0*F2/L0) */
812 greg 2.1 {
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 greg 2.9 x = (altitude - M_PI/4.0)/(M_PI/4.0);
825 greg 2.1 nsc = nf[i=4];
826     while (i--)
827     nsc = nsc*x + nf[i];
828    
829     return(nsc);
830     }
831    
832    
833    
834 greg 2.14
835    
836 greg 2.9 void printhead(int ac, char** av) /* print command header */
837 greg 2.1 {
838     putchar('#');
839     while (ac--) {
840     putchar(' ');
841     fputs(*av++, stdout);
842     }
843     putchar('\n');
844     }
845    
846    
847    
848    
849 greg 2.14
850    
851 greg 2.1 /* 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 greg 2.12
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 greg 2.1 /* 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 greg 2.9 c[category_number]*cos(sunzenith*M_PI/180) + d[category_number]*log(skybrightness);
930 greg 2.1
931     return(value);
932     }
933    
934    
935 greg 2.14
936    
937 greg 2.1 /* 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 greg 2.12 check_parametrization();
945 greg 2.1
946 greg 2.12
947     /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
948 greg 2.14 fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */
949 greg 2.12
950 greg 2.1 /* initialize category bounds (clearness index bounds) */
951    
952     category_total_number = 8;
953    
954 greg 2.12 //XXX: category_bounds > 0.1
955 greg 2.1 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 greg 2.9 category_number = -1;
1006 greg 2.1 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 greg 2.9 if (category_number == -1) {
1013     if (suppress_warnings==0)
1014 greg 2.14 fprintf(stderr, "Warning: sky clearness (= %.3f) too high, printing error sky\n", skyclearness);
1015 greg 2.9 print_error_sky();
1016 greg 2.14 exit(0);
1017 greg 2.9 }
1018    
1019    
1020     value = a[category_number] + b[category_number]*atm_preci_water + c[category_number]*cos(sunzenith*M_PI/180) +
1021 greg 2.1 d[category_number]*log(skybrightness);
1022    
1023     return(value);
1024 greg 2.12
1025 greg 2.1 }
1026    
1027    
1028 greg 2.12
1029 greg 2.14
1030    
1031    
1032 greg 2.1 /* 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 greg 2.14 /*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 greg 2.1
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 greg 2.9 value = a[category_number] + b[category_number]*atm_preci_water + c[category_number]*exp(5.73*sunzenith*M_PI/180 - 5) + d[category_number]*skybrightness;
1107 greg 2.1
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 greg 2.14
1118 greg 2.9 if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup)
1119 greg 2.1 {
1120 greg 2.9
1121     /* limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */
1122 greg 2.14
1123 greg 2.9 if (skyclearness<skyclearinf){
1124 greg 2.14 /* if (suppress_warnings==0)
1125     fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */
1126 greg 2.9 skyclearness=skyclearinf;
1127     }
1128     if (skyclearness>skyclearsup){
1129 greg 2.14 /* if (suppress_warnings==0)
1130     fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */
1131     skyclearness=skyclearsup-0.001;
1132 greg 2.9 }
1133     if (skybrightness<skybriginf){
1134 greg 2.14 /* if (suppress_warnings==0)
1135     fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */
1136 greg 2.9 skybrightness=skybriginf;
1137     }
1138     if (skybrightness>skybrigsup){
1139 greg 2.14 /* if (suppress_warnings==0)
1140     fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */
1141 greg 2.9 skybrightness=skybrigsup;
1142 greg 2.1 }
1143 greg 2.9
1144     return; }
1145 greg 2.1 else return;
1146     }
1147    
1148    
1149 greg 2.14
1150    
1151    
1152 greg 2.9 /* validity of the direct and diffuse components */
1153 greg 2.1 void check_illuminances()
1154     {
1155 greg 2.9 if (directilluminance < 0) {
1156 greg 2.14 if(suppress_warnings==0)
1157     { fprintf(stderr,"Warning: direct illuminance < 0. Using 0.0\n"); }
1158 greg 2.9 directilluminance = 0.0;
1159     }
1160     if (diffuseilluminance < 0) {
1161 greg 2.14 if(suppress_warnings==0)
1162     { fprintf(stderr,"Warning: diffuse illuminance < 0. Using 0.0\n"); }
1163 greg 2.9 diffuseilluminance = 0.0;
1164     }
1165 greg 2.14
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 greg 2.1 }
1179     }
1180    
1181    
1182     void check_irradiances()
1183     {
1184 greg 2.9 if (directirradiance < 0) {
1185 greg 2.14 if(suppress_warnings==0)
1186     { fprintf(stderr,"Warning: direct irradiance < 0. Using 0.0\n"); }
1187 greg 2.9 directirradiance = 0.0;
1188     }
1189     if (diffuseirradiance < 0) {
1190 greg 2.14 if(suppress_warnings==0)
1191     { fprintf(stderr,"Warning: diffuse irradiance < 0. Using 0.0\n"); }
1192 greg 2.9 diffuseirradiance = 0.0;
1193     }
1194 greg 2.14
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 greg 2.9 if (directirradiance > solar_constant_e) {
1203 greg 2.14 if(suppress_warnings==0)
1204     { fprintf(stderr,"Warning: direct irradiance exceeds solar constant\n"); }
1205     print_error_sky();
1206     exit(0);
1207 greg 2.9 }
1208 greg 2.1 }
1209    
1210    
1211    
1212     /* Perez sky's brightness */
1213     double sky_brightness()
1214     {
1215     double value;
1216    
1217 greg 2.9 value = diffuseirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
1218 greg 2.1
1219     return(value);
1220     }
1221    
1222    
1223     /* Perez sky's clearness */
1224     double sky_clearness()
1225     {
1226 greg 2.9 double value;
1227 greg 2.1
1228 greg 2.9 value = ( (diffuseirradiance + directirradiance)/(diffuseirradiance) + 1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180 ) / (1 + 1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) ;
1229 greg 2.1
1230 greg 2.9 return(value);
1231 greg 2.1 }
1232    
1233    
1234    
1235     /* diffus horizontal irradiance from Perez sky's brightness */
1236 greg 2.9 double diffuse_irradiance_from_sky_brightness()
1237 greg 2.1 {
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 greg 2.9 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 greg 2.1
1254     return(value);
1255     }
1256    
1257    
1258 greg 2.12
1259    
1260 greg 2.9 void illu_to_irra_index()
1261 greg 2.1 {
1262 greg 2.9 double test1=0.1, test2=0.1, d_eff;
1263 greg 2.1 int counter=0;
1264    
1265 greg 2.14 diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l);
1266     directirradiance = directilluminance*solar_constant_e/(solar_constant_l);
1267 greg 2.1 skyclearness = sky_clearness();
1268     skybrightness = sky_brightness();
1269 greg 2.9 check_parametrization();
1270 greg 2.12
1271    
1272 greg 2.9 while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
1273 greg 2.12 || (!(skyclearness<skyclearinf || skyclearness>skyclearsup))
1274     || (!(skybrightness<skybriginf || skybrightness>skybrigsup)) )
1275     && !(counter==9) )
1276 greg 2.1 {
1277 greg 2.12
1278 greg 2.9 test1=diffuseirradiance;
1279 greg 2.1 test2=directirradiance;
1280     counter++;
1281    
1282 greg 2.9 diffuseirradiance = diffuseilluminance/glob_h_diffuse_effi_PEREZ();
1283     d_eff = direct_n_effi_PEREZ();
1284 greg 2.12
1285    
1286 greg 2.9 if (d_eff < 0.1)
1287     directirradiance = 0;
1288 greg 2.12 else
1289 greg 2.9 directirradiance = directilluminance/d_eff;
1290 greg 2.1
1291     skybrightness = sky_brightness();
1292     skyclearness = sky_clearness();
1293 greg 2.9 check_parametrization();
1294 greg 2.12
1295 greg 2.1 }
1296    
1297    
1298     return;
1299     }
1300    
1301 greg 2.9 static int get_numlin(float epsilon)
1302 greg 2.1 {
1303 greg 2.9 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 greg 2.1 }
1319    
1320     /* sky luminance perez model */
1321 greg 2.9 double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[])
1322 greg 2.1 {
1323 greg 2.12
1324 greg 2.1 float x[5][4];
1325     int i,j,num_lin;
1326     double c_perez[5];
1327    
1328     if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) )
1329     {
1330 greg 2.14 fprintf(stderr,"Error: epsilon out of range in function calc_rel_lum_perez!\n");
1331 greg 2.1 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 greg 2.12
1340    
1341 greg 2.9 num_lin = get_numlin(epsilon);
1342 greg 2.12
1343 greg 2.1 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 greg 2.12 /* fprintf(stderr,"x %d %d vaut %f\n",i,j,x[i][j]); */
1348 greg 2.1 }
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 greg 2.9 void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[])
1375 greg 2.1 {
1376     float x[5][4];
1377     int i,j,num_lin;
1378    
1379     if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) )
1380     {
1381 greg 2.14 fprintf(stderr,"Error: epsilon out of range in function coeff_lum_perez!\n");
1382 greg 2.1 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 greg 2.12
1391    
1392 greg 2.9 num_lin = get_numlin(epsilon);
1393    
1394 greg 2.12 /*fprintf(stderr,"numlin %d\n", num_lin);*/
1395 greg 2.1
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 greg 2.12
1427 greg 2.1 /* degrees into radians */
1428     double radians(double degres)
1429     {
1430 greg 2.9 return degres*M_PI/180.0;
1431 greg 2.1 }
1432    
1433 greg 2.12
1434 greg 2.1 /* radian into degrees */
1435     double degres(double radians)
1436     {
1437 greg 2.9 return radians/M_PI*180.0;
1438 greg 2.1 }
1439    
1440 greg 2.12
1441 greg 2.1 /* 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 greg 2.14 exit(1);
1451 greg 2.1 }
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 greg 2.12
1463 greg 2.1 for (i=0;i<145;i++)
1464 greg 2.12 {
1465 greg 2.1 buffer += (*(lv+i))*cos(radians(*(theta+i)));
1466 greg 2.12 }
1467    
1468 greg 2.9 return buffer*2*M_PI/144;
1469 greg 2.1 }
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 greg 2.9 day_angle = 2*M_PI*(daynumber - 1)/365;
1481 greg 2.1 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 greg 2.14 if(suppress_warnings==0)
1495     { fprintf(stderr, "Warning: air mass has reached the maximal value\n"); }
1496     sunzenith=90;
1497 greg 2.1 }
1498 greg 2.9 m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1499 greg 2.1 return(m);
1500     }
1501    
1502    
1503