ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gendaylit.c
Revision: 2.12
Committed: Fri Aug 9 16:44:19 2013 UTC (10 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.11: +307 -152 lines
Log Message:
Bug fixes and new -E and -i options (from Wendelin Sprenger and Jan Wienold)

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