ViewVC Help

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

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

Comparing ray/src/gen/gendaylit.c (file contents):
Revision 2.3 by greg, Sat Jun 20 21:34:34 2009 UTC vs.
Revision 2.17 by greg, Fri Aug 31 16:01:45 2018 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	+	*  print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018
8		*/
9
10	<
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	<
10	>	#define  _USE_MATH_DEFINES
11		#include  <stdio.h>
12		#include  <string.h>
13		#include  <math.h>
14		#include  <stdlib.h>
24	–	#include  <ctype.h>
15
26	–	#include  "rtio.h"
27	–	#include  "fvect.h"
16		#include  "color.h"
17	+	#include  "sun.h"
18		#include  "paths.h"
19
20	<	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);
20	>	#define  DOT(v1,v2)     (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
21
22		double  normsc();
23
24	<	#define DATFILE         "coeff_perez.dat"
24	>	/*static        char *rcsid="$Header$";*/
25
26	+	float coeff_perez[] = {
27	+	1.3525,-0.2576,-0.2690,-1.4366,-0.7670,0.0007,1.2734,-0.1233,2.8000,0.6004,1.2375,1.000,1.8734,0.6297,
28	+	0.9738,0.2809,0.0356,-0.1246,-0.5718,0.9938,-1.2219,-0.7730,1.4148,1.1016,-0.2054,0.0367,-3.9128,0.9156,
29	+	6.9750,0.1774,6.4477,-0.1239,-1.5798,-0.5081,-1.7812,0.1080,0.2624,0.0672,-0.2190,-0.4285,-1.1000,-0.2515,
30	+	0.8952,0.0156,0.2782,-0.1812,-4.5000,1.1766,24.7219,-13.0812,-37.7000,34.8438,-5.0000,1.5218,3.9229,
31	+	-2.6204,-0.0156,0.1597,0.4199,-0.5562,-0.5484,-0.6654,-0.2672,0.7117,0.7234,-0.6219,-5.6812,2.6297,
32	+	33.3389,-18.3000,-62.2500,52.0781,-3.5000,0.0016,1.1477,0.1062,0.4659,-0.3296,-0.0876,-0.0329,-0.6000,
33	+	-0.3566,-2.5000,2.3250,0.2937,0.0496,-5.6812,1.8415,21.0000,-4.7656,-21.5906,7.2492,-3.5000,-0.1554,
34	+	1.4062,0.3988,0.0032,0.0766,-0.0656,-0.1294,-1.0156,-0.3670,1.0078,1.4051,0.2875,-0.5328,-3.8500,3.3750,
35	+	14.0000,-0.9999,-7.1406,7.5469,-3.4000,-0.1078,-1.0750,1.5702,-0.0672,0.4016,0.3017,-0.4844,-1.0000,
36	+	0.0211,0.5025,-0.5119,-0.3000,0.1922,0.7023,-1.6317,19.0000,-5.0000,1.2438,-1.9094,-4.0000,0.0250,0.3844,
37	+	0.2656,1.0468,-0.3788,-2.4517,1.4656,-1.0500,0.0289,0.4260,0.3590,-0.3250,0.1156,0.7781,0.0025,31.0625,
38	+	-14.5000,-46.1148,55.3750,-7.2312,0.4050,13.3500,0.6234,1.5000,-0.6426,1.8564,0.5636};
39
40
41	<	/* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */
41	>	float defangle_theta[] = {
42	>	84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84,
43	>	84, 84, 84, 84, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72,
44	>	72, 72, 72, 72, 72, 72, 72, 72, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60,
45	>	60, 60, 60, 60, 60, 60, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48,
46	>	48, 48, 48, 48, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 24, 24, 24, 24,
47	>	24, 24, 24, 24, 24, 24, 24, 24, 12, 12, 12, 12, 12, 12, 0};
48	>
49	>	float defangle_phi[] = {
50	>	0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264,
51	>	276, 288, 300, 312, 324, 336, 348, 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180,
52	>	192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 0, 15, 30, 45, 60, 75, 90, 105,
53	>	120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 15, 30, 45, 60, 75,
54	>	90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 20, 40, 60,
55	>	80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 0, 30, 60, 90, 120, 150, 180, 210,
56	>	240, 270, 300, 330, 0, 60, 120, 180, 240, 300, 0};
57	>	/* default values for Berlin */
58	>	float   locus[] = {
59	>	-4.843e9,2.5568e6,0.24282e3,0.23258,-4.843e9,2.5568e6,0.24282e3,0.23258,-1.2848,1.7519,-0.093786};
60	>
61	>
62	>
63	>	/* Perez sky parametrization: epsilon and delta calculations from the direct and diffuse irradiances */
64		double sky_brightness();
65		double sky_clearness();
66
67		/* calculation of the direct and diffuse components from the Perez parametrization */
68	<	double  diffus_irradiance_from_sky_brightness();
68	>	double  diffuse_irradiance_from_sky_brightness();
69		double  direct_irradiance_from_sky_clearness();
70
71	+	/* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : */
72	+	/* input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */
73
52	–	/* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */
74		double  glob_h_effi_PEREZ();
75		double  glob_h_diffuse_effi_PEREZ();
76		double  direct_n_effi_PEREZ();
77	+
78		/*likelihood check of the epsilon, delta, direct and diffuse components*/
79		void    check_parametrization();
80		void    check_irradiances();
81		void    check_illuminances();
82		void    illu_to_irra_index();
83	+	void    print_error_sky();
84
85	<
86	<	/* 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);
85	>	double  calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]);
86	>	void    coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]);
87		double  radians(double degres);
88		double  degres(double radians);
89		void    theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z);
90		double  integ_lv(float *lv,float *theta);
73	–	float   *theta_ordered(char *filename);
74	–	float   *phi_ordered(char *filename);
75	–	void    skip_comments(FILE *fp);
91
92	+	void printdefaults();
93	+	void check_sun_position();
94	+	void computesky();
95	+	void printhead(int ac, char** av);
96	+	void usage_error(char* msg);
97	+	void printsky();
98
99	+	FILE * frlibopen(char* fname);
100
101		/* astronomy and geometry*/
102		double  get_eccentricity();
103		double  air_mass();
82	–	double  get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut);
104
105	<
106	<	/* date*/
105	>	double  solar_sunset(int month, int day);
106	>	double  solar_sunrise(int month, int day);
107	>	double  stadj();
108		int     jdate(int month, int day);
109
88	–
89	–
90	–
91	–
92	–	/* sun calculation constants */
93	–	extern double  s_latitude;
94	–	extern double  s_longitude;
95	–	extern double  s_meridian;
96	–
110		const double    AU = 149597890E3;
111		const double    solar_constant_e = 1367;    /* solar constant W/m^2 */
112	<	const double    solar_constant_l = 127.5;   /* solar constant klux */
112	>	const double    solar_constant_l = 127500;   /* solar constant lux */
113
114		const double    half_sun_angle = 0.2665;
115		const double    half_direct_angle = 2.85;
116
117	<	const double    skyclearinf = 1.000;    /* limitations for the variation of the Perez parameters */
118	<	const double    skyclearsup = 12.1;
117	>	const double    skyclearinf = 1.0;          /* limitations for the variation of the Perez parameters */
118	>	const double    skyclearsup = 12.01;
119		const double    skybriginf = 0.01;
120		const double    skybrigsup = 0.6;
121
#	Line 117 | Line 130 | double  altitude, azimuth;                     /* or solar angles */
130
131
132		/* definition of the sky conditions through the Perez parametrization */
133	<	double  skyclearness, skybrightness;
134	<	double  solarradiance;  /*radiance of the sun disk and of the circumsolar area*/
135	<	double  diffusilluminance, directilluminance, diffusirradiance, directirradiance;
136	<	double  sunzenith, daynumber=150, atm_preci_water=2;
133	>	double  skyclearness = 0;
134	>	double  skybrightness = 0;
135	>	double  solarradiance;
136	>	double  diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance;
137	>	double  sunzenith, daynumber, atm_preci_water=2;
138
139	<	double  diffnormalization, dirnormalization;
139	>	/*double  sunaltitude_border = 0;*/
140	>	double  diffnormalization = 0;
141	>	double  dirnormalization = 0;
142		double  *c_perez;
143
144	<	int     output=0;       /*define the unit of the output (sky luminance or radiance): visible watt=0, solar watt=1, lumen=2*/
145	<	int     input=0;        /*define the input for the calulation*/
144	>	int     output=0;       /* define the unit of the output (sky luminance or radiance): */
145	>	/* visible watt=0, solar watt=1, lumen=2 */
146	>	int     input=0;        /* define the input for the calulation */
147	>	int     color_output=0;
148	>	int     suppress_warnings=0;
149
150		/* default values */
151	<	int  cloudy = 0;                                /* 1=standard, 2=uniform */
152	<	int  dosun = 1;
151	>	int     cloudy = 0;                             /* 1=standard, 2=uniform */
152	>	int     dosun = 1;
153		double  zenithbr = -1.0;
154		double  betaturbidity = 0.1;
155		double  gprefl = 0.2;
156		int     S_INTER=0;
157
158	+
159		/* computed values */
160		double  sundir[3];
161	<	double  groundbr;
161	>	double  groundbr = 0;
162		double  F2;
163		double  solarbr = 0.0;
164		int     u_solar = 0;                            /* -1=irradiance, 1=radiance */
165	+	float   timeinterval = 0;
166
167	<	char  *progname;
168	<	char  errmsg[128];
167	>	char    *progname;
168	>	char    errmsg[128];
169
170	+	double  st;
171
172	<	main(argc, argv)
173	<	int  argc;
152	<	char  *argv[];
172	>
173	>	int main(int argc, char** argv)
174		{
175		int  i;
176
177		progname = argv[0];
178		if (argc == 2 && !strcmp(argv[1], "-defaults")) {
179		printdefaults();
180	<	exit(0);
180	>	return 0;
181		}
182		if (argc < 4)
183	<	userror("arg count");
183	>	usage_error("arg count");
184		if (!strcmp(argv[1], "-ang")) {
185		altitude = atof(argv[2]) * (M_PI/180);
186		azimuth = atof(argv[3]) * (M_PI/180);
#	Line 167 | Line 188 | char  *argv[];
188		} else {
189		month = atoi(argv[1]);
190		if (month < 1 || month > 12)
191	<	userror("bad month");
191	>	usage_error("bad month");
192		day = atoi(argv[2]);
193		if (day < 1 || day > 31)
194	<	userror("bad day");
194	>	usage_error("bad day");
195		hour = atof(argv[3]);
196		if (hour < 0 || hour >= 24)
197	<	userror("bad hour");
197	>	usage_error("bad hour");
198		tsolar = argv[3][0] == '+';
199		}
200		for (i = 4; i < argc; i++)
#	Line 183 | Line 204 | char  *argv[];
204		cloudy = 0;
205		dosun = argv[i][0] == '+';
206		break;
186	–	case 'r':
207		case 'R':
208		u_solar = argv[i][1] == 'R' ? -1 : 1;
209		solarbr = atof(argv[++i]);
#	Line 192 | Line 212 | char  *argv[];
212		cloudy = argv[i][0] == '+' ? 2 : 1;
213		dosun = 0;
214		break;
215	+	case 'C':
216	+	if (argv[i][2] == 'I' && argv[i][3] == 'E' ) {
217	+	locus[0] = -4.607e9;
218	+	locus[1] = 2.9678e6;
219	+	locus[2] = 0.09911e3;
220	+	locus[3] = 0.244063;
221	+	locus[4] = -2.0064e9;
222	+	locus[5] = 1.9018e6;
223	+	locus[6] = 0.24748e3;
224	+	locus[7] = 0.23704;
225	+	locus[8] = -3.0;
226	+	locus[9] = 2.87;
227	+	locus[10] = -0.275;
228	+	}else{ color_output = 1;
229	+	}
230	+	break;
231	+	case 'l':
232	+	locus[0] = atof(argv[++i]);
233	+	locus[1] = atof(argv[++i]);
234	+	locus[2] = atof(argv[++i]);
235	+	locus[3] = atof(argv[++i]);
236	+	locus[4] = locus[0];
237	+	locus[5] = locus[1];
238	+	locus[6] = locus[2];
239	+	locus[7] = locus[3];
240	+	locus[8] = atof(argv[++i]);
241	+	locus[9] = atof(argv[++i]);
242	+	locus[10] = atof(argv[++i]);
243	+	break;
244	+
245		case 't':
246		betaturbidity = atof(argv[++i]);
247		break;
248	+	case 'w':
249	+	suppress_warnings = 1;
250	+	break;
251		case 'b':
252		zenithbr = atof(argv[++i]);
253		break;
#	Line 210 | Line 263 | char  *argv[];
263		case 'm':
264		s_meridian = atof(argv[++i]) * (M_PI/180);
265		break;
213	–
266
267		case 'O':
268	<	output = atof(argv[++i]);       /*define the unit of the output of the program :
269	<	sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm)
270	<	default is set to 0*/
268	>	output = atof(argv[++i]);       /*define the unit of the output of the program:
269	>	sky and sun luminance/radiance
270	>	(0==W visible, 1==W solar radiation, 2==lm) */
271		break;
272
273		case 'P':
#	Line 227 | Line 279 | char  *argv[];
279		case 'W':                                       /* direct normal Irradiance [W/m^2] */
280		input = 1;                              /* diffuse horizontal Irrad. [W/m^2] */
281		directirradiance = atof(argv[++i]);
282	<	diffusirradiance = atof(argv[++i]);
282	>	diffuseirradiance = atof(argv[++i]);
283		break;
284
285		case 'L':                                       /* direct normal Illuminance [Lux] */
286		input = 2;                              /* diffuse horizontal Ill. [Lux] */
287		directilluminance = atof(argv[++i]);
288	<	diffusilluminance = atof(argv[++i]);
288	>	diffuseilluminance = atof(argv[++i]);
289		break;
290
291		case 'G':                                       /* direct horizontal Irradiance [W/m^2] */
292		input = 3;                              /* diffuse horizontal Irrad. [W/m^2] */
293		directirradiance = atof(argv[++i]);
294	<	diffusirradiance = atof(argv[++i]);
294	>	diffuseirradiance = atof(argv[++i]);
295		break;
244	–
296
297	+	case 'E':                                       /* Erbs model based on the */
298	+	input = 4;                              /* global-horizontal irradiance [W/m^2] */
299	+	globalirradiance = atof(argv[++i]);
300	+	break;
301	+
302	+	case 'i':
303	+	timeinterval = atof(argv[++i]);
304	+	break;
305	+
306	+
307		default:
308		sprintf(errmsg, "unknown option: %s", argv[i]);
309	<	userror(errmsg);
309	>	usage_error(errmsg);
310		}
311		else
312	<	userror("bad option");
312	>	usage_error("bad option");
313
314	<	if (fabs(s_meridian-s_longitude) > 30*M_PI/180)
315	<	fprintf(stderr,
255	<	"%s: warning: %.1f hours btwn. standard meridian and longitude\n",
314	>	if (month && !tsolar && fabs(s_meridian-s_longitude) > 45*M_PI/180)
315	>	fprintf(stderr,"%s: warning: %.1f hours btwn. standard meridian and longitude\n",
316		progname, (s_longitude-s_meridian)*12/M_PI);
317
318
319	<	/* allocation dynamique de memoire pour les pointeurs */
320	<	if ( (c_perez = malloc(5*sizeof(double))) == NULL )
321	<	{
262	<	fprintf(stderr,"Out of memory error in function main !");
263	<	exit(1);
264	<	}
319	>	/* dynamic memory allocation for the pointers */
320	>	if ( (c_perez = calloc(5, sizeof(double))) == NULL )
321	>	{ fprintf(stderr,"Out of memory error in function main"); return 1; }
322
323	<
323	>
324		printhead(argc, argv);
268	–
325		computesky();
326		printsky();
327	+	return 0;
328
272	–	exit(0);
329		}
330
331
332	<	computesky()                    /* compute sky parameters */
332	>
333	>
334	>
335	>	void computesky()
336		{
337
338	<	/* new variables */
339	<	int     j, i;
340	<	float   *lv_mod;  /* 145 luminance values*/
341	<	/* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */
283	<	float   *theta_o, *phi_o, *coeff_perez;
338	>	int     j;
339	>
340	>	float   *lv_mod;  /* 145 luminance values */
341	>	float   *theta_o, *phi_o;
342		double  dzeta, gamma;
285	–	double  diffusion;
343		double  normfactor;
344	+	double  erbs_s0, erbs_kt;
345
346
289	–
347		/* compute solar direction */
348	<
348	>
349		if (month) {                    /* from date and time */
350		int  jd;
351	<	double  sd, st;
351	>	double  sd;
352
353		jd = jdate(month, day);         /* Julian date */
354		sd = sdec(jd);                  /* solar declination */
#	Line 299 | Line 356 | computesky()                   /* compute sky parameters */
356		st = hour;
357		else
358		st = hour + stadj(jd);
359	+
360	+
361	+	if(timeinterval) {
362	+
363	+	if(timeinterval<0) {
364	+	fprintf(stderr, "time interval negative\n");
365	+	exit(1);
366	+	}
367	+
368	+	if(fabs(solar_sunrise(month,day)-st)<=timeinterval/120) {
369	+	st= (st+timeinterval/120+solar_sunrise(month,day))/2;
370	+	if(suppress_warnings==0)
371	+	{ fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); }
372	+	}
373	+
374	+	if(fabs(solar_sunset(month,day)-st)<timeinterval/120) {
375	+	st= (st-timeinterval/120+solar_sunset(month,day))/2;
376	+	if(suppress_warnings==0)
377	+	{ fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); }
378	+	}
379	+
380	+	if((st<solar_sunrise(month,day)-timeinterval/120) || (st>solar_sunset(month,day)+timeinterval/120)) {
381	+	if(suppress_warnings==0)
382	+	{ fprintf(stderr, "Warning: sun position too low, printing error sky at %d %d %.3f\n",month,day,hour); }
383	+	altitude = salt(sd, st);
384	+	azimuth = sazi(sd, st);
385	+	print_error_sky();
386	+	exit(0);
387	+	}
388	+	}
389	+	else
390	+
391	+	if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) {
392	+	if(suppress_warnings==0)
393	+	{ fprintf(stderr, "Warning: sun altitude below zero at time step %i %i %.2f, printing error sky\n",month,day,hour); }
394	+	altitude = salt(sd, st);
395	+	azimuth = sazi(sd, st);
396	+	print_error_sky();
397	+	exit(0);
398	+	}
399	+
400		altitude = salt(sd, st);
401		azimuth = sazi(sd, st);
402
403		daynumber = (double)jdate(month, day);
404	<
404	>
405		}
406	+
407	+
408	+
409	+
410	+
411		if (!cloudy && altitude > 87.*M_PI/180.) {
412	<	fprintf(stderr,
412	>
413	>	if (suppress_warnings==0) {
414	>	fprintf(stderr,
415		"%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n",
416		progname);
417	<	printf(
313	<	"# warning - sun too close to zenith, reducing altitude to 87 degrees\n");
417	>	}
418		altitude = 87.*M_PI/180.;
419		}
420	+
421	+
422	+
423		sundir[0] = -sin(azimuth)*cos(altitude);
424		sundir[1] = -cos(azimuth)*cos(altitude);
425		sundir[2] = sin(altitude);
#	Line 320 | Line 427 | computesky()                   /* compute sky parameters */
427
428		/* calculation for the new functions */
429		sunzenith = 90 - altitude*180/M_PI;
430	<
324	<
430	>
431
432	<	/* compute the inputs for the calculation of the light distribution over the sky*/
433	<	if (input==0)
432	>	/* compute the inputs for the calculation of the light distribution over the sky*/
433	>	if (input==0)           /* P */
434		{
435		check_parametrization();
436	<	diffusirradiance = diffus_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
436	>	diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
437		directirradiance = direct_irradiance_from_sky_clearness();
438		check_irradiances();
439
440		if (output==0 || output==2)
441		{
442	<	diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
442	>	diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
443		directilluminance = directirradiance*direct_n_effi_PEREZ();
444		check_illuminances();
445		}
446		}
447
448
449	<	else if (input==1)
449	>	else if (input==1)      /* W */
450		{
451		check_irradiances();
452		skybrightness = sky_brightness();
453		skyclearness =  sky_clearness();
454	+
455		check_parametrization();
456	<
456	>
457		if (output==0 || output==2)
458		{
459	<	diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
459	>	diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
460		directilluminance = directirradiance*direct_n_effi_PEREZ();
461		check_illuminances();
462		}
#	Line 357 | Line 464 | computesky()                   /* compute sky parameters */
464		}
465
466
467	<	else if (input==2)
467	>	else if (input==2)      /* L */
468		{
469		check_illuminances();
470		illu_to_irra_index();
#	Line 365 | Line 472 | computesky()                   /* compute sky parameters */
472		}
473
474
475	<	else if (input==3)
475	>	else if (input==3)      /* G */
476		{
477		if (altitude<=0)
478		{
479	<	fprintf(stderr, "solar zenith angle larger than 90� \n the models used are not more valid\n");
480	<	exit(1);
479	>	if (suppress_warnings==0)
480	>	fprintf(stderr, "Warning: sun altitude < 0, proceed with irradiance values of zero\n");
481	>	directirradiance = 0;
482	>	diffuseirradiance = 0;
483	>	} else {
484	>
485	>	directirradiance=directirradiance/sin(altitude);
486		}
487	<
376	<	directirradiance=directirradiance/sin(altitude);
487	>
488		check_irradiances();
489		skybrightness = sky_brightness();
490		skyclearness =  sky_clearness();
#	Line 381 | Line 492 | computesky()                   /* compute sky parameters */
492
493		if (output==0 || output==2)
494		{
495	<	diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
495	>	diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
496		directilluminance = directirradiance*direct_n_effi_PEREZ();
497		check_illuminances();
498		}
499
500		}
501
391	–
392	–	else    {fprintf(stderr,"error in giving the input arguments"); exit(1);}
502
503	+	else if (input==4)      /* E */         /* Implementation of the Erbs model. W.Sprenger (04/13) */
504	+	{
505	+
506	+	if (altitude<=0)
507	+	{
508	+	if (suppress_warnings==0 && globalirradiance > 50)
509	+	fprintf(stderr, "Warning: global irradiance higher than 50 W/m^2 while the sun altitude is lower than zero\n");
510	+	globalirradiance = 0; diffuseirradiance = 0; directirradiance = 0;
511	+
512	+	} else {
513	+
514	+	erbs_s0 = solar_constant_e*get_eccentricity()*sin(altitude);
515	+
516	+	if (globalirradiance>erbs_s0)
517	+	{
518	+	if (suppress_warnings==0)
519	+	fprintf(stderr, "Warning: global irradiance is higher than the time-dependent solar constant s0\n");
520	+	globalirradiance=erbs_s0*0.999;
521	+	}
522	+
523	+	erbs_kt=globalirradiance/erbs_s0;
524	+
525	+	if (erbs_kt<=0.22)      diffuseirradiance=globalirradiance*(1-0.09*erbs_kt);
526	+	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));
527	+	else if (erbs_kt<1)     diffuseirradiance=globalirradiance*(0.165);
528	+
529	+	directirradiance=globalirradiance-diffuseirradiance;
530	+
531	+	printf("# erbs_s0, erbs_kt, irr_dir_h, irr_diff: %.3f %.3f %.3f %.3f\n", erbs_s0, erbs_kt, directirradiance, diffuseirradiance);
532	+	printf("# WARNING: the -E option is only recommended for a rough estimation!\n");
533	+
534	+	directirradiance=directirradiance/sin(altitude);
535	+
536	+	}
537	+
538	+	check_irradiances();
539	+	skybrightness = sky_brightness();
540	+	skyclearness =  sky_clearness();
541	+	check_parametrization();
542
543	+	if (output==0 || output==2)
544	+	{
545	+	diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
546	+	directilluminance = directirradiance*direct_n_effi_PEREZ();
547	+	check_illuminances();
548	+	}
549	+
550	+	}
551	+
552	+
553	+
554
555	<	/* normalization factor for the relative sky luminance distribution, diffuse part*/
555	>	else    { fprintf(stderr,"error at the input arguments"); exit(1); }
556
398	–	/* allocation dynamique de memoire pour les pointeurs */
399	–	if ( (coeff_perez = malloc(8*20*sizeof(float))) == NULL )
400	–	{
401	–	fprintf(stderr,"Out of memory error in function main !");
402	–	exit(1);
403	–	}
557
558	<	/* read the coefficients for the Perez sky luminance model */
559	<	if (lect_coeff_perez(DATFILE, &coeff_perez) > 0)
560	<	{
408	<	fprintf(stderr,"lect_coeff_perez does not work\n");
409	<	exit(2);
410	<	}
411	<
558	>
559	>	/* normalization factor for the relative sky luminance distribution, diffuse part*/
560	>
561		if ( (lv_mod = malloc(145*sizeof(float))) == NULL)
562		{
563		fprintf(stderr,"Out of memory in function main");
#	Line 416 | Line 565 | computesky()                   /* compute sky parameters */
565		}
566
567		/* read the angles */
568	<	theta_o = theta_ordered("defangle.dat");
569	<	phi_o = phi_ordered("defangle.dat");
568	>	theta_o = defangle_theta;
569	>	phi_o = defangle_phi;
570	>
571
572	<	/* parameters for the perez model */
572	>	/* parameters for the perez model */
573		coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez);
574
575	<	/*calculation of the modelled luminance */
575	>
576	>
577	>	/*calculation of the modelled luminance */
578		for (j=0;j<145;j++)
579		{
580		theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith));
581	+
582		*(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez);
583	<	/*printf("theta, phi, lv_mod %lf\t %lf\t %lf\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j));*/
583	>
584	>	/* fprintf(stderr,"theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); */
585		}
586	<
586	>
587		/* integration of luminance for the normalization factor, diffuse part of the sky*/
588	+
589		diffnormalization = integ_lv(lv_mod, theta_o);
435	–	/*printf("perez integration %lf\n", diffnormalization);*/
590
437	–
591
592
593	<	/*normalization coefficient in lumen or in watt*/
593	>	/*normalization coefficient in lumen or in watt*/
594		if (output==0)
595		{
596	<	diffnormalization = diffusilluminance/diffnormalization/WHTEFFICACY;
596	>	diffnormalization = diffuseilluminance/diffnormalization/WHTEFFICACY;
597		}
598		else if (output==1)
599		{
600	<	diffnormalization = diffusirradiance/diffnormalization;
600	>	diffnormalization = diffuseirradiance/diffnormalization;
601		}
602		else if (output==2)
603		{
604	<	diffnormalization = diffusilluminance/diffnormalization;
604	>	diffnormalization = diffuseilluminance/diffnormalization;
605		}
606
607	<	else    {fprintf(stderr,"output argument : wrong number"); exit(1);}
607	>	else    {fprintf(stderr,"Wrong output specification.\n"); exit(1);}
608
609
610
611
612	<	/* calculation for the solar source */
612	>	/* calculation for the solar source */
613		if (output==0)
614		solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY;
615
#	Line 466 | Line 619 | computesky()                   /* compute sky parameters */
619		else
620		solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
621
469	–
622
623
624	<	/* Compute the ground radiance */
625	<	zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
626	<	zenithbr*=diffnormalization;
627	<	/*
628	<	fprintf(stderr, "gendaylit : the actual zenith radiance(W/m^2/sr) or luminance(cd/m^2) is : %.0lf\n", zenithbr);
477	<	*/
478	<
479	<	if (skyclearness==1)
624	>	/* Compute the ground radiance */
625	>	zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
626	>	zenithbr*=diffnormalization;
627	>
628	>	if (skyclearness==1)
629		normfactor = 0.777778;
630
631	<	if (skyclearness>=6)
631	>	if (skyclearness>=6)
632		{
633		F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]);
634		normfactor = normsc()/F2/M_PI;
635		}
636
637	<	if ( (skyclearness>1) && (skyclearness<6) )
637	>	if ( (skyclearness>1) && (skyclearness<6) )
638		{
639		S_INTER=1;
640		F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude));
641		normfactor = normsc()/F2/M_PI;
642		}
643
644	<	groundbr = zenithbr*normfactor;
496	<	printf("# Ground ambient level: %.1f\n", groundbr);
644	>	groundbr = zenithbr*normfactor;
645
646	<	if (dosun&&(skyclearness>1))
647	<	groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];
646	>	if (dosun&&(skyclearness>1))
647	>	groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];
648
649	<	groundbr *= gprefl;
649	>	groundbr *= gprefl;
650
651
652	+
653	+	if(*(c_perez+1)>0)
654	+	{
655	+	if(suppress_warnings==0)
656	+	{  fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1));}
657	+	print_error_sky();
658	+	exit(0);
659	+	}
660
661	+
662		return;
663		}
664
#	Line 509 | Line 666 | return;
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	<	printsky()                      /* print out sky */
678	>
679	>
680	>	double solar_sunrise(int month,int day)
681		{
682	+	float W;
683	+	extern double s_latitude;
684	+	W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day))));
685	+	return(12-(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15));
686	+	}
687	+
688	+
689	+
690	+
691	+	void printsky()
692	+	{
693	+
694	+	printf("# Local solar time: %.2f\n", st);
695	+	printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI);
696	+
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	<	}
525	<
526	<	if (dosun&&(skyclearness==1))
527	<	{
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	>	/* print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018 */
715	>	if  (color_output==1 && skyclearness < 4.5 && skyclearness >1.065 )
716	>	{
717	>	fprintf(stderr, "       warning: sky clearness(epsilon)= %f \n",skyclearness);
718	>	fprintf(stderr, "       warning: intermediate sky!! \n");
719	>	fprintf(stderr, "       warning: color model for intermediate sky pending  \n");
720	>	fprintf(stderr, "       warning: no color output ! \n");
721	>	color_output=0;
722	>	}
723	>	if (color_output==1)
724	>	{
725	>	printf("\nvoid colorfunc skyfunc\n");
726	>	printf("4 skybright_r skybright_g skybright_b perezlum_c.cal\n");
727	>	printf("0\n");
728	>	printf("22 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n", diffnormalization, groundbr,
729	>	*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
730	>	sundir[0], sundir[1], sundir[2],skyclearness,locus[0],locus[1],locus[2],locus[3],locus[4],locus[5],locus[6],locus[7],locus[8],locus[9],locus[10]);
731	>	}else{
732	>	printf("\nvoid brightfunc skyfunc\n");
733	>	printf("2 skybright perezlum.cal\n");
734	>	printf("0\n");
735	>	printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr,
736	>	*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
737	>	sundir[0], sundir[1], sundir[2]);
738	>	}
739
740	+	}
741
742	+
743	+
744	+	void print_error_sky()
745	+	{
746	+
747	+
748	+	sundir[0] = -sin(azimuth)*cos(altitude);
749	+	sundir[1] = -cos(azimuth)*cos(altitude);
750	+	sundir[2] = sin(altitude);
751	+
752	+	printf("# Local solar time: %.2f\n", st);
753	+	printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI);
754	+
755		printf("\nvoid brightfunc skyfunc\n");
756		printf("2 skybright perezlum.cal\n");
757		printf("0\n");
758	<	printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr,
541	<	*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
542	<	sundir[0], sundir[1], sundir[2]);
758	>	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]);
759		}
760	+
761
762
763	<	printdefaults()                 /* print default values */
763	>
764	>
765	>	void printdefaults()                    /* print default values */
766		{
767		printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
768		if (zenithbr > 0.0)
#	Line 556 | Line 775 | printdefaults()                        /* print default values */
775		}
776
777
778	<	userror(msg)                    /* print usage error and quit */
779	<	char  *msg;
778	>
779	>
780	>	void usage_error(char* msg)                     /* print usage error and quit */
781		{
782		if (msg != NULL)
783	<	fprintf(stderr, "%s: Use error - %s\n", progname, msg);
784	<	fprintf(stderr, "Usage: %s month day hour [-P|-W|-L] direct_value diffus_value [options]\n", progname);
785	<	fprintf(stderr, "or   : %s -ang altitude azimuth [-P|-W|-L] direct_value diffus_value [options]\n", progname);
783	>	fprintf(stderr, "%s: Use error - %s\n\n", progname, msg);
784	>	fprintf(stderr, "Usage: %s      month day hour    [...]\n", progname);
785	>	fprintf(stderr, "   or: %s -ang altitude azimuth  [...]\n", progname);
786	>	fprintf(stderr, "               followed by:      -P          epsilon delta [options]\n");
787	>	fprintf(stderr, "                        or:      [-W|-L|-G]  direct_value diffuse_value [options]\n");
788	>	fprintf(stderr, "                        or:      -E          global_irradiance [options]\n\n");
789	>	fprintf(stderr, "       Description:\n");
790		fprintf(stderr, "       -P epsilon delta  (these are the Perez parameters) \n");
791		fprintf(stderr, "       -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
792		fprintf(stderr, "       -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n");
793		fprintf(stderr, "       -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
794	+	fprintf(stderr, "       -E global-horizontal-irradiance (W/m^2)\n\n");
795	+	fprintf(stderr, "       Output specification with option:\n");
796		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");
797	+	fprintf(stderr, "       gendaylit version 2.5 (2018/04/18)  \n\n");
798		exit(1);
799		}
800
801
802
803	<	double
804	<	normsc()                        /* compute normalization factor (E0*F2/L0) */
803	>
804	>	double normsc()           /* compute normalization factor (E0*F2/L0) */
805		{
806		static double  nfc[2][5] = {
807		/* clear sky approx. */
#	Line 597 | Line 824 | normsc()                       /* compute normalization factor (E0*F2/L0)
824
825
826
827	<	printhead(ac, av)               /* print command header */
828	<	register int  ac;
829	<	register char  **av;
827	>
828	>
829	>	void printhead(int ac, char** av)               /* print command header */
830		{
831		putchar('#');
832		while (ac--) {
#	Line 612 | Line 839 | register char  **av;
839
840
841
615	–	void
616	–	skip_comments(FILE *fp)         /* skip comments in file */
617	–	{
618	–	int     c;
619	–
620	–	while ((c = getc(fp)) != EOF)
621	–	if (c == '#') {
622	–	while ((c = getc(fp)) != EOF)
623	–	if (c == '\n')
624	–	break;
625	–	} else if (!isspace(c)) {
626	–	ungetc(c, fp);
627	–	break;
628	–	}
629	–	}
842
843
632	–
844		/* Perez models */
845
846		/* Perez global horizontal luminous efficacy model */
#	Line 639 | Line 850 | double glob_h_effi_PEREZ()
850		double  value;
851		double  category_bounds[10], a[10], b[10], c[10], d[10];
852		int     category_total_number, category_number, i;
853	<
854	<
855	<	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
856	<	fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
857	<
853	>
854	>	check_parametrization();
855	>
856	>
857	>	/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
858	>	fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n"); */
859	>
860	>
861		/* initialize category bounds (clearness index bounds) */
862
863		category_total_number = 8;
#	Line 698 | Line 912 | fprintf(stderr, "Warning : skyclearness or skybrightne
912
913
914
701	–
915		for (i=1; i<=category_total_number; i++)
916		{
917		if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
#	Line 712 | Line 925 | fprintf(stderr, "Warning : skyclearness or skybrightne
925		}
926
927
928	+
929	+
930		/* global horizontal diffuse efficacy model, according to PEREZ */
931		double glob_h_diffuse_effi_PEREZ()
932		{
#	Line 719 | Line 934 | double glob_h_diffuse_effi_PEREZ()
934		double  category_bounds[10], a[10], b[10], c[10], d[10];
935		int     category_total_number, category_number, i;
936
937	+	check_parametrization();
938
939	<	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
940	<	fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
941	<
939	>
940	>	/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
941	>	fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */
942	>
943		/* initialize category bounds (clearness index bounds) */
944
945		category_total_number = 8;
946
947	+	//XXX:  category_bounds > 0.1
948		category_bounds[1] = 1;
949		category_bounds[2] = 1.065;
950		category_bounds[3] = 1.230;
#	Line 777 | Line 995 | fprintf(stderr, "Warning : skyclearness or skybrightne
995
996
997
998	<
998	>	category_number = -1;
999		for (i=1; i<=category_total_number; i++)
1000		{
1001		if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
1002		category_number = i;
1003		}
1004
1005	+	if (category_number == -1) {
1006	+	if (suppress_warnings==0)
1007	+	fprintf(stderr, "Warning: sky clearness (= %.3f) too high, printing error sky\n", skyclearness);
1008	+	print_error_sky();
1009	+	exit(0);
1010	+	}
1011	+
1012	+
1013		value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*cos(sunzenith*M_PI/180) +
1014		d[category_number]*log(skybrightness);
1015
1016		return(value);
1017	+
1018		}
1019
1020
1021	+
1022	+
1023	+
1024	+
1025		/* direct normal efficacy model, according to PEREZ */
1026
1027		double direct_n_effi_PEREZ()
#	Line 801 | Line 1032 | double         category_bounds[10], a[10], b[10], c[10], d[10
1032		int     category_total_number, category_number, i;
1033
1034
1035	<	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
1036	<	fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
1035	>	/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
1036	>	fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");*/
1037
1038
1039		/* initialize category bounds (clearness index bounds) */
#	Line 876 | Line 1107 | return(value);
1107		/*check the range of epsilon and delta indexes of the perez parametrization*/
1108		void check_parametrization()
1109		{
1110	<	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
1110	>
1111	>	if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup)
1112		{
1113	<	fprintf(stderr,"sky clearness or sky brightness out of range %lf\t %lf\n", skyclearness, skybrightness);
1114	<	exit(1);
1113	>
1114	>	/*  limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */
1115	>
1116	>	if (skyclearness<skyclearinf){
1117	>	/* if (suppress_warnings==0)
1118	>	fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */
1119	>	skyclearness=skyclearinf;
1120		}
1121	+	if (skyclearness>skyclearsup){
1122	+	/* if (suppress_warnings==0)
1123	+	fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */
1124	+	skyclearness=skyclearsup-0.001;
1125	+	}
1126	+	if (skybrightness<skybriginf){
1127	+	/* if (suppress_warnings==0)
1128	+	fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */
1129	+	skybrightness=skybriginf;
1130	+	}
1131	+	if (skybrightness>skybrigsup){
1132	+	/* if (suppress_warnings==0)
1133	+	fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */
1134	+	skybrightness=skybrigsup;
1135	+	}
1136	+
1137	+	return; }
1138		else return;
1139		}
1140
1141
1142	<	/* likelihood of the direct and diffuse components */
1142	>
1143	>
1144	>
1145	>	/* validity of the direct and diffuse components */
1146		void    check_illuminances()
1147		{
1148	<	if (!( (directilluminance>=0) && (directilluminance<=solar_constant_l*1000) && (diffusilluminance>0) ))
1149	<	{
1150	<	fprintf(stderr,"direct or diffuse illuminances out of range\n");
1151	<	exit(1);
1148	>	if (directilluminance < 0) {
1149	>	if(suppress_warnings==0)
1150	>	{ fprintf(stderr,"Warning: direct illuminance < 0. Using 0.0\n"); }
1151	>	directilluminance = 0.0;
1152		}
1153	<	return;
1153	>	if (diffuseilluminance < 0) {
1154	>	if(suppress_warnings==0)
1155	>	{ fprintf(stderr,"Warning: diffuse illuminance < 0. Using 0.0\n"); }
1156	>	diffuseilluminance = 0.0;
1157	>	}
1158	>
1159	>	if (directilluminance+diffuseilluminance==0 && altitude > 0) {
1160	>	if(suppress_warnings==0)
1161	>	{ fprintf(stderr,"Warning: zero illuminance at sun altitude > 0, printing error sky\n"); }
1162	>	print_error_sky();
1163	>	exit(0);
1164	>	}
1165	>
1166	>	if (directilluminance > solar_constant_l) {
1167	>	if(suppress_warnings==0)
1168	>	{ fprintf(stderr,"Warning: direct illuminance exceeds solar constant\n"); }
1169	>	print_error_sky();
1170	>	exit(0);
1171	>	}
1172		}
1173
1174
1175		void    check_irradiances()
1176		{
1177	<	if (!( (directirradiance>=0) && (directirradiance<=solar_constant_e) && (diffusirradiance>0) ))
1178	<	{
1179	<	fprintf(stderr,"direct or diffuse irradiances out of range\n");
1180	<	exit(1);
1181	<	}
1182	<	return;
1177	>	if (directirradiance < 0) {
1178	>	if(suppress_warnings==0)
1179	>	{ fprintf(stderr,"Warning: direct irradiance < 0. Using 0.0\n"); }
1180	>	directirradiance = 0.0;
1181	>	}
1182	>	if (diffuseirradiance < 0) {
1183	>	if(suppress_warnings==0)
1184	>	{ fprintf(stderr,"Warning: diffuse irradiance < 0. Using 0.0\n"); }
1185	>	diffuseirradiance = 0.0;
1186	>	}
1187	>
1188	>	if (directirradiance+diffuseirradiance==0 && altitude > 0) {
1189	>	if(suppress_warnings==0)
1190	>	{ fprintf(stderr,"Warning: zero irradiance at sun altitude > 0, printing error sky\n"); }
1191	>	print_error_sky();
1192	>	exit(0);
1193	>	}
1194	>
1195	>	if (directirradiance > solar_constant_e) {
1196	>	if(suppress_warnings==0)
1197	>	{ fprintf(stderr,"Warning: direct irradiance exceeds solar constant\n"); }
1198	>	print_error_sky();
1199	>	exit(0);
1200	>	}
1201		}
1202
1203
#	Line 914 | Line 1207 | double sky_brightness()
1207		{
1208		double value;
1209
1210	<	value = diffusirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
1210	>	value = diffuseirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
1211
1212		return(value);
1213		}
#	Line 923 | Line 1216 | return(value);
1216		/* Perez sky's clearness */
1217		double sky_clearness()
1218		{
1219	<	double value;
1219	>	double value;
1220
1221	<	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) ;
1221	>	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) ;
1222
1223	<	return(value);
1223	>	return(value);
1224		}
1225
1226
1227
1228		/* diffus horizontal irradiance from Perez sky's brightness */
1229	<	double diffus_irradiance_from_sky_brightness()
1229	>	double diffuse_irradiance_from_sky_brightness()
1230		{
1231		double value;
1232
#	Line 948 | Line 1241 | double direct_irradiance_from_sky_clearness()
1241		{
1242		double value;
1243
1244	<	value = diffus_irradiance_from_sky_brightness();
1244	>	value = diffuse_irradiance_from_sky_brightness();
1245		value = value * ( (skyclearness-1) * (1+1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) );
1246
1247		return(value);
1248		}
1249
1250
1251	<	void illu_to_irra_index(void)
1251	>
1252	>
1253	>	void illu_to_irra_index()
1254		{
1255	<	double  test1=0.1, test2=0.1;
1255	>	double  test1=0.1, test2=0.1, d_eff;
1256		int     counter=0;
1257
1258	<	diffusirradiance = diffusilluminance*solar_constant_e/(solar_constant_l*1000);
1259	<	directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000);
1258	>	diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l);
1259	>	directirradiance = directilluminance*solar_constant_e/(solar_constant_l);
1260		skyclearness =  sky_clearness();
1261		skybrightness = sky_brightness();
1262	<	if (skyclearness>12) skyclearness=12;
968	<	if (skybrightness<0.05) skybrightness=0.01;
969	<
970	<
971	<	while ( ((fabs(diffusirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
972	<	|| skyclearness>skyclearinf || skyclearness<skyclearsup
973	<	|| skybrightness>skybriginf || skybrightness<skybrigsup )
974	<	&& !(counter==5) )
975	<	{
976	<	/*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
1262	>	check_parametrization();
1263
1264	<	test1=diffusirradiance;
1264	>
1265	>	while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
1266	>	|| (!(skyclearness<skyclearinf || skyclearness>skyclearsup))
1267	>	|| (!(skybrightness<skybriginf || skybrightness>skybrigsup)) )
1268	>	&& !(counter==9) )
1269	>	{
1270	>
1271	>	test1=diffuseirradiance;
1272		test2=directirradiance;
1273		counter++;
1274
1275	<	diffusirradiance = diffusilluminance/glob_h_diffuse_effi_PEREZ();
1276	<	directirradiance = directilluminance/direct_n_effi_PEREZ();
984	<	/*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
1275	>	diffuseirradiance = diffuseilluminance/glob_h_diffuse_effi_PEREZ();
1276	>	d_eff = direct_n_effi_PEREZ();
1277
1278	+
1279	+	if (d_eff < 0.1)
1280	+	directirradiance = 0;
1281	+	else
1282	+	directirradiance = directilluminance/d_eff;
1283	+
1284		skybrightness = sky_brightness();
1285		skyclearness =  sky_clearness();
1286	<	if (skyclearness>12) skyclearness=12;
1287	<	if (skybrightness<0.05) skybrightness=0.01;
990	<
991	<	/*fprintf(stderr, "%lf\t %lf\n", skybrightness, skyclearness);*/
992	<
1286	>	check_parametrization();
1287	>
1288		}
1289
1290
1291		return;
1292		}
1293
1294	<
1000	<	int lect_coeff_perez(char *filename,float **coeff_perez)
1294	>	static int get_numlin(float epsilon)
1295		{
1296	<	FILE *fcoeff_perez;
1297	<	float temp;
1298	<	int i,j;
1299	<
1300	<	if ((fcoeff_perez = frlibopen(filename)) == NULL)
1301	<	{
1302	<	fprintf(stderr,"file %s cannot be opened\n", filename);
1303	<	return 1; /* il y a un probleme de fichier */
1304	<	}
1305	<	else
1306	<	{
1307	<	/*printf("file %s  open\n", filename);*/
1308	<	}
1309	<
1310	<	skip_comments(fcoeff_perez);
1017	<
1018	<	for (i=0;i<8;i++)
1019	<	for (j=0;j<20;j++)
1020	<	{
1021	<	fscanf(fcoeff_perez,"%f",&temp);
1022	<	*(*coeff_perez+i*20+j) = temp;
1023	<	}
1024	<	fclose(fcoeff_perez);
1025	<
1026	<	return 0; /* tout est OK */
1296	>	if (epsilon < 1.065)
1297	>	return 0;
1298	>	else if (epsilon < 1.230)
1299	>	return 1;
1300	>	else if (epsilon < 1.500)
1301	>	return 2;
1302	>	else if (epsilon < 1.950)
1303	>	return 3;
1304	>	else if (epsilon < 2.800)
1305	>	return 4;
1306	>	else if (epsilon < 4.500)
1307	>	return 5;
1308	>	else if (epsilon < 6.200)
1309	>	return 6;
1310	>	return 7;
1311		}
1312
1029	–
1030	–
1313		/* sky luminance perez model */
1314	<	double calc_rel_lum_perez(double dzeta,double gamma,double Z,
1033	<	double epsilon,double Delta,float *coeff_perez)
1314	>	double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[])
1315		{
1316	+
1317		float x[5][4];
1318		int i,j,num_lin;
1319		double c_perez[5];
1320
1321		if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
1322		{
1323	<	fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n");
1323	>	fprintf(stderr,"Error: epsilon out of range in function calc_rel_lum_perez!\n");
1324		exit(1);
1325		}
1326
#	Line 1047 | Line 1329 | double epsilon,double Delta,float *coeff_perez)
1329		{
1330		if ( Delta < 0.2 ) Delta = 0.2;
1331		}
1332	<
1333	<	if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0;
1334	<	if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
1335	<	if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2;
1054	<	if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3;
1055	<	if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4;
1056	<	if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5;
1057	<	if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6;
1058	<	if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7;
1059	<
1332	>
1333	>
1334	>	num_lin = get_numlin(epsilon);
1335	>
1336		for (i=0;i<5;i++)
1337		for (j=0;j<4;j++)
1338		{
1339		x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j);
1340	<	/* printf("x %d %d vaut %f\n",i,j,x[i][j]); */
1340	>	/* fprintf(stderr,"x %d %d vaut %f\n",i,j,x[i][j]); */
1341		}
1342
1343
#	Line 1088 | Line 1364 | double epsilon,double Delta,float *coeff_perez)
1364
1365
1366		/* coefficients for the sky luminance perez model */
1367	<	void coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez)
1367	>	void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[])
1368		{
1369		float x[5][4];
1370		int i,j,num_lin;
1371
1372		if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
1373		{
1374	<	fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n");
1374	>	fprintf(stderr,"Error: epsilon out of range in function coeff_lum_perez!\n");
1375		exit(1);
1376		}
1377
#	Line 1104 | Line 1380 | void coeff_lum_perez(double Z, double epsilon, double
1380		{
1381		if ( Delta < 0.2 ) Delta = 0.2;
1382		}
1383	+
1384	+
1385	+	num_lin = get_numlin(epsilon);
1386
1387	<	if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0;
1109	<	if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
1110	<	if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2;
1111	<	if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3;
1112	<	if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4;
1113	<	if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5;
1114	<	if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6;
1115	<	if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7;
1387	>	/*fprintf(stderr,"numlin %d\n", num_lin);*/
1388
1389		for (i=0;i<5;i++)
1390		for (j=0;j<4;j++)
#	Line 1144 | Line 1416 | void coeff_lum_perez(double Z, double epsilon, double
1416		}
1417
1418
1419	+
1420		/* degrees into radians */
1421		double radians(double degres)
1422		{
1423		return degres*M_PI/180.0;
1424		}
1425
1426	+
1427		/* radian into degrees */
1428		double degres(double radians)
1429		{
1430		return radians/M_PI*180.0;
1431		}
1432
1433	+
1434		/* calculation of the angles dzeta and gamma */
1435		void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z)
1436		{
#	Line 1165 | Line 1440 | void theta_phi_to_dzeta_gamma(double theta,double phi,
1440		else if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1.1 )
1441		{
1442		printf("error in calculation of gamma (angle between point and sun");
1443	<	exit(3);
1443	>	exit(1);
1444		}
1445		else
1446		*gamma = acos(cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi));
1447		}
1448
1449
1175	–	/********************************************************************************/
1176	–	/*      Fonction: theta_ordered                                                 */
1177	–	/*                                                                              */
1178	–	/*      In: char *filename                                                      */
1179	–	/*                                                                              */
1180	–	/*      Out: float *                                                            */
1181	–	/*                                                                              */
1182	–	/*      Update: 29/08/93                                                        */
1183	–	/*                                                                              */
1184	–	/*      Rem: theta en degres                                                    */
1185	–	/*                                                                              */
1186	–	/*      But: fournit les valeurs de theta du fichier d'entree a la memoire      */
1187	–	/*                                                                              */
1188	–	/********************************************************************************/
1189	–	float *theta_ordered(char *filename)
1190	–	{
1191	–	int i;
1192	–	float buffer,*ptr;
1193	–	FILE *file_in;
1450
1195	–	if ( (file_in = frlibopen(filename)) == NULL )
1196	–	{
1197	–	fprintf(stderr,"Cannot open file %s in function theta_ordered\n",filename);
1198	–	exit(1);
1199	–	}
1200	–
1201	–	skip_comments(file_in);
1202	–
1203	–	if ( (ptr = malloc(145*sizeof(float))) == NULL )
1204	–	{
1205	–	fprintf(stderr,"Out of memory in function theta_ordered\n");
1206	–	exit(1);
1207	–	}
1208	–
1209	–	for (i=0;i<145;i++)
1210	–	{
1211	–	fscanf(file_in,"%f",&buffer);
1212	–	*(ptr+i) = buffer;
1213	–	fscanf(file_in,"%f",&buffer);
1214	–	}
1215	–
1216	–	fclose(file_in);
1217	–	return ptr;
1218	–	}
1219	–
1220	–
1221	–	/********************************************************************************/
1222	–	/*      Fonction: phi_ordered                                                   */
1223	–	/*                                                                              */
1224	–	/*      In: char *filename                                                      */
1225	–	/*                                                                              */
1226	–	/*      Out: float *                                                            */
1227	–	/*                                                                              */
1228	–	/*      Update: 29/08/93                                                        */
1229	–	/*                                                                              */
1230	–	/*      Rem: valeurs de Phi en DEGRES                                           */
1231	–	/*                                                                              */
1232	–	/*      But: mettre les angles contenus dans le fichier d'entree dans la memoire */
1233	–	/*                                                                              */
1234	–	/********************************************************************************/
1235	–	float *phi_ordered(char *filename)
1236	–	{
1237	–	int i;
1238	–	float buffer,*ptr;
1239	–	FILE *file_in;
1240	–
1241	–	if ( (file_in = frlibopen(filename)) == NULL )
1242	–	{
1243	–	fprintf(stderr,"Cannot open file %s in function phi_ordered\n",filename);
1244	–	exit(1);
1245	–	}
1246	–
1247	–	skip_comments(file_in);
1248	–
1249	–	if ( (ptr = malloc(145*sizeof(float))) == NULL )
1250	–	{
1251	–	fprintf(stderr,"Out of memory in function phi_ordered");
1252	–	exit(1);
1253	–	}
1254	–
1255	–	for (i=0;i<145;i++)
1256	–	{
1257	–	fscanf(file_in,"%f",&buffer);
1258	–	fscanf(file_in,"%f",&buffer);
1259	–	*(ptr+i) = buffer;
1260	–	}
1261	–
1262	–	fclose(file_in);
1263	–	return ptr;
1264	–	}
1265	–
1266	–
1267	–	/********************************************************************************/
1268	–	/*      Fonction: integ_lv                                                      */
1269	–	/*                                                                              */
1270	–	/*      In: float *lv,*theta                                                    */
1271	–	/*          int sun_pos                                                         */
1272	–	/*                                                                              */
1273	–	/*      Out: double                                                             */
1274	–	/*                                                                              */
1275	–	/*      Update: 29/08/93                                                        */
1276	–	/*                                                                              */
1277	–	/*      Rem:                                                                    */
1278	–	/*                                                                              */
1279	–	/*      But: calcul l'integrale de luminance relative sans la dir. du soleil    */
1280	–	/*                                                                              */
1281	–	/********************************************************************************/
1451		double integ_lv(float *lv,float *theta)
1452		{
1453		int i;
1454		double buffer=0.0;
1455	<
1455	>
1456		for (i=0;i<145;i++)
1457	+	{
1458		buffer += (*(lv+i))*cos(radians(*(theta+i)));
1459	<
1459	>	}
1460	>
1461		return buffer*2*M_PI/144;
1291	–
1462		}
1463
1464
1465
1296	–
1297	–
1298	–
1466		/* enter day number(double), return E0 = square(R0/R):  eccentricity correction factor  */
1467
1468		double get_eccentricity()
#	Line 1308 | Line 1475 | double get_eccentricity()
1475		0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle);
1476
1477		return (E0);
1311	–
1478		}
1479
1480
#	Line 1316 | Line 1482 | double get_eccentricity()
1482		double  air_mass()
1483		{
1484		double  m;
1319	–
1485		if (sunzenith>90)
1486		{
1487	<	fprintf(stderr, "solar zenith angle larger than 90� in fuction air_mass():\n the models used are not more valid\n");
1488	<	exit(1);
1487	>	if(suppress_warnings==0)
1488	>	{ fprintf(stderr, "Warning: air mass has reached the maximal value\n"); }
1489	>	sunzenith=90;
1490		}
1325	–
1491		m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1492		return(m);
1493		}
1329	–
1330	–
1331	–	double get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut)
1332	–
1333	–	{
1334	–
1335	–	double angle;
1336	–
1337	–
1338	–	if (sun_zenith == 0)
1339	–	puts("WARNING: zenith_angle = 0 in function get_angle_sun_vert_plan");
1340	–
1341	–	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) );
1342	–	angle = angle*180/M_PI;
1343	–	return(angle);
1344	–	}
1345	–
1346	–
1347	–
1348	–
1349	–
1494
1495
1496

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)