ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/radiance/ray/src/gen/gendaylit.c

Comparing ray/src/gen/gendaylit.c (file contents):
Revision 2.5 by greg, Wed Aug 10 22:30:31 2011 UTC vs.
Revision 2.12 by greg, Fri Aug 9 16:44:19 2013 UTC

#	Line 1 | Line 1
1	<	#ifndef lint
2	<	static const char RCSid[] = "$Id$";
3	<	#endif
4	<	/*        Copyright (c) 1994    *Fraunhofer Institut for Solar Energy Systems
5	<	*                              Oltmannstr 5, D-79100 Freiburg, Germany
1	>	/*      Copyright (c) 1994,2006 *Fraunhofer Institut for Solar Energy Systems
2	>	*                              Heidenhofstr. 2, D-79110 Freiburg, Germany
3		*                              *Agence de l'Environnement et de la Maitrise de l'Energie
4		*                              Centre de Valbonne, 500 route des Lucioles, 06565 Sophia Antipolis Cedex, France
5		*                              *BOUYGUES
6		*                              1 Avenue Eugene Freyssinet, Saint-Quentin-Yvelines, France
7		*/
8
12	–
13	–
14	–	/*
15	–	*  gendaylit.c         program to generate the angular distribution of the daylight.
16	–	*                      Our zenith is along the Z-axis, the X-axis
17	–	*                      points east, and the Y-axis points north.
18	–	*/
19	–
9		#include  <stdio.h>
10		#include  <string.h>
11		#include  <math.h>
12		#include  <stdlib.h>
24	–	#include  <ctype.h>
13
26	–	#include  "rtio.h"
27	–	#include  "fvect.h"
14		#include  "color.h"
15	+	#include  "sun.h"
16		#include  "paths.h"
17
18	<	extern int jdate(int month, int day);
32	<	extern double stadj(int  jd);
33	<	extern double sdec(int  jd);
34	<	extern double salt(double sd, double st);
35	<	extern double sazi(double sd, double st);
18	>	#define  DOT(v1,v2)     (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
19
20		double  normsc();
21
22	<	#define DATFILE         "coeff_perez.dat"
22	>	/*static        char *rcsid="$Header$";*/
23
24	+	float coeff_perez[] = {
25	+	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	+
56	+
57	+
58		/* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */
59		double sky_brightness();
60		double sky_clearness();
46	–	void computesky();
61
62		/* calculation of the direct and diffuse components from the Perez parametrization */
63	<	double  diffus_irradiance_from_sky_brightness();
63	>	double  diffuse_irradiance_from_sky_brightness();
64		double  direct_irradiance_from_sky_clearness();
65
66	+	/* 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
53	–	/* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */
69		double  glob_h_effi_PEREZ();
70		double  glob_h_diffuse_effi_PEREZ();
71		double  direct_n_effi_PEREZ();
72	+
73		/*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	+	void    print_error_sky();
79
80	<
81	<	/* Perez sky luminance model */
65	<	int     lect_coeff_perez(char *filename,float **coeff_perez);
66	<	double  calc_rel_lum_perez(double dzeta,double gamma,double Z,
67	<	double epsilon,double Delta,float *coeff_perez);
68	<	/* coefficients for the sky luminance perez model */
69	<	void    coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez);
80	>	double  calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]);
81	>	void    coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]);
82		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);
74	–	float   *theta_ordered(char *filename);
75	–	float   *phi_ordered(char *filename);
76	–	void    skip_comments(FILE *fp);
86
87	+	void printdefaults();
88	+	void check_sun_position();
89	+	void computesky();
90	+	void printhead(int ac, char** av);
91	+	void userror(char* msg);
92	+	void printsky();
93
94	+	FILE * frlibopen(char* fname);
95
96		/* astronomy and geometry*/
97		double  get_eccentricity();
98		double  air_mass();
83	–	double  get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut);
99
100	<
101	<	/* date*/
100	>	double  solar_sunset();
101	>	double  solar_sunrise();
102	>	double  stadj();
103		int     jdate(int month, int day);
104
105
106
91	–
92	–
107		/* sun calculation constants */
108	<	extern double  s_latitude;
109	<	extern double  s_longitude;
110	<	extern double  s_meridian;
108	>	extern double   s_latitude;
109	>	extern double   s_longitude;
110	>	extern double   s_meridian;
111
112		const double    AU = 149597890E3;
113		const double    solar_constant_e = 1367;    /* solar constant W/m^2 */
#	Line 102 | Line 116 | const double   solar_constant_l = 127.5;   /* solar co
116		const double    half_sun_angle = 0.2665;
117		const double    half_direct_angle = 2.85;
118
119	<	const double    skyclearinf = 1.000;    /* limitations for the variation of the Perez parameters */
119	>	const double    skyclearinf = 1.0;      /* limitations for the variation of the Perez parameters */
120		const double    skyclearsup = 12.1;
121		const double    skybriginf = 0.01;
122		const double    skybrigsup = 0.6;
#	Line 118 | Line 132 | double  altitude, azimuth;                     /* or solar angles */
132
133
134		/* definition of the sky conditions through the Perez parametrization */
135	<	double  skyclearness, skybrightness;
136	<	double  solarradiance;  /*radiance of the sun disk and of the circumsolar area*/
137	<	double  diffusilluminance, directilluminance, diffusirradiance, directirradiance;
138	<	double  sunzenith, daynumber=150, atm_preci_water=2;
135	>	double  skyclearness = 0;
136	>	double  skybrightness = 0;
137	>	double  solarradiance;
138	>	double  diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance;
139	>	double  sunzenith, daynumber, atm_preci_water=2;
140
141	<	double  diffnormalization, dirnormalization;
141	>	/*double  sunaltitude_border = 0;*/
142	>	double  diffnormalization = 0;
143	>	double  dirnormalization = 0;
144		double  *c_perez;
145
146	<	int     output=0;       /*define the unit of the output (sky luminance or radiance): visible watt=0, solar watt=1, lumen=2*/
147	<	int     input=0;        /*define the input for the calulation*/
146	>	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
150	+	int     suppress_warnings=0;
151	+
152		/* default values */
153	<	int  cloudy = 0;                                /* 1=standard, 2=uniform */
154	<	int  dosun = 1;
153	>	int     cloudy = 0;                             /* 1=standard, 2=uniform */
154	>	int     dosun = 1;
155		double  zenithbr = -1.0;
156		double  betaturbidity = 0.1;
157		double  gprefl = 0.2;
#	Line 139 | Line 159 | int    S_INTER=0;
159
160		/* computed values */
161		double  sundir[3];
162	<	double  groundbr;
162	>	double  groundbr = 0;
163		double  F2;
164		double  solarbr = 0.0;
165		int     u_solar = 0;                            /* -1=irradiance, 1=radiance */
166	+	float   timeinterval = 0;
167
168	<	char  *progname;
169	<	char  errmsg[128];
168	>	char    *progname;
169	>	char    errmsg[128];
170
171
172	<	main(argc, argv)
173	<	int  argc;
174	<	char  *argv[];
172	>
173	>
174	>	int main(int argc, char** argv)
175		{
176		int  i;
177
178		progname = argv[0];
179		if (argc == 2 && !strcmp(argv[1], "-defaults")) {
180		printdefaults();
181	<	exit(0);
181	>	return 0;
182		}
183		if (argc < 4)
184		userror("arg count");
#	Line 184 | Line 205 | char  *argv[];
205		cloudy = 0;
206		dosun = argv[i][0] == '+';
207		break;
187	–	case 'r':
208		case 'R':
209		u_solar = argv[i][1] == 'R' ? -1 : 1;
210		solarbr = atof(argv[++i]);
#	Line 196 | Line 216 | char  *argv[];
216		case 't':
217		betaturbidity = atof(argv[++i]);
218		break;
219	+	case 'w':
220	+	suppress_warnings = 1;
221	+	break;
222		case 'b':
223		zenithbr = atof(argv[++i]);
224		break;
#	Line 211 | Line 234 | char  *argv[];
234		case 'm':
235		s_meridian = atof(argv[++i]) * (M_PI/180);
236		break;
214	–
237
238		case 'O':
239		output = atof(argv[++i]);       /*define the unit of the output of the program :
240	<	sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm)
219	<	default is set to 0*/
240	>	sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm) */
241		break;
242
243		case 'P':
#	Line 228 | Line 249 | char  *argv[];
249		case 'W':                                       /* direct normal Irradiance [W/m^2] */
250		input = 1;                              /* diffuse horizontal Irrad. [W/m^2] */
251		directirradiance = atof(argv[++i]);
252	<	diffusirradiance = atof(argv[++i]);
252	>	diffuseirradiance = atof(argv[++i]);
253		break;
254
255		case 'L':                                       /* direct normal Illuminance [Lux] */
256		input = 2;                              /* diffuse horizontal Ill. [Lux] */
257		directilluminance = atof(argv[++i]);
258	<	diffusilluminance = atof(argv[++i]);
258	>	diffuseilluminance = atof(argv[++i]);
259		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	<	diffusirradiance = atof(argv[++i]);
264	>	diffuseirradiance = atof(argv[++i]);
265		break;
266	+
267	+	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	+	sunaltitude_border = atof(argv[++i]);
275	+	break;
276	+	*/
277
278	+	case 'i':
279	+	timeinterval = atof(argv[++i]);
280	+	break;
281
282	+
283		default:
284		sprintf(errmsg, "unknown option: %s", argv[i]);
285		userror(errmsg);
#	Line 257 | Line 293 | char  *argv[];
293		progname, (s_longitude-s_meridian)*12/M_PI);
294
295
296	<	/* allocation dynamique de memoire pour les pointeurs */
297	<	if ( (c_perez = malloc(5*sizeof(double))) == NULL )
296	>	/* dynamic memory allocation for the pointers */
297	>
298	>	if ( (c_perez = calloc(5, sizeof(double))) == NULL )
299		{
300	<	fprintf(stderr,"Out of memory error in function main !");
301	<	exit(1);
300	>	fprintf(stderr,"Out of memory error in function main");
301	>	return 1;
302		}
303
267	–
304		printhead(argc, argv);
269	–
305		computesky();
306	+
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		printsky();
315	<
273	<	exit(0);
315	>	return 0;
316		}
317
318
319	<	void
320	<	computesky()                    /* compute sky parameters */
319	>
320	>
321	>
322	>
323	>
324	>
325	>	void computesky()
326		{
327
328	<	/* new variables */
329	<	int     j, i;
330	<	float   *lv_mod;  /* 145 luminance values*/
331	<	/* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */
285	<	float   *theta_o, *phi_o, *coeff_perez;
328	>	int     j;
329	>
330	>	float   *lv_mod;  /* 145 luminance values */
331	>	float   *theta_o, *phi_o;
332		double  dzeta, gamma;
287	–	double  diffusion;
333		double  normfactor;
334	+	double  erbs_s0, erbs_kt;
335
336
291	–
337		/* compute solar direction */
338	<
338	>
339		if (month) {                    /* from date and time */
340		int  jd;
341		double  sd, st;
#	Line 301 | Line 346 | computesky()                   /* compute sky parameters */
346		st = hour;
347		else
348		st = hour + stadj(jd);
349	+
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		altitude = salt(sd, st);
378		azimuth = sazi(sd, st);
379
380		daynumber = (double)jdate(month, day);
381	+
382	+	}
383	+
384	+
385
386	+
387	+	/* if loop for the -l option. W.Sprenger (01/2013) */
388	+	/*
389	+	if (altitude*180/M_PI < sunaltitude_border) {
390	+
391	+	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	+	print_error_sky();
395	+	exit(1);
396		}
397	+	*/
398	+
399	+
400		if (!cloudy && altitude > 87.*M_PI/180.) {
401	<	fprintf(stderr,
401	>
402	>	if (suppress_warnings==0) {
403	>	fprintf(stderr,
404		"%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n",
405		progname);
406	<	printf(
315	<	"# warning - sun too close to zenith, reducing altitude to 87 degrees\n");
406	>	}
407		altitude = 87.*M_PI/180.;
408		}
409	+
410		sundir[0] = -sin(azimuth)*cos(altitude);
411		sundir[1] = -cos(azimuth)*cos(altitude);
412		sundir[2] = sin(altitude);
#	Line 325 | Line 417 | computesky()                   /* compute sky parameters */
417
418
419
420	<	/* compute the inputs for the calculation of the light distribution over the sky*/
421	<	if (input==0)
420	>	/* compute the inputs for the calculation of the light distribution over the sky*/
421	>	if (input==0)           /* P */
422		{
423		check_parametrization();
424	<	diffusirradiance = diffus_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
424	>	diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
425		directirradiance = direct_irradiance_from_sky_clearness();
426		check_irradiances();
427
428		if (output==0 || output==2)
429		{
430	<	diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
430	>	diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
431		directilluminance = directirradiance*direct_n_effi_PEREZ();
432		check_illuminances();
433		}
434		}
435
436
437	<	else if (input==1)
437	>	else if (input==1)      /* W */
438		{
439		check_irradiances();
440		skybrightness = sky_brightness();
#	Line 351 | Line 443 | computesky()                   /* compute sky parameters */
443
444		if (output==0 || output==2)
445		{
446	<	diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
446	>	diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
447		directilluminance = directirradiance*direct_n_effi_PEREZ();
448		check_illuminances();
449		}
#	Line 359 | Line 451 | computesky()                   /* compute sky parameters */
451		}
452
453
454	<	else if (input==2)
454	>	else if (input==2)      /* L */
455		{
456		check_illuminances();
457		illu_to_irra_index();
#	Line 367 | Line 459 | computesky()                   /* compute sky parameters */
459		}
460
461
462	<	else if (input==3)
462	>	else if (input==3)      /* G */
463		{
464		if (altitude<=0)
465		{
466	<	fprintf(stderr, "solar zenith angle larger than 90� \n the models used are not more valid\n");
467	<	exit(1);
466	>	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	>
472	>	directirradiance=directirradiance/sin(altitude);
473		}
474	<
378	<	directirradiance=directirradiance/sin(altitude);
474	>
475		check_irradiances();
476		skybrightness = sky_brightness();
477		skyclearness =  sky_clearness();
#	Line 383 | Line 479 | computesky()                   /* compute sky parameters */
479
480		if (output==0 || output==2)
481		{
482	<	diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
482	>	diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
483		directilluminance = directirradiance*direct_n_effi_PEREZ();
484		check_illuminances();
485		}
486
487		}
488
393	–
394	–	else    {fprintf(stderr,"error in giving the input arguments"); exit(1);}
489
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
543	<	/* normalization factor for the relative sky luminance distribution, diffuse part*/
543	>	else    {fprintf(stderr,"error at the input arguments"); exit(1);}
544
400	–	/* allocation dynamique de memoire pour les pointeurs */
401	–	if ( (coeff_perez = malloc(8*20*sizeof(float))) == NULL )
402	–	{
403	–	fprintf(stderr,"Out of memory error in function main !");
404	–	exit(1);
405	–	}
545
546	<	/* read the coefficients for the Perez sky luminance model */
547	<	if (lect_coeff_perez(DATFILE, &coeff_perez) > 0)
548	<	{
410	<	fprintf(stderr,"lect_coeff_perez does not work\n");
411	<	exit(2);
412	<	}
413	<
546	>
547	>	/* normalization factor for the relative sky luminance distribution, diffuse part*/
548	>
549		if ( (lv_mod = malloc(145*sizeof(float))) == NULL)
550		{
551		fprintf(stderr,"Out of memory in function main");
#	Line 418 | Line 553 | computesky()                   /* compute sky parameters */
553		}
554
555		/* read the angles */
556	<	theta_o = theta_ordered("defangle.dat");
557	<	phi_o = phi_ordered("defangle.dat");
556	>	theta_o = defangle_theta;
557	>	phi_o = defangle_phi;
558	>
559
560	<	/* parameters for the perez model */
560	>	/* parameters for the perez model */
561		coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez);
562
563	<	/*calculation of the modelled luminance */
563	>
564	>
565	>
566	>
567	>	/*calculation of the modelled luminance */
568		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	+
572		*(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez);
573	<	/*printf("theta, phi, lv_mod %lf\t %lf\t %lf\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j));*/
573	>
574	>	/* fprintf(stderr,"theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); */
575		}
576	<
576	>
577		/* integration of luminance for the normalization factor, diffuse part of the sky*/
578	+
579		diffnormalization = integ_lv(lv_mod, theta_o);
437	–	/*printf("perez integration %lf\n", diffnormalization);*/
580
439	–
581
582
583	<	/*normalization coefficient in lumen or in watt*/
583	>	/*normalization coefficient in lumen or in watt*/
584		if (output==0)
585		{
586	<	diffnormalization = diffusilluminance/diffnormalization/WHTEFFICACY;
586	>	diffnormalization = diffuseilluminance/diffnormalization/WHTEFFICACY;
587		}
588		else if (output==1)
589		{
590	<	diffnormalization = diffusirradiance/diffnormalization;
590	>	diffnormalization = diffuseirradiance/diffnormalization;
591		}
592		else if (output==2)
593		{
594	<	diffnormalization = diffusilluminance/diffnormalization;
594	>	diffnormalization = diffuseilluminance/diffnormalization;
595		}
596
597	<	else    {fprintf(stderr,"output argument : wrong number"); exit(1);}
597	>	else    {fprintf(stderr,"Wrong output specification.\n"); exit(1);}
598
599
600
601
602	<	/* calculation for the solar source */
602	>	/* calculation for the solar source */
603		if (output==0)
604		solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY;
605
#	Line 474 | Line 615 | computesky()                   /* compute sky parameters */
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	<	/*
478	<	fprintf(stderr, "gendaylit : the actual zenith radiance(W/m^2/sr) or luminance(cd/m^2) is : %.0lf\n", zenithbr);
479	<	*/
480	<
618	>
619		if (skyclearness==1)
620		normfactor = 0.777778;
621
#	Line 495 | Line 633 | if ( (skyclearness>1) && (skyclearness<6) )
633		}
634
635		groundbr = zenithbr*normfactor;
498	–	printf("# Ground ambient level: %.1f\n", groundbr);
636
637		if (dosun&&(skyclearness>1))
638	<	groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];
638	>	groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];
639
640		groundbr *= gprefl;
641
#	Line 510 | Line 647 | return;
647
648
649
650	+	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
663
664
665	<	printsky()                      /* print out sky */
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	+	void printsky()                 /* print out sky */
696	+	{
697		if (dosun&&(skyclearness>1))
698	<	{
698	>	{
699		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	<	}
527	<
528	<	if (dosun&&(skyclearness==1))
529	<	{
705	>	} else if (dosun) {
706		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	<	}
712	>	}
713
714
715		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	<	*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
720	<	sundir[0], sundir[1], sundir[2]);
719	>	*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
720	>	sundir[0], sundir[1], sundir[2]);
721	>
722		}
723
724
725	<	printdefaults()                 /* print default values */
725	>	void printdefaults()                    /* print default values */
726		{
727		printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
728		if (zenithbr > 0.0)
#	Line 558 | Line 735 | printdefaults()                        /* print default values */
735		}
736
737
738	<	userror(msg)                    /* print usage error and quit */
562	<	char  *msg;
738	>	void userror(char* msg)                 /* print usage error and quit */
739		{
740		if (msg != NULL)
741	<	fprintf(stderr, "%s: Use error - %s\n", progname, msg);
742	<	fprintf(stderr, "Usage: %s month day hour [-P|-W|-L] direct_value diffus_value [options]\n", progname);
743	<	fprintf(stderr, "or   : %s -ang altitude azimuth [-P|-W|-L] direct_value diffus_value [options]\n", progname);
741	>	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		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	+	fprintf(stderr, "       -E global-horizontal-irradiance (W/m^2)\n\n");
753	+	fprintf(stderr, "       Output specification with option:\n");
754		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	+	fprintf(stderr, "       gendaylit version 2.3 (2013/08/08)  \n\n");
756		exit(1);
757		}
758
759
760
761	<	double
579	<	normsc()                        /* compute normalization factor (E0*F2/L0) */
761	>	double normsc()           /* compute normalization factor (E0*F2/L0) */
762		{
763		static double  nfc[2][5] = {
764		/* clear sky approx. */
#	Line 599 | Line 781 | normsc()                       /* compute normalization factor (E0*F2/L0)
781
782
783
784	<	printhead(ac, av)               /* print command header */
603	<	register int  ac;
604	<	register char  **av;
784	>	void printhead(int ac, char** av)               /* print command header */
785		{
786		putchar('#');
787		while (ac--) {
#	Line 614 | Line 794 | register char  **av;
794
795
796
617	–	void
618	–	skip_comments(FILE *fp)         /* skip comments in file */
619	–	{
620	–	int     c;
621	–
622	–	while ((c = getc(fp)) != EOF)
623	–	if (c == '#') {
624	–	while ((c = getc(fp)) != EOF)
625	–	if (c == '\n')
626	–	break;
627	–	} else if (!isspace(c)) {
628	–	ungetc(c, fp);
629	–	break;
630	–	}
631	–	}
632	–
633	–
634	–
797		/* Perez models */
798
799		/* Perez global horizontal luminous efficacy model */
#	Line 641 | Line 803 | double glob_h_effi_PEREZ()
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	<
807	<
808	<	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
809	<	fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
810	<
806	>
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		/* initialize category bounds (clearness index bounds) */
815
816		category_total_number = 8;
#	Line 721 | Line 886 | double glob_h_diffuse_effi_PEREZ()
886		double  category_bounds[10], a[10], b[10], c[10], d[10];
887		int     category_total_number, category_number, i;
888
889	+
890	+
891
892	<	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
893	<	fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
892	>	check_parametrization();
893	>
894	>
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
900		/* initialize category bounds (clearness index bounds) */
901
902		category_total_number = 8;
903
904	+	//XXX:  category_bounds > 0.1
905		category_bounds[1] = 1;
906		category_bounds[2] = 1.065;
907		category_bounds[3] = 1.230;
#	Line 780 | Line 953 | fprintf(stderr, "Warning : skyclearness or skybrightne
953
954
955
956	+	category_number = -1;
957		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	+	if (category_number == -1) {
964	+	if (suppress_warnings==0)
965	+	fprintf(stderr, "ERROR: Model parameters out of range, skyclearness = %lf \n", skyclearness);
966	+	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		d[category_number]*log(skybrightness);
973
974		return(value);
975	+
976		}
977
978
979	+
980		/* direct normal efficacy model, according to PEREZ */
981
982		double direct_n_effi_PEREZ()
#	Line 803 | Line 987 | double         category_bounds[10], a[10], b[10], c[10], d[10
987		int     category_total_number, category_number, i;
988
989
990	<	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
991	<	fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
990	>	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
993
994		/* initialize category bounds (clearness index bounds) */
#	Line 878 | Line 1062 | return(value);
1062		/*check the range of epsilon and delta indexes of the perez parametrization*/
1063		void check_parametrization()
1064		{
1065	<	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
1065	>	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup)
1066		{
1067	<	fprintf(stderr,"sky clearness or sky brightness out of range %lf\t %lf\n", skyclearness, skybrightness);
1068	<	exit(1);
1067	>
1068	>	/*  limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */
1069	>	if (skyclearness<skyclearinf){
1070	>	if (suppress_warnings==0)
1071	>	/* fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */
1072	>	skyclearness=skyclearinf;
1073		}
1074	+	if (skyclearness>skyclearsup){
1075	+	if (suppress_warnings==0)
1076	+	/* fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */
1077	+	skyclearness=skyclearsup-0.1;
1078	+	}
1079	+	if (skybrightness<skybriginf){
1080	+	if (suppress_warnings==0)
1081	+	/* fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */
1082	+	skybrightness=skybriginf;
1083	+	}
1084	+	if (skybrightness>skybrigsup){
1085	+	if (suppress_warnings==0)
1086	+	/* fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */
1087	+	skybrightness=skybrigsup;
1088	+	}
1089	+
1090	+	return; }
1091		else return;
1092		}
1093
1094
1095	<	/* likelihood of the direct and diffuse components */
1095	>	/* validity of the direct and diffuse components */
1096		void    check_illuminances()
1097		{
1098	<	if (!( (directilluminance>=0) && (directilluminance<=solar_constant_l*1000) && (diffusilluminance>0) ))
1099	<	{
1100	<	fprintf(stderr,"direct or diffuse illuminances out of range\n");
896	<	exit(1);
1098	>	if (directilluminance < 0) {
1099	>	fprintf(stderr,"WARNING: direct illuminance < 0. Using 0.0\n");
1100	>	directilluminance = 0.0;
1101		}
1102	<	return;
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	>	}
1110		}
1111
1112
1113		void    check_irradiances()
1114		{
1115	<	if (!( (directirradiance>=0) && (directirradiance<=solar_constant_e) && (diffusirradiance>0) ))
1116	<	{
1117	<	fprintf(stderr,"direct or diffuse irradiances out of range\n");
1118	<	exit(1);
1119	<	}
1120	<	return;
1115	>	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		}
1128
1129
#	Line 916 | Line 1133 | double sky_brightness()
1133		{
1134		double value;
1135
1136	<	value = diffusirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
1136	>	value = diffuseirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
1137
1138		return(value);
1139		}
#	Line 925 | Line 1142 | return(value);
1142		/* Perez sky's clearness */
1143		double sky_clearness()
1144		{
1145	<	double value;
1145	>	double value;
1146
1147	<	value = ( (diffusirradiance + directirradiance)/(diffusirradiance) + 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) ;
1147	>	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
1149	<	return(value);
1149	>	return(value);
1150		}
1151
1152
1153
1154		/* diffus horizontal irradiance from Perez sky's brightness */
1155	<	double diffus_irradiance_from_sky_brightness()
1155	>	double diffuse_irradiance_from_sky_brightness()
1156		{
1157		double value;
1158
#	Line 950 | Line 1167 | double direct_irradiance_from_sky_clearness()
1167		{
1168		double value;
1169
1170	<	value = diffus_irradiance_from_sky_brightness();
1170	>	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
1173		return(value);
1174		}
1175
1176
1177	<	void illu_to_irra_index(void)
1177	>
1178	>
1179	>	void illu_to_irra_index()
1180		{
1181	<	double  test1=0.1, test2=0.1;
1181	>	double  test1=0.1, test2=0.1, d_eff;
1182		int     counter=0;
1183
1184	<	diffusirradiance = diffusilluminance*solar_constant_e/(solar_constant_l*1000);
1184	>	diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l*1000);
1185		directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000);
1186		skyclearness =  sky_clearness();
1187		skybrightness = sky_brightness();
1188	<	if (skyclearness>12) skyclearness=12;
970	<	if (skybrightness<0.05) skybrightness=0.01;
971	<
972	<
973	<	while ( ((fabs(diffusirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
974	<	|| skyclearness>skyclearinf || skyclearness<skyclearsup
975	<	|| skybrightness>skybriginf || skybrightness<skybrigsup )
976	<	&& !(counter==5) )
977	<	{
978	<	/*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
1188	>	check_parametrization();
1189
1190	<	test1=diffusirradiance;
1190	>
1191	>	while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
1192	>	|| (!(skyclearness<skyclearinf || skyclearness>skyclearsup))
1193	>	|| (!(skybrightness<skybriginf || skybrightness>skybrigsup)) )
1194	>	&& !(counter==9) )
1195	>	{
1196	>
1197	>	test1=diffuseirradiance;
1198		test2=directirradiance;
1199		counter++;
1200
1201	<	diffusirradiance = diffusilluminance/glob_h_diffuse_effi_PEREZ();
1202	<	directirradiance = directilluminance/direct_n_effi_PEREZ();
986	<	/*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
1201	>	diffuseirradiance = diffuseilluminance/glob_h_diffuse_effi_PEREZ();
1202	>	d_eff = direct_n_effi_PEREZ();
1203
1204	+
1205	+	if (d_eff < 0.1)
1206	+	directirradiance = 0;
1207	+	else
1208	+	directirradiance = directilluminance/d_eff;
1209	+
1210		skybrightness = sky_brightness();
1211		skyclearness =  sky_clearness();
1212	<	if (skyclearness>12) skyclearness=12;
1213	<	if (skybrightness<0.05) skybrightness=0.01;
1214	<
1215	<	/*fprintf(stderr, "%lf\t %lf\n", skybrightness, skyclearness);*/
994	<
1212	>	check_parametrization();
1213	>
1214	>	/*fprintf(stderr,"skyclearness = %lf, skybrightness = %lf, directirradiance = %lf, diffuseirradiance = %lf\n",skyclearness, skybrightness, directirradiance, diffuseirradiance);*/
1215	>
1216		}
1217
1218
1219		return;
1220		}
1221
1222	<
1002	<	int lect_coeff_perez(char *filename,float **coeff_perez)
1222	>	static int get_numlin(float epsilon)
1223		{
1224	<	FILE *fcoeff_perez;
1225	<	float temp;
1226	<	int i,j;
1227	<
1228	<	if ((fcoeff_perez = frlibopen(filename)) == NULL)
1229	<	{
1230	<	fprintf(stderr,"file %s cannot be opened\n", filename);
1231	<	return 1; /* il y a un probleme de fichier */
1232	<	}
1233	<	else
1234	<	{
1235	<	/*printf("file %s  open\n", filename);*/
1236	<	}
1237	<
1238	<	skip_comments(fcoeff_perez);
1019	<
1020	<	for (i=0;i<8;i++)
1021	<	for (j=0;j<20;j++)
1022	<	{
1023	<	fscanf(fcoeff_perez,"%f",&temp);
1024	<	*(*coeff_perez+i*20+j) = temp;
1025	<	}
1026	<	fclose(fcoeff_perez);
1027	<
1028	<	return 0; /* tout est OK */
1224	>	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		}
1240
1031	–
1032	–
1241		/* sky luminance perez model */
1242	<	double calc_rel_lum_perez(double dzeta,double gamma,double Z,
1035	<	double epsilon,double Delta,float *coeff_perez)
1242	>	double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[])
1243		{
1244	+
1245		float x[5][4];
1246		int i,j,num_lin;
1247		double c_perez[5];
1248
1249		if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
1250		{
1251	<	fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n");
1251	>	fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez!\n");
1252		exit(1);
1253		}
1254
#	Line 1049 | Line 1257 | double epsilon,double Delta,float *coeff_perez)
1257		{
1258		if ( Delta < 0.2 ) Delta = 0.2;
1259		}
1260	<
1261	<	if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0;
1262	<	if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
1263	<	if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2;
1056	<	if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3;
1057	<	if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4;
1058	<	if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5;
1059	<	if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6;
1060	<	if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7;
1061	<
1260	>
1261	>
1262	>	num_lin = get_numlin(epsilon);
1263	>
1264		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	<	/* printf("x %d %d vaut %f\n",i,j,x[i][j]); */
1268	>	/* fprintf(stderr,"x %d %d vaut %f\n",i,j,x[i][j]); */
1269		}
1270
1271
#	Line 1090 | Line 1292 | double epsilon,double Delta,float *coeff_perez)
1292
1293
1294		/* coefficients for the sky luminance perez model */
1295	<	void coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez)
1295	>	void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[])
1296		{
1297		float x[5][4];
1298		int i,j,num_lin;
#	Line 1106 | Line 1308 | void coeff_lum_perez(double Z, double epsilon, double
1308		{
1309		if ( Delta < 0.2 ) Delta = 0.2;
1310		}
1311	+
1312	+
1313	+	num_lin = get_numlin(epsilon);
1314
1315	<	if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0;
1111	<	if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
1112	<	if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2;
1113	<	if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3;
1114	<	if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4;
1115	<	if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5;
1116	<	if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6;
1117	<	if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7;
1315	>	/*fprintf(stderr,"numlin %d\n", num_lin);*/
1316
1317		for (i=0;i<5;i++)
1318		for (j=0;j<4;j++)
#	Line 1146 | Line 1344 | void coeff_lum_perez(double Z, double epsilon, double
1344		}
1345
1346
1347	+
1348		/* degrees into radians */
1349		double radians(double degres)
1350		{
1351		return degres*M_PI/180.0;
1352		}
1353
1354	+
1355		/* radian into degrees */
1356		double degres(double radians)
1357		{
1358		return radians/M_PI*180.0;
1359		}
1360
1361	+
1362		/* 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		{
#	Line 1174 | Line 1375 | void theta_phi_to_dzeta_gamma(double theta,double phi,
1375		}
1376
1377
1177	–	/********************************************************************************/
1178	–	/*      Fonction: theta_ordered                                                 */
1179	–	/*                                                                              */
1180	–	/*      In: char *filename                                                      */
1181	–	/*                                                                              */
1182	–	/*      Out: float *                                                            */
1183	–	/*                                                                              */
1184	–	/*      Update: 29/08/93                                                        */
1185	–	/*                                                                              */
1186	–	/*      Rem: theta en degres                                                    */
1187	–	/*                                                                              */
1188	–	/*      But: fournit les valeurs de theta du fichier d'entree a la memoire      */
1189	–	/*                                                                              */
1190	–	/********************************************************************************/
1191	–	float *theta_ordered(char *filename)
1192	–	{
1193	–	int i;
1194	–	float buffer,*ptr;
1195	–	FILE *file_in;
1378
1197	–	if ( (file_in = frlibopen(filename)) == NULL )
1198	–	{
1199	–	fprintf(stderr,"Cannot open file %s in function theta_ordered\n",filename);
1200	–	exit(1);
1201	–	}
1202	–
1203	–	skip_comments(file_in);
1204	–
1205	–	if ( (ptr = malloc(145*sizeof(float))) == NULL )
1206	–	{
1207	–	fprintf(stderr,"Out of memory in function theta_ordered\n");
1208	–	exit(1);
1209	–	}
1210	–
1211	–	for (i=0;i<145;i++)
1212	–	{
1213	–	fscanf(file_in,"%f",&buffer);
1214	–	*(ptr+i) = buffer;
1215	–	fscanf(file_in,"%f",&buffer);
1216	–	}
1217	–
1218	–	fclose(file_in);
1219	–	return ptr;
1220	–	}
1221	–
1222	–
1223	–	/********************************************************************************/
1224	–	/*      Fonction: phi_ordered                                                   */
1225	–	/*                                                                              */
1226	–	/*      In: char *filename                                                      */
1227	–	/*                                                                              */
1228	–	/*      Out: float *                                                            */
1229	–	/*                                                                              */
1230	–	/*      Update: 29/08/93                                                        */
1231	–	/*                                                                              */
1232	–	/*      Rem: valeurs de Phi en DEGRES                                           */
1233	–	/*                                                                              */
1234	–	/*      But: mettre les angles contenus dans le fichier d'entree dans la memoire */
1235	–	/*                                                                              */
1236	–	/********************************************************************************/
1237	–	float *phi_ordered(char *filename)
1238	–	{
1239	–	int i;
1240	–	float buffer,*ptr;
1241	–	FILE *file_in;
1242	–
1243	–	if ( (file_in = frlibopen(filename)) == NULL )
1244	–	{
1245	–	fprintf(stderr,"Cannot open file %s in function phi_ordered\n",filename);
1246	–	exit(1);
1247	–	}
1248	–
1249	–	skip_comments(file_in);
1250	–
1251	–	if ( (ptr = malloc(145*sizeof(float))) == NULL )
1252	–	{
1253	–	fprintf(stderr,"Out of memory in function phi_ordered");
1254	–	exit(1);
1255	–	}
1256	–
1257	–	for (i=0;i<145;i++)
1258	–	{
1259	–	fscanf(file_in,"%f",&buffer);
1260	–	fscanf(file_in,"%f",&buffer);
1261	–	*(ptr+i) = buffer;
1262	–	}
1263	–
1264	–	fclose(file_in);
1265	–	return ptr;
1266	–	}
1267	–
1268	–
1269	–	/********************************************************************************/
1270	–	/*      Fonction: integ_lv                                                      */
1271	–	/*                                                                              */
1272	–	/*      In: float *lv,*theta                                                    */
1273	–	/*          int sun_pos                                                         */
1274	–	/*                                                                              */
1275	–	/*      Out: double                                                             */
1276	–	/*                                                                              */
1277	–	/*      Update: 29/08/93                                                        */
1278	–	/*                                                                              */
1279	–	/*      Rem:                                                                    */
1280	–	/*                                                                              */
1281	–	/*      But: calcul l'integrale de luminance relative sans la dir. du soleil    */
1282	–	/*                                                                              */
1283	–	/********************************************************************************/
1379		double integ_lv(float *lv,float *theta)
1380		{
1381		int i;
1382		double buffer=0.0;
1383	<
1383	>
1384		for (i=0;i<145;i++)
1385	+	{
1386		buffer += (*(lv+i))*cos(radians(*(theta+i)));
1387	<
1387	>	}
1388	>
1389		return buffer*2*M_PI/144;
1293	–
1390		}
1391
1392
1393
1298	–
1299	–
1300	–
1394		/* enter day number(double), return E0 = square(R0/R):  eccentricity correction factor  */
1395
1396		double get_eccentricity()
#	Line 1310 | Line 1403 | double get_eccentricity()
1403		0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle);
1404
1405		return (E0);
1313	–
1406		}
1407
1408
#	Line 1318 | Line 1410 | double get_eccentricity()
1410		double  air_mass()
1411		{
1412		double  m;
1321	–
1413		if (sunzenith>90)
1414		{
1415	<	fprintf(stderr, "solar zenith angle larger than 90� in fuction air_mass():\n the models used are not more valid\n");
1415	>	fprintf(stderr, "Solar zenith angle larger than 90 degrees in function air_mass()\n");
1416		exit(1);
1417		}
1327	–
1418		m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1419		return(m);
1420		}
1331	–
1332	–
1333	–	double get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut)
1334	–
1335	–	{
1336	–
1337	–	double angle;
1338	–
1339	–
1340	–	if (sun_zenith == 0)
1341	–	puts("WARNING: zenith_angle = 0 in function get_angle_sun_vert_plan");
1342	–
1343	–	angle = acos( cos(sun_zenith*M_PI/180)*cos(direction_zenith*M_PI/180) + sin(sun_zenith*M_PI/180)*sin(direction_zenith*M_PI/180)*cos((sun_azimut-direction_azimut)*M_PI/180) );
1344	–	angle = angle*180/M_PI;
1345	–	return(angle);
1346	–	}
1347	–
1348	–
1349	–
1350	–
1351	–
1421
1422
1423

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines