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.4 by greg, Wed Aug 10 22:28:45 2011 UTC vs.
Revision 2.13 by greg, Wed Aug 14 17:11:43 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
9	<
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	>	#define _USE_MATH_DEFINES
10		#include  <stdio.h>
11		#include  <string.h>
12		#include  <math.h>
13		#include  <stdlib.h>
24	–	#include  <ctype.h>
14
26	–	#include  "rtio.h"
27	–	#include  "fvect.h"
15		#include  "color.h"
16	+	#include  "sun.h"
17		#include  "paths.h"
18
19	<	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);
19	>	#define  DOT(v1,v2)     (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
20
21		double  normsc();
22
23	<	#define DATFILE         "coeff_perez.dat"
23	>	/*static        char *rcsid="$Header$";*/
24
25	+	float coeff_perez[] = {
26	+	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	+
57	+
58	+
59		/* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */
60		double sky_brightness();
61		double sky_clearness();
62
63		/* calculation of the direct and diffuse components from the Perez parametrization */
64	<	double  diffus_irradiance_from_sky_brightness();
64	>	double  diffuse_irradiance_from_sky_brightness();
65		double  direct_irradiance_from_sky_clearness();
66
67	+	/* 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
52	–	/* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */
70		double  glob_h_effi_PEREZ();
71		double  glob_h_diffuse_effi_PEREZ();
72		double  direct_n_effi_PEREZ();
73	+
74		/*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	+	void    print_error_sky();
80
81	<
82	<	/* Perez sky luminance model */
64	<	int     lect_coeff_perez(char *filename,float **coeff_perez);
65	<	double  calc_rel_lum_perez(double dzeta,double gamma,double Z,
66	<	double epsilon,double Delta,float *coeff_perez);
67	<	/* coefficients for the sky luminance perez model */
68	<	void    coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez);
81	>	double  calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]);
82	>	void    coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]);
83		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);
73	–	float   *theta_ordered(char *filename);
74	–	float   *phi_ordered(char *filename);
75	–	void    skip_comments(FILE *fp);
87
88	+	void printdefaults();
89	+	void check_sun_position();
90	+	void computesky();
91	+	void printhead(int ac, char** av);
92	+	void userror(char* msg);
93	+	void printsky();
94
95	+	FILE * frlibopen(char* fname);
96
97		/* astronomy and geometry*/
98		double  get_eccentricity();
99		double  air_mass();
82	–	double  get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut);
100
101	<
102	<	/* date*/
101	>	double  solar_sunset();
102	>	double  solar_sunrise();
103	>	double  stadj();
104		int     jdate(int month, int day);
105
106
107
90	–
91	–
108		/* sun calculation constants */
109	<	extern double  s_latitude;
110	<	extern double  s_longitude;
111	<	extern double  s_meridian;
109	>	extern double   s_latitude;
110	>	extern double   s_longitude;
111	>	extern double   s_meridian;
112
113		const double    AU = 149597890E3;
114		const double    solar_constant_e = 1367;    /* solar constant W/m^2 */
#	Line 101 | Line 117 | const double   solar_constant_l = 127.5;   /* solar co
117		const double    half_sun_angle = 0.2665;
118		const double    half_direct_angle = 2.85;
119
120	<	const double    skyclearinf = 1.000;    /* limitations for the variation of the Perez parameters */
120	>	const double    skyclearinf = 1.0;      /* limitations for the variation of the Perez parameters */
121		const double    skyclearsup = 12.1;
122		const double    skybriginf = 0.01;
123		const double    skybrigsup = 0.6;
#	Line 117 | Line 133 | double  altitude, azimuth;                     /* or solar angles */
133
134
135		/* definition of the sky conditions through the Perez parametrization */
136	<	double  skyclearness, skybrightness;
137	<	double  solarradiance;  /*radiance of the sun disk and of the circumsolar area*/
138	<	double  diffusilluminance, directilluminance, diffusirradiance, directirradiance;
139	<	double  sunzenith, daynumber=150, atm_preci_water=2;
136	>	double  skyclearness = 0;
137	>	double  skybrightness = 0;
138	>	double  solarradiance;
139	>	double  diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance;
140	>	double  sunzenith, daynumber, atm_preci_water=2;
141
142	<	double  diffnormalization, dirnormalization;
142	>	/*double  sunaltitude_border = 0;*/
143	>	double  diffnormalization = 0;
144	>	double  dirnormalization = 0;
145		double  *c_perez;
146
147	<	int     output=0;       /*define the unit of the output (sky luminance or radiance): visible watt=0, solar watt=1, lumen=2*/
148	<	int     input=0;        /*define the input for the calulation*/
147	>	int     output=0;       /* define the unit of the output (sky luminance or radiance): */
148	>	/* visible watt=0, solar watt=1, lumen=2 */
149	>	int     input=0;        /* define the input for the calulation */
150
151	+	int     suppress_warnings=0;
152	+
153		/* default values */
154	<	int  cloudy = 0;                                /* 1=standard, 2=uniform */
155	<	int  dosun = 1;
154	>	int     cloudy = 0;                             /* 1=standard, 2=uniform */
155	>	int     dosun = 1;
156		double  zenithbr = -1.0;
157		double  betaturbidity = 0.1;
158		double  gprefl = 0.2;
#	Line 138 | Line 160 | int    S_INTER=0;
160
161		/* computed values */
162		double  sundir[3];
163	<	double  groundbr;
163	>	double  groundbr = 0;
164		double  F2;
165		double  solarbr = 0.0;
166		int     u_solar = 0;                            /* -1=irradiance, 1=radiance */
167	+	float   timeinterval = 0;
168
169	<	char  *progname;
170	<	char  errmsg[128];
169	>	char    *progname;
170	>	char    errmsg[128];
171
172
173	<	main(argc, argv)
174	<	int  argc;
175	<	char  *argv[];
173	>
174	>
175	>	int main(int argc, char** argv)
176		{
177		int  i;
178
179		progname = argv[0];
180		if (argc == 2 && !strcmp(argv[1], "-defaults")) {
181		printdefaults();
182	<	exit(0);
182	>	return 0;
183		}
184		if (argc < 4)
185		userror("arg count");
#	Line 183 | Line 206 | char  *argv[];
206		cloudy = 0;
207		dosun = argv[i][0] == '+';
208		break;
186	–	case 'r':
209		case 'R':
210		u_solar = argv[i][1] == 'R' ? -1 : 1;
211		solarbr = atof(argv[++i]);
#	Line 195 | Line 217 | char  *argv[];
217		case 't':
218		betaturbidity = atof(argv[++i]);
219		break;
220	+	case 'w':
221	+	suppress_warnings = 1;
222	+	break;
223		case 'b':
224		zenithbr = atof(argv[++i]);
225		break;
#	Line 210 | Line 235 | char  *argv[];
235		case 'm':
236		s_meridian = atof(argv[++i]) * (M_PI/180);
237		break;
213	–
238
239		case 'O':
240		output = atof(argv[++i]);       /*define the unit of the output of the program :
241	<	sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm)
218	<	default is set to 0*/
241	>	sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm) */
242		break;
243
244		case 'P':
#	Line 227 | Line 250 | char  *argv[];
250		case 'W':                                       /* direct normal Irradiance [W/m^2] */
251		input = 1;                              /* diffuse horizontal Irrad. [W/m^2] */
252		directirradiance = atof(argv[++i]);
253	<	diffusirradiance = atof(argv[++i]);
253	>	diffuseirradiance = atof(argv[++i]);
254		break;
255
256		case 'L':                                       /* direct normal Illuminance [Lux] */
257		input = 2;                              /* diffuse horizontal Ill. [Lux] */
258		directilluminance = atof(argv[++i]);
259	<	diffusilluminance = atof(argv[++i]);
259	>	diffuseilluminance = atof(argv[++i]);
260		break;
261
262		case 'G':                                       /* direct horizontal Irradiance [W/m^2] */
263		input = 3;                              /* diffuse horizontal Irrad. [W/m^2] */
264		directirradiance = atof(argv[++i]);
265	<	diffusirradiance = atof(argv[++i]);
265	>	diffuseirradiance = atof(argv[++i]);
266		break;
267	+
268	+	case 'E':                                       /* Erbs model based on the */
269	+	input = 4;                              /* global-horizontal irradiance [W/m^2] */
270	+	globalirradiance = atof(argv[++i]);
271	+	break;
272	+
273	+	/*
274	+	case 'l':
275	+	sunaltitude_border = atof(argv[++i]);
276	+	break;
277	+	*/
278
279	+	case 'i':
280	+	timeinterval = atof(argv[++i]);
281	+	break;
282
283	+
284		default:
285		sprintf(errmsg, "unknown option: %s", argv[i]);
286		userror(errmsg);
#	Line 256 | Line 294 | char  *argv[];
294		progname, (s_longitude-s_meridian)*12/M_PI);
295
296
297	<	/* allocation dynamique de memoire pour les pointeurs */
298	<	if ( (c_perez = malloc(5*sizeof(double))) == NULL )
297	>	/* dynamic memory allocation for the pointers */
298	>
299	>	if ( (c_perez = calloc(5, sizeof(double))) == NULL )
300		{
301	<	fprintf(stderr,"Out of memory error in function main !");
302	<	exit(1);
301	>	fprintf(stderr,"Out of memory error in function main");
302	>	return 1;
303		}
304
266	–
305		printhead(argc, argv);
268	–
306		computesky();
307	+
308	+	if(*(c_perez+1)>0)
309	+	{
310	+	fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1));
311	+	print_error_sky();
312	+	exit(1);
313	+	}
314	+
315		printsky();
316	<
272	<	exit(0);
316	>	return 0;
317		}
318
319
320	<	void
321	<	computesky()                    /* compute sky parameters */
320	>
321	>
322	>
323	>
324	>
325	>
326	>	void computesky()
327		{
328
329	<	/* new variables */
330	<	int     j, i;
331	<	float   *lv_mod;  /* 145 luminance values*/
332	<	/* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */
284	<	float   *theta_o, *phi_o, *coeff_perez;
329	>	int     j;
330	>
331	>	float   *lv_mod;  /* 145 luminance values */
332	>	float   *theta_o, *phi_o;
333		double  dzeta, gamma;
286	–	double  diffusion;
334		double  normfactor;
335	+	double  erbs_s0, erbs_kt;
336
337
290	–
338		/* compute solar direction */
339	<
339	>
340		if (month) {                    /* from date and time */
341		int  jd;
342		double  sd, st;
#	Line 300 | Line 347 | computesky()                   /* compute sky parameters */
347		st = hour;
348		else
349		st = hour + stadj(jd);
350	+
351	+
352	+	if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) {
353	+	print_error_sky();
354	+	exit(1);
355	+	}
356	+
357	+
358	+	if(timeinterval) {
359	+
360	+	if(timeinterval<0) {
361	+	fprintf(stderr, "time interval negative\n");
362	+	exit(1);
363	+	}
364	+
365	+	if(fabs(solar_sunrise(month,day)-st)<timeinterval/60) {
366	+
367	+	fprintf(stderr, "Solar position corrected at %d %d %.3f\n",month,day,hour);
368	+	st= (st+timeinterval/120+solar_sunrise(month,day))/2;
369	+	}
370	+
371	+	if(fabs(solar_sunset(month,day)-st)<timeinterval/60) {
372	+	fprintf(stderr, "Solar position corrected at %d %d %.3f\n",month,day,hour);
373	+	st= (st-timeinterval/120+solar_sunset(month,day))/2;
374	+	}
375	+	}
376	+
377	+
378		altitude = salt(sd, st);
379		azimuth = sazi(sd, st);
380
381		daynumber = (double)jdate(month, day);
382	+
383	+	}
384	+
385	+
386
387	+
388	+	/* if loop for the -l option. W.Sprenger (01/2013) */
389	+	/*
390	+	if (altitude*180/M_PI < sunaltitude_border) {
391	+
392	+	if (suppress_warnings==0) {
393	+	fprintf(stderr, "Warning: sun altitude (%.3f degrees) below the border (%.3f degrees)\n",altitude*180/M_PI,sunaltitude_border);
394	+	}
395	+	print_error_sky();
396	+	exit(1);
397		}
398	+	*/
399	+
400	+
401		if (!cloudy && altitude > 87.*M_PI/180.) {
402	<	fprintf(stderr,
402	>
403	>	if (suppress_warnings==0) {
404	>	fprintf(stderr,
405		"%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n",
406		progname);
407	<	printf(
314	<	"# warning - sun too close to zenith, reducing altitude to 87 degrees\n");
407	>	}
408		altitude = 87.*M_PI/180.;
409		}
410	+
411		sundir[0] = -sin(azimuth)*cos(altitude);
412		sundir[1] = -cos(azimuth)*cos(altitude);
413		sundir[2] = sin(altitude);
#	Line 324 | Line 418 | computesky()                   /* compute sky parameters */
418
419
420
421	<	/* compute the inputs for the calculation of the light distribution over the sky*/
422	<	if (input==0)
421	>	/* compute the inputs for the calculation of the light distribution over the sky*/
422	>	if (input==0)           /* P */
423		{
424		check_parametrization();
425	<	diffusirradiance = diffus_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
425	>	diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
426		directirradiance = direct_irradiance_from_sky_clearness();
427		check_irradiances();
428
429		if (output==0 || output==2)
430		{
431	<	diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
431	>	diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
432		directilluminance = directirradiance*direct_n_effi_PEREZ();
433		check_illuminances();
434		}
435		}
436
437
438	<	else if (input==1)
438	>	else if (input==1)      /* W */
439		{
440		check_irradiances();
441		skybrightness = sky_brightness();
#	Line 350 | Line 444 | computesky()                   /* compute sky parameters */
444
445		if (output==0 || output==2)
446		{
447	<	diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
447	>	diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
448		directilluminance = directirradiance*direct_n_effi_PEREZ();
449		check_illuminances();
450		}
#	Line 358 | Line 452 | computesky()                   /* compute sky parameters */
452		}
453
454
455	<	else if (input==2)
455	>	else if (input==2)      /* L */
456		{
457		check_illuminances();
458		illu_to_irra_index();
#	Line 366 | Line 460 | computesky()                   /* compute sky parameters */
460		}
461
462
463	<	else if (input==3)
463	>	else if (input==3)      /* G */
464		{
465		if (altitude<=0)
466		{
467	<	fprintf(stderr, "solar zenith angle larger than 90� \n the models used are not more valid\n");
468	<	exit(1);
467	>	if (suppress_warnings==0)
468	>	fprintf(stderr, "Warning: solar zenith angle larger than 90 degrees; using zero irradiance to proceed\n");
469	>	directirradiance = 0;
470	>	diffuseirradiance = 0;
471	>	} else {
472	>
473	>	directirradiance=directirradiance/sin(altitude);
474		}
475	<
377	<	directirradiance=directirradiance/sin(altitude);
475	>
476		check_irradiances();
477		skybrightness = sky_brightness();
478		skyclearness =  sky_clearness();
#	Line 382 | Line 480 | computesky()                   /* compute sky parameters */
480
481		if (output==0 || output==2)
482		{
483	<	diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
483	>	diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
484		directilluminance = directirradiance*direct_n_effi_PEREZ();
485		check_illuminances();
486		}
487
488		}
489
392	–
393	–	else    {fprintf(stderr,"error in giving the input arguments"); exit(1);}
490
491	+	else if (input==4)      /* E */         /* Implementation of the Erbs model. W.Sprenger (04/13) */
492	+	{
493	+
494	+	if (altitude<=0)
495	+	{
496	+	if (suppress_warnings==0 && globalirradiance > 50)
497	+	fprintf(stderr, "Warning: global irradiance higher than 50 W/m^2 while the sun altitude is lower than zero\n");
498	+	globalirradiance = 0; diffuseirradiance = 0; directirradiance = 0;
499	+
500	+	} else {
501	+
502	+	erbs_s0 = solar_constant_e*get_eccentricity()*sin(altitude);
503	+
504	+	if (globalirradiance>erbs_s0)
505	+	{
506	+	if (suppress_warnings==0)
507	+	fprintf(stderr, "Warning: global irradiance is higher than the time-dependent solar constant s0\n");
508	+
509	+	globalirradiance=erbs_s0*0.999;
510	+	}
511	+
512	+	erbs_kt=globalirradiance/erbs_s0;
513	+
514	+	if (erbs_kt<=0.22)      diffuseirradiance=globalirradiance*(1-0.09*erbs_kt);
515	+	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));
516	+	else if (erbs_kt<1)     diffuseirradiance=globalirradiance*(0.165);
517	+
518	+	directirradiance=globalirradiance-diffuseirradiance;
519	+
520	+	printf("# erbs_s0, erbs_kt, irr_dir_h, irr_diff: %.3f %.3f %.3f %.3f\n", erbs_s0, erbs_kt, directirradiance, diffuseirradiance);
521	+	printf("# WARNING: the -E option is only recommended for a rough estimation!");
522	+
523	+	directirradiance=directirradiance/sin(altitude);
524	+
525	+	}
526	+
527	+	check_irradiances();
528	+	skybrightness = sky_brightness();
529	+	skyclearness =  sky_clearness();
530	+	check_parametrization();
531
532	+	if (output==0 || output==2)
533	+	{
534	+	diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
535	+	directilluminance = directirradiance*direct_n_effi_PEREZ();
536	+	check_illuminances();
537	+	}
538	+
539	+	}
540	+
541	+
542	+
543
544	<	/* normalization factor for the relative sky luminance distribution, diffuse part*/
544	>	else    {fprintf(stderr,"error at the input arguments"); exit(1);}
545
399	–	/* allocation dynamique de memoire pour les pointeurs */
400	–	if ( (coeff_perez = malloc(8*20*sizeof(float))) == NULL )
401	–	{
402	–	fprintf(stderr,"Out of memory error in function main !");
403	–	exit(1);
404	–	}
546
547	<	/* read the coefficients for the Perez sky luminance model */
548	<	if (lect_coeff_perez(DATFILE, &coeff_perez) > 0)
549	<	{
409	<	fprintf(stderr,"lect_coeff_perez does not work\n");
410	<	exit(2);
411	<	}
412	<
547	>
548	>	/* normalization factor for the relative sky luminance distribution, diffuse part*/
549	>
550		if ( (lv_mod = malloc(145*sizeof(float))) == NULL)
551		{
552		fprintf(stderr,"Out of memory in function main");
#	Line 417 | Line 554 | computesky()                   /* compute sky parameters */
554		}
555
556		/* read the angles */
557	<	theta_o = theta_ordered("defangle.dat");
558	<	phi_o = phi_ordered("defangle.dat");
557	>	theta_o = defangle_theta;
558	>	phi_o = defangle_phi;
559	>
560
561	<	/* parameters for the perez model */
561	>	/* parameters for the perez model */
562		coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez);
563
564	<	/*calculation of the modelled luminance */
564	>
565	>
566	>
567	>
568	>	/*calculation of the modelled luminance */
569		for (j=0;j<145;j++)
570		{
571		theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith));
572	+
573		*(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez);
574	<	/*printf("theta, phi, lv_mod %lf\t %lf\t %lf\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j));*/
574	>
575	>	/* fprintf(stderr,"theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); */
576		}
577	<
577	>
578		/* integration of luminance for the normalization factor, diffuse part of the sky*/
579	+
580		diffnormalization = integ_lv(lv_mod, theta_o);
436	–	/*printf("perez integration %lf\n", diffnormalization);*/
581
438	–
582
583
584	<	/*normalization coefficient in lumen or in watt*/
584	>	/*normalization coefficient in lumen or in watt*/
585		if (output==0)
586		{
587	<	diffnormalization = diffusilluminance/diffnormalization/WHTEFFICACY;
587	>	diffnormalization = diffuseilluminance/diffnormalization/WHTEFFICACY;
588		}
589		else if (output==1)
590		{
591	<	diffnormalization = diffusirradiance/diffnormalization;
591	>	diffnormalization = diffuseirradiance/diffnormalization;
592		}
593		else if (output==2)
594		{
595	<	diffnormalization = diffusilluminance/diffnormalization;
595	>	diffnormalization = diffuseilluminance/diffnormalization;
596		}
597
598	<	else    {fprintf(stderr,"output argument : wrong number"); exit(1);}
598	>	else    {fprintf(stderr,"Wrong output specification.\n"); exit(1);}
599
600
601
602
603	<	/* calculation for the solar source */
603	>	/* calculation for the solar source */
604		if (output==0)
605		solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY;
606
#	Line 473 | Line 616 | computesky()                   /* compute sky parameters */
616		/* Compute the ground radiance */
617		zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
618		zenithbr*=diffnormalization;
619	<	/*
477	<	fprintf(stderr, "gendaylit : the actual zenith radiance(W/m^2/sr) or luminance(cd/m^2) is : %.0lf\n", zenithbr);
478	<	*/
479	<
619	>
620		if (skyclearness==1)
621		normfactor = 0.777778;
622
#	Line 494 | Line 634 | if ( (skyclearness>1) && (skyclearness<6) )
634		}
635
636		groundbr = zenithbr*normfactor;
497	–	printf("# Ground ambient level: %.1f\n", groundbr);
637
638		if (dosun&&(skyclearness>1))
639	<	groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];
639	>	groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];
640
641		groundbr *= gprefl;
642
#	Line 509 | Line 648 | return;
648
649
650
651	+	void print_error_sky()
652	+	{
653	+	sundir[0] = -sin(azimuth)*cos(altitude);
654	+	sundir[1] = -cos(azimuth)*cos(altitude);
655	+	sundir[2] = sin(altitude);
656	+
657	+	printf("\nvoid brightfunc skyfunc\n");
658	+	printf("2 skybright perezlum.cal\n");
659	+	printf("0\n");
660	+	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]);
661	+	}
662	+
663
664
665
666	<	printsky()                      /* print out sky */
666	>
667	>
668	>
669	>	double solar_sunset(int month,int day)
670		{
671	+	float W;
672	+	extern double s_latitude;
673	+	W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day))));
674	+	return(12+(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15));
675	+	}
676	+
677	+
678	+	double solar_sunrise(int month,int day)
679	+	{
680	+	float W;
681	+	extern double s_latitude;
682	+	W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day))));
683	+	return(12-(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15));
684	+	}
685	+
686	+
687	+
688	+
689	+
690	+
691	+
692	+
693	+
694	+
695	+
696	+	void printsky()                 /* print out sky */
697	+	{
698		if (dosun&&(skyclearness>1))
699	<	{
699	>	{
700		printf("\nvoid light solar\n");
701		printf("0\n0\n");
702		printf("3 %.3e %.3e %.3e\n", solarradiance, solarradiance, solarradiance);
703		printf("\nsolar source sun\n");
704		printf("0\n0\n");
705		printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
706	<	}
526	<
527	<	if (dosun&&(skyclearness==1))
528	<	{
706	>	} else if (dosun) {
707		printf("\nvoid light solar\n");
708		printf("0\n0\n");
709		printf("3 0.0 0.0 0.0\n");
710		printf("\nsolar source sun\n");
711		printf("0\n0\n");
712		printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
713	<	}
713	>	}
714
715
716		printf("\nvoid brightfunc skyfunc\n");
717		printf("2 skybright perezlum.cal\n");
718		printf("0\n");
719		printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr,
720	<	*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
721	<	sundir[0], sundir[1], sundir[2]);
720	>	*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
721	>	sundir[0], sundir[1], sundir[2]);
722	>
723		}
724
725
726	<	printdefaults()                 /* print default values */
726	>	void printdefaults()                    /* print default values */
727		{
728		printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
729		if (zenithbr > 0.0)
#	Line 557 | Line 736 | printdefaults()                        /* print default values */
736		}
737
738
739	<	userror(msg)                    /* print usage error and quit */
561	<	char  *msg;
739	>	void userror(char* msg)                 /* print usage error and quit */
740		{
741		if (msg != NULL)
742	<	fprintf(stderr, "%s: Use error - %s\n", progname, msg);
743	<	fprintf(stderr, "Usage: %s month day hour [-P|-W|-L] direct_value diffus_value [options]\n", progname);
744	<	fprintf(stderr, "or   : %s -ang altitude azimuth [-P|-W|-L] direct_value diffus_value [options]\n", progname);
742	>	fprintf(stderr, "%s: Use error - %s\n\n", progname, msg);
743	>	fprintf(stderr, "Usage: %s      month day hour    [...]\n", progname);
744	>	fprintf(stderr, "   or: %s -ang altitude azimuth  [...]\n", progname);
745	>	fprintf(stderr, "               followed by:      -P          epsilon delta [options]\n");
746	>	fprintf(stderr, "                        or:      [-W|-L|-G]  direct_value diffuse_value [options]\n");
747	>	fprintf(stderr, "                        or:      -E          global_irradiance [options]\n\n");
748	>	fprintf(stderr, "       Description:\n");
749		fprintf(stderr, "       -P epsilon delta  (these are the Perez parameters) \n");
750		fprintf(stderr, "       -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
751		fprintf(stderr, "       -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n");
752		fprintf(stderr, "       -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
753	+	fprintf(stderr, "       -E global-horizontal-irradiance (W/m^2)\n\n");
754	+	fprintf(stderr, "       Output specification with option:\n");
755		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");
756	+	fprintf(stderr, "       gendaylit version 2.3 (2013/08/08)  \n\n");
757		exit(1);
758		}
759
760
761
762	<	double
578	<	normsc()                        /* compute normalization factor (E0*F2/L0) */
762	>	double normsc()           /* compute normalization factor (E0*F2/L0) */
763		{
764		static double  nfc[2][5] = {
765		/* clear sky approx. */
#	Line 598 | Line 782 | normsc()                       /* compute normalization factor (E0*F2/L0)
782
783
784
785	<	printhead(ac, av)               /* print command header */
602	<	register int  ac;
603	<	register char  **av;
785	>	void printhead(int ac, char** av)               /* print command header */
786		{
787		putchar('#');
788		while (ac--) {
#	Line 613 | Line 795 | register char  **av;
795
796
797
616	–	void
617	–	skip_comments(FILE *fp)         /* skip comments in file */
618	–	{
619	–	int     c;
620	–
621	–	while ((c = getc(fp)) != EOF)
622	–	if (c == '#') {
623	–	while ((c = getc(fp)) != EOF)
624	–	if (c == '\n')
625	–	break;
626	–	} else if (!isspace(c)) {
627	–	ungetc(c, fp);
628	–	break;
629	–	}
630	–	}
631	–
632	–
633	–
798		/* Perez models */
799
800		/* Perez global horizontal luminous efficacy model */
#	Line 640 | Line 804 | double glob_h_effi_PEREZ()
804		double  value;
805		double  category_bounds[10], a[10], b[10], c[10], d[10];
806		int     category_total_number, category_number, i;
807	<
808	<
809	<	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
810	<	fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
811	<
807	>
808	>	check_parametrization();
809	>
810	>
811	>	/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
812	>	fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n"); */
813	>
814	>
815		/* initialize category bounds (clearness index bounds) */
816
817		category_total_number = 8;
#	Line 720 | Line 887 | double glob_h_diffuse_effi_PEREZ()
887		double  category_bounds[10], a[10], b[10], c[10], d[10];
888		int     category_total_number, category_number, i;
889
890	+
891	+
892
893	<	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
894	<	fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
893	>	check_parametrization();
894	>
895	>
896	>	/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
897	>	fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */
898	>
899	>
900
901		/* initialize category bounds (clearness index bounds) */
902
903		category_total_number = 8;
904
905	+	//XXX:  category_bounds > 0.1
906		category_bounds[1] = 1;
907		category_bounds[2] = 1.065;
908		category_bounds[3] = 1.230;
#	Line 779 | Line 954 | fprintf(stderr, "Warning : skyclearness or skybrightne
954
955
956
957	+	category_number = -1;
958		for (i=1; i<=category_total_number; i++)
959		{
960		if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
961		category_number = i;
962		}
963
964	+	if (category_number == -1) {
965	+	if (suppress_warnings==0)
966	+	fprintf(stderr, "ERROR: Model parameters out of range, skyclearness = %lf \n", skyclearness);
967	+	print_error_sky();
968	+	exit(1);
969	+	}
970	+
971	+
972		value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*cos(sunzenith*M_PI/180) +
973		d[category_number]*log(skybrightness);
974
975		return(value);
976	+
977		}
978
979
980	+
981		/* direct normal efficacy model, according to PEREZ */
982
983		double direct_n_effi_PEREZ()
#	Line 802 | Line 988 | double         category_bounds[10], a[10], b[10], c[10], d[10
988		int     category_total_number, category_number, i;
989
990
991	<	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
992	<	fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
991	>	if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
992	>	fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");
993
994
995		/* initialize category bounds (clearness index bounds) */
#	Line 877 | Line 1063 | return(value);
1063		/*check the range of epsilon and delta indexes of the perez parametrization*/
1064		void check_parametrization()
1065		{
1066	<	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
1066	>	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup)
1067		{
1068	<	fprintf(stderr,"sky clearness or sky brightness out of range %lf\t %lf\n", skyclearness, skybrightness);
1069	<	exit(1);
1068	>
1069	>	/*  limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */
1070	>	if (skyclearness<skyclearinf){
1071	>	if (suppress_warnings==0)
1072	>	/* fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */
1073	>	skyclearness=skyclearinf;
1074		}
1075	+	if (skyclearness>skyclearsup){
1076	+	if (suppress_warnings==0)
1077	+	/* fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */
1078	+	skyclearness=skyclearsup-0.1;
1079	+	}
1080	+	if (skybrightness<skybriginf){
1081	+	if (suppress_warnings==0)
1082	+	/* fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */
1083	+	skybrightness=skybriginf;
1084	+	}
1085	+	if (skybrightness>skybrigsup){
1086	+	if (suppress_warnings==0)
1087	+	/* fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */
1088	+	skybrightness=skybrigsup;
1089	+	}
1090	+
1091	+	return; }
1092		else return;
1093		}
1094
1095
1096	<	/* likelihood of the direct and diffuse components */
1096	>	/* validity of the direct and diffuse components */
1097		void    check_illuminances()
1098		{
1099	<	if (!( (directilluminance>=0) && (directilluminance<=solar_constant_l*1000) && (diffusilluminance>0) ))
1100	<	{
1101	<	fprintf(stderr,"direct or diffuse illuminances out of range\n");
895	<	exit(1);
1099	>	if (directilluminance < 0) {
1100	>	fprintf(stderr,"WARNING: direct illuminance < 0. Using 0.0\n");
1101	>	directilluminance = 0.0;
1102		}
1103	<	return;
1103	>	if (diffuseilluminance < 0) {
1104	>	fprintf(stderr,"WARNING: diffuse illuminance < 0. Using 0.0\n");
1105	>	diffuseilluminance = 0.0;
1106	>	}
1107	>	if (directilluminance > solar_constant_l*1000.0) {
1108	>	fprintf(stderr,"ERROR: direct illuminance exceeds solar constant\n");
1109	>	exit(1);
1110	>	}
1111		}
1112
1113
1114		void    check_irradiances()
1115		{
1116	<	if (!( (directirradiance>=0) && (directirradiance<=solar_constant_e) && (diffusirradiance>0) ))
1117	<	{
1118	<	fprintf(stderr,"direct or diffuse irradiances out of range\n");
1119	<	exit(1);
1120	<	}
1121	<	return;
1116	>	if (directirradiance < 0) {
1117	>	fprintf(stderr,"WARNING: direct irradiance < 0. Using 0.0\n");
1118	>	directirradiance = 0.0;
1119	>	}
1120	>	if (diffuseirradiance < 0) {
1121	>	fprintf(stderr,"WARNING: diffuse irradiance < 0. Using 0.0\n");
1122	>	diffuseirradiance = 0.0;
1123	>	}
1124	>	if (directirradiance > solar_constant_e) {
1125	>	fprintf(stderr,"ERROR: direct irradiance exceeds solar constant\n");
1126	>	exit(1);
1127	>	}
1128		}
1129
1130
#	Line 915 | Line 1134 | double sky_brightness()
1134		{
1135		double value;
1136
1137	<	value = diffusirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
1137	>	value = diffuseirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
1138
1139		return(value);
1140		}
#	Line 924 | Line 1143 | return(value);
1143		/* Perez sky's clearness */
1144		double sky_clearness()
1145		{
1146	<	double value;
1146	>	double value;
1147
1148	<	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) ;
1148	>	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) ;
1149
1150	<	return(value);
1150	>	return(value);
1151		}
1152
1153
1154
1155		/* diffus horizontal irradiance from Perez sky's brightness */
1156	<	double diffus_irradiance_from_sky_brightness()
1156	>	double diffuse_irradiance_from_sky_brightness()
1157		{
1158		double value;
1159
#	Line 949 | Line 1168 | double direct_irradiance_from_sky_clearness()
1168		{
1169		double value;
1170
1171	<	value = diffus_irradiance_from_sky_brightness();
1171	>	value = diffuse_irradiance_from_sky_brightness();
1172		value = value * ( (skyclearness-1) * (1+1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) );
1173
1174		return(value);
1175		}
1176
1177
1178	<	void illu_to_irra_index(void)
1178	>
1179	>
1180	>	void illu_to_irra_index()
1181		{
1182	<	double  test1=0.1, test2=0.1;
1182	>	double  test1=0.1, test2=0.1, d_eff;
1183		int     counter=0;
1184
1185	<	diffusirradiance = diffusilluminance*solar_constant_e/(solar_constant_l*1000);
1185	>	diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l*1000);
1186		directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000);
1187		skyclearness =  sky_clearness();
1188		skybrightness = sky_brightness();
1189	<	if (skyclearness>12) skyclearness=12;
969	<	if (skybrightness<0.05) skybrightness=0.01;
970	<
971	<
972	<	while ( ((fabs(diffusirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
973	<	|| skyclearness>skyclearinf || skyclearness<skyclearsup
974	<	|| skybrightness>skybriginf || skybrightness<skybrigsup )
975	<	&& !(counter==5) )
976	<	{
977	<	/*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
1189	>	check_parametrization();
1190
1191	<	test1=diffusirradiance;
1191	>
1192	>	while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
1193	>	|| (!(skyclearness<skyclearinf || skyclearness>skyclearsup))
1194	>	|| (!(skybrightness<skybriginf || skybrightness>skybrigsup)) )
1195	>	&& !(counter==9) )
1196	>	{
1197	>
1198	>	test1=diffuseirradiance;
1199		test2=directirradiance;
1200		counter++;
1201
1202	<	diffusirradiance = diffusilluminance/glob_h_diffuse_effi_PEREZ();
1203	<	directirradiance = directilluminance/direct_n_effi_PEREZ();
985	<	/*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
1202	>	diffuseirradiance = diffuseilluminance/glob_h_diffuse_effi_PEREZ();
1203	>	d_eff = direct_n_effi_PEREZ();
1204
1205	+
1206	+	if (d_eff < 0.1)
1207	+	directirradiance = 0;
1208	+	else
1209	+	directirradiance = directilluminance/d_eff;
1210	+
1211		skybrightness = sky_brightness();
1212		skyclearness =  sky_clearness();
1213	<	if (skyclearness>12) skyclearness=12;
1214	<	if (skybrightness<0.05) skybrightness=0.01;
1215	<
1216	<	/*fprintf(stderr, "%lf\t %lf\n", skybrightness, skyclearness);*/
993	<
1213	>	check_parametrization();
1214	>
1215	>	/*fprintf(stderr,"skyclearness = %lf, skybrightness = %lf, directirradiance = %lf, diffuseirradiance = %lf\n",skyclearness, skybrightness, directirradiance, diffuseirradiance);*/
1216	>
1217		}
1218
1219
1220		return;
1221		}
1222
1223	<
1001	<	int lect_coeff_perez(char *filename,float **coeff_perez)
1223	>	static int get_numlin(float epsilon)
1224		{
1225	<	FILE *fcoeff_perez;
1226	<	float temp;
1227	<	int i,j;
1228	<
1229	<	if ((fcoeff_perez = frlibopen(filename)) == NULL)
1230	<	{
1231	<	fprintf(stderr,"file %s cannot be opened\n", filename);
1232	<	return 1; /* il y a un probleme de fichier */
1233	<	}
1234	<	else
1235	<	{
1236	<	/*printf("file %s  open\n", filename);*/
1237	<	}
1238	<
1239	<	skip_comments(fcoeff_perez);
1018	<
1019	<	for (i=0;i<8;i++)
1020	<	for (j=0;j<20;j++)
1021	<	{
1022	<	fscanf(fcoeff_perez,"%f",&temp);
1023	<	*(*coeff_perez+i*20+j) = temp;
1024	<	}
1025	<	fclose(fcoeff_perez);
1026	<
1027	<	return 0; /* tout est OK */
1225	>	if (epsilon < 1.065)
1226	>	return 0;
1227	>	else if (epsilon < 1.230)
1228	>	return 1;
1229	>	else if (epsilon < 1.500)
1230	>	return 2;
1231	>	else if (epsilon < 1.950)
1232	>	return 3;
1233	>	else if (epsilon < 2.800)
1234	>	return 4;
1235	>	else if (epsilon < 4.500)
1236	>	return 5;
1237	>	else if (epsilon < 6.200)
1238	>	return 6;
1239	>	return 7;
1240		}
1241
1030	–
1031	–
1242		/* sky luminance perez model */
1243	<	double calc_rel_lum_perez(double dzeta,double gamma,double Z,
1034	<	double epsilon,double Delta,float *coeff_perez)
1243	>	double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[])
1244		{
1245	+
1246		float x[5][4];
1247		int i,j,num_lin;
1248		double c_perez[5];
1249
1250		if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
1251		{
1252	<	fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n");
1252	>	fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez!\n");
1253		exit(1);
1254		}
1255
#	Line 1048 | Line 1258 | double epsilon,double Delta,float *coeff_perez)
1258		{
1259		if ( Delta < 0.2 ) Delta = 0.2;
1260		}
1261	<
1262	<	if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0;
1263	<	if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
1264	<	if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2;
1055	<	if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3;
1056	<	if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4;
1057	<	if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5;
1058	<	if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6;
1059	<	if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7;
1060	<
1261	>
1262	>
1263	>	num_lin = get_numlin(epsilon);
1264	>
1265		for (i=0;i<5;i++)
1266		for (j=0;j<4;j++)
1267		{
1268		x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j);
1269	<	/* printf("x %d %d vaut %f\n",i,j,x[i][j]); */
1269	>	/* fprintf(stderr,"x %d %d vaut %f\n",i,j,x[i][j]); */
1270		}
1271
1272
#	Line 1089 | Line 1293 | double epsilon,double Delta,float *coeff_perez)
1293
1294
1295		/* coefficients for the sky luminance perez model */
1296	<	void coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez)
1296	>	void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[])
1297		{
1298		float x[5][4];
1299		int i,j,num_lin;
#	Line 1105 | Line 1309 | void coeff_lum_perez(double Z, double epsilon, double
1309		{
1310		if ( Delta < 0.2 ) Delta = 0.2;
1311		}
1312	+
1313	+
1314	+	num_lin = get_numlin(epsilon);
1315
1316	<	if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0;
1110	<	if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
1111	<	if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2;
1112	<	if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3;
1113	<	if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4;
1114	<	if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5;
1115	<	if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6;
1116	<	if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7;
1316	>	/*fprintf(stderr,"numlin %d\n", num_lin);*/
1317
1318		for (i=0;i<5;i++)
1319		for (j=0;j<4;j++)
#	Line 1145 | Line 1345 | void coeff_lum_perez(double Z, double epsilon, double
1345		}
1346
1347
1348	+
1349		/* degrees into radians */
1350		double radians(double degres)
1351		{
1352		return degres*M_PI/180.0;
1353		}
1354
1355	+
1356		/* radian into degrees */
1357		double degres(double radians)
1358		{
1359		return radians/M_PI*180.0;
1360		}
1361
1362	+
1363		/* calculation of the angles dzeta and gamma */
1364		void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z)
1365		{
#	Line 1173 | Line 1376 | void theta_phi_to_dzeta_gamma(double theta,double phi,
1376		}
1377
1378
1176	–	/********************************************************************************/
1177	–	/*      Fonction: theta_ordered                                                 */
1178	–	/*                                                                              */
1179	–	/*      In: char *filename                                                      */
1180	–	/*                                                                              */
1181	–	/*      Out: float *                                                            */
1182	–	/*                                                                              */
1183	–	/*      Update: 29/08/93                                                        */
1184	–	/*                                                                              */
1185	–	/*      Rem: theta en degres                                                    */
1186	–	/*                                                                              */
1187	–	/*      But: fournit les valeurs de theta du fichier d'entree a la memoire      */
1188	–	/*                                                                              */
1189	–	/********************************************************************************/
1190	–	float *theta_ordered(char *filename)
1191	–	{
1192	–	int i;
1193	–	float buffer,*ptr;
1194	–	FILE *file_in;
1379
1196	–	if ( (file_in = frlibopen(filename)) == NULL )
1197	–	{
1198	–	fprintf(stderr,"Cannot open file %s in function theta_ordered\n",filename);
1199	–	exit(1);
1200	–	}
1201	–
1202	–	skip_comments(file_in);
1203	–
1204	–	if ( (ptr = malloc(145*sizeof(float))) == NULL )
1205	–	{
1206	–	fprintf(stderr,"Out of memory in function theta_ordered\n");
1207	–	exit(1);
1208	–	}
1209	–
1210	–	for (i=0;i<145;i++)
1211	–	{
1212	–	fscanf(file_in,"%f",&buffer);
1213	–	*(ptr+i) = buffer;
1214	–	fscanf(file_in,"%f",&buffer);
1215	–	}
1216	–
1217	–	fclose(file_in);
1218	–	return ptr;
1219	–	}
1220	–
1221	–
1222	–	/********************************************************************************/
1223	–	/*      Fonction: phi_ordered                                                   */
1224	–	/*                                                                              */
1225	–	/*      In: char *filename                                                      */
1226	–	/*                                                                              */
1227	–	/*      Out: float *                                                            */
1228	–	/*                                                                              */
1229	–	/*      Update: 29/08/93                                                        */
1230	–	/*                                                                              */
1231	–	/*      Rem: valeurs de Phi en DEGRES                                           */
1232	–	/*                                                                              */
1233	–	/*      But: mettre les angles contenus dans le fichier d'entree dans la memoire */
1234	–	/*                                                                              */
1235	–	/********************************************************************************/
1236	–	float *phi_ordered(char *filename)
1237	–	{
1238	–	int i;
1239	–	float buffer,*ptr;
1240	–	FILE *file_in;
1241	–
1242	–	if ( (file_in = frlibopen(filename)) == NULL )
1243	–	{
1244	–	fprintf(stderr,"Cannot open file %s in function phi_ordered\n",filename);
1245	–	exit(1);
1246	–	}
1247	–
1248	–	skip_comments(file_in);
1249	–
1250	–	if ( (ptr = malloc(145*sizeof(float))) == NULL )
1251	–	{
1252	–	fprintf(stderr,"Out of memory in function phi_ordered");
1253	–	exit(1);
1254	–	}
1255	–
1256	–	for (i=0;i<145;i++)
1257	–	{
1258	–	fscanf(file_in,"%f",&buffer);
1259	–	fscanf(file_in,"%f",&buffer);
1260	–	*(ptr+i) = buffer;
1261	–	}
1262	–
1263	–	fclose(file_in);
1264	–	return ptr;
1265	–	}
1266	–
1267	–
1268	–	/********************************************************************************/
1269	–	/*      Fonction: integ_lv                                                      */
1270	–	/*                                                                              */
1271	–	/*      In: float *lv,*theta                                                    */
1272	–	/*          int sun_pos                                                         */
1273	–	/*                                                                              */
1274	–	/*      Out: double                                                             */
1275	–	/*                                                                              */
1276	–	/*      Update: 29/08/93                                                        */
1277	–	/*                                                                              */
1278	–	/*      Rem:                                                                    */
1279	–	/*                                                                              */
1280	–	/*      But: calcul l'integrale de luminance relative sans la dir. du soleil    */
1281	–	/*                                                                              */
1282	–	/********************************************************************************/
1380		double integ_lv(float *lv,float *theta)
1381		{
1382		int i;
1383		double buffer=0.0;
1384	<
1384	>
1385		for (i=0;i<145;i++)
1386	+	{
1387		buffer += (*(lv+i))*cos(radians(*(theta+i)));
1388	<
1388	>	}
1389	>
1390		return buffer*2*M_PI/144;
1292	–
1391		}
1392
1393
1394
1297	–
1298	–
1299	–
1395		/* enter day number(double), return E0 = square(R0/R):  eccentricity correction factor  */
1396
1397		double get_eccentricity()
#	Line 1309 | Line 1404 | double get_eccentricity()
1404		0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle);
1405
1406		return (E0);
1312	–
1407		}
1408
1409
#	Line 1317 | Line 1411 | double get_eccentricity()
1411		double  air_mass()
1412		{
1413		double  m;
1320	–
1414		if (sunzenith>90)
1415		{
1416	<	fprintf(stderr, "solar zenith angle larger than 90� in fuction air_mass():\n the models used are not more valid\n");
1416	>	fprintf(stderr, "Solar zenith angle larger than 90 degrees in function air_mass()\n");
1417		exit(1);
1418		}
1326	–
1419		m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1420		return(m);
1421		}
1330	–
1331	–
1332	–	double get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut)
1333	–
1334	–	{
1335	–
1336	–	double angle;
1337	–
1338	–
1339	–	if (sun_zenith == 0)
1340	–	puts("WARNING: zenith_angle = 0 in function get_angle_sun_vert_plan");
1341	–
1342	–	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) );
1343	–	angle = angle*180/M_PI;
1344	–	return(angle);
1345	–	}
1346	–
1347	–
1348	–
1349	–
1350	–
1422
1423
1424

Diff Legend

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