ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gendaylit.c
Revision: 2.16
Committed: Wed Jul 30 17:30:27 2014 UTC (9 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad5R0, rad5R1, rad4R2, rad4R2P1
Changes since 2.15: +1 -7 lines
Log Message:
Improved use of sun.h

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