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 |
+ |
* version 2.6 (2021/01/29): dew point dependency added according to Perez publication 1990 (W -> atm_preci_water=exp(0.07*Td-0.075) ). by J. Wienold, EPFL |
9 |
|
*/ |
10 |
|
|
11 |
< |
|
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 |
< |
|
11 |
> |
#define _USE_MATH_DEFINES |
12 |
|
#include <stdio.h> |
13 |
|
#include <string.h> |
14 |
|
#include <math.h> |
15 |
|
#include <stdlib.h> |
16 |
|
|
25 |
– |
#include "rtio.h" |
26 |
– |
#include "fvect.h" |
17 |
|
#include "color.h" |
18 |
+ |
#include "sun.h" |
19 |
|
#include "paths.h" |
20 |
|
|
21 |
< |
extern int jdate(int month, int day); |
31 |
< |
extern double stadj(int jd); |
32 |
< |
extern double sdec(int jd); |
33 |
< |
extern double salt(double sd, double st); |
34 |
< |
extern double sazi(double sd, double st); |
21 |
> |
#define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]) |
22 |
|
|
23 |
|
double normsc(); |
24 |
|
|
25 |
< |
#define DATFILE "coeff_perez.dat" |
25 |
> |
/*static char *rcsid="$Header$";*/ |
26 |
|
|
27 |
+ |
float coeff_perez[] = { |
28 |
+ |
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, |
29 |
+ |
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, |
30 |
+ |
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, |
31 |
+ |
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, |
32 |
+ |
-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, |
33 |
+ |
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, |
34 |
+ |
-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, |
35 |
+ |
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, |
36 |
+ |
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, |
37 |
+ |
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, |
38 |
+ |
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, |
39 |
+ |
-14.5000,-46.1148,55.3750,-7.2312,0.4050,13.3500,0.6234,1.5000,-0.6426,1.8564,0.5636}; |
40 |
|
|
41 |
|
|
42 |
< |
/* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */ |
42 |
> |
float defangle_theta[] = { |
43 |
> |
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, |
44 |
> |
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, |
45 |
> |
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, |
46 |
> |
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, |
47 |
> |
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, |
48 |
> |
24, 24, 24, 24, 24, 24, 24, 24, 12, 12, 12, 12, 12, 12, 0}; |
49 |
> |
|
50 |
> |
float defangle_phi[] = { |
51 |
> |
0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264, |
52 |
> |
276, 288, 300, 312, 324, 336, 348, 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, |
53 |
> |
192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 0, 15, 30, 45, 60, 75, 90, 105, |
54 |
> |
120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 15, 30, 45, 60, 75, |
55 |
> |
90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 20, 40, 60, |
56 |
> |
80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 0, 30, 60, 90, 120, 150, 180, 210, |
57 |
> |
240, 270, 300, 330, 0, 60, 120, 180, 240, 300, 0}; |
58 |
> |
/* default values for Berlin */ |
59 |
> |
float locus[] = { |
60 |
> |
-4.843e9,2.5568e6,0.24282e3,0.23258,-4.843e9,2.5568e6,0.24282e3,0.23258,-1.2848,1.7519,-0.093786}; |
61 |
> |
|
62 |
> |
|
63 |
> |
|
64 |
> |
/* Perez sky parametrization: epsilon and delta calculations from the direct and diffuse irradiances */ |
65 |
|
double sky_brightness(); |
66 |
|
double sky_clearness(); |
67 |
|
|
68 |
|
/* calculation of the direct and diffuse components from the Perez parametrization */ |
69 |
< |
double diffus_irradiance_from_sky_brightness(); |
69 |
> |
double diffuse_irradiance_from_sky_brightness(); |
70 |
|
double direct_irradiance_from_sky_clearness(); |
71 |
|
|
72 |
+ |
/* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : */ |
73 |
+ |
/* input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */ |
74 |
|
|
51 |
– |
/* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */ |
75 |
|
double glob_h_effi_PEREZ(); |
76 |
|
double glob_h_diffuse_effi_PEREZ(); |
77 |
|
double direct_n_effi_PEREZ(); |
78 |
+ |
|
79 |
|
/*likelihood check of the epsilon, delta, direct and diffuse components*/ |
80 |
|
void check_parametrization(); |
81 |
|
void check_irradiances(); |
82 |
|
void check_illuminances(); |
83 |
|
void illu_to_irra_index(); |
84 |
+ |
void print_error_sky(); |
85 |
|
|
86 |
< |
|
87 |
< |
/* Perez sky luminance model */ |
63 |
< |
int lect_coeff_perez(char *filename,float **coeff_perez); |
64 |
< |
double calc_rel_lum_perez(double dzeta,double gamma,double Z, |
65 |
< |
double epsilon,double Delta,float *coeff_perez); |
66 |
< |
/* coefficients for the sky luminance perez model */ |
67 |
< |
void coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez); |
86 |
> |
double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]); |
87 |
> |
void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]); |
88 |
|
double radians(double degres); |
89 |
|
double degres(double radians); |
90 |
|
void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z); |
91 |
|
double integ_lv(float *lv,float *theta); |
72 |
– |
float *theta_ordered(char *filename); |
73 |
– |
float *phi_ordered(char *filename); |
92 |
|
|
93 |
+ |
void printdefaults(); |
94 |
+ |
void check_sun_position(); |
95 |
+ |
void computesky(); |
96 |
+ |
void printhead(int ac, char** av); |
97 |
+ |
void usage_error(char* msg); |
98 |
+ |
void printsky(); |
99 |
|
|
100 |
+ |
FILE * frlibopen(char* fname); |
101 |
|
|
102 |
|
/* astronomy and geometry*/ |
103 |
|
double get_eccentricity(); |
104 |
|
double air_mass(); |
80 |
– |
double get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut); |
105 |
|
|
106 |
+ |
double solar_sunset(int month, int day); |
107 |
+ |
double solar_sunrise(int month, int day); |
108 |
|
|
83 |
– |
/* date*/ |
84 |
– |
int jdate(int month, int day); |
85 |
– |
|
86 |
– |
|
87 |
– |
|
88 |
– |
|
89 |
– |
|
90 |
– |
/* sun calculation constants */ |
91 |
– |
extern double s_latitude; |
92 |
– |
extern double s_longitude; |
93 |
– |
extern double s_meridian; |
94 |
– |
|
109 |
|
const double AU = 149597890E3; |
110 |
|
const double solar_constant_e = 1367; /* solar constant W/m^2 */ |
111 |
< |
const double solar_constant_l = 127.5; /* solar constant klux */ |
111 |
> |
const double solar_constant_l = 127500; /* solar constant lux */ |
112 |
|
|
113 |
|
const double half_sun_angle = 0.2665; |
114 |
|
const double half_direct_angle = 2.85; |
115 |
|
|
116 |
< |
const double skyclearinf = 1.000; /* limitations for the variation of the Perez parameters */ |
117 |
< |
const double skyclearsup = 12.1; |
116 |
> |
const double skyclearinf = 1.0; /* limitations for the variation of the Perez parameters */ |
117 |
> |
const double skyclearsup = 12.01; |
118 |
|
const double skybriginf = 0.01; |
119 |
|
const double skybrigsup = 0.6; |
120 |
|
|
121 |
|
|
122 |
|
|
123 |
|
/* required values */ |
124 |
+ |
int year = 0; /* year (optional) */ |
125 |
|
int month, day; /* date */ |
126 |
|
double hour; /* time */ |
127 |
|
int tsolar; /* 0=standard, 1=solar */ |
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, Td=10.97353115; |
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; |
150 |
< |
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); |
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++) |
201 |
|
if (argv[i][0] == '-' || argv[i][0] == '+') |
202 |
|
switch (argv[i][1]) { |
203 |
+ |
case 'd': |
204 |
+ |
Td = atof(argv[++i]); |
205 |
+ |
if (Td < -40 || Td > 40) { |
206 |
+ |
Td=10.97353115; } |
207 |
+ |
break; |
208 |
|
case 's': |
209 |
|
cloudy = 0; |
210 |
|
dosun = argv[i][0] == '+'; |
211 |
|
break; |
212 |
< |
case 'r': |
212 |
> |
case 'y': |
213 |
> |
year = atoi(argv[++i]); |
214 |
> |
break; |
215 |
|
case 'R': |
216 |
|
u_solar = argv[i][1] == 'R' ? -1 : 1; |
217 |
|
solarbr = atof(argv[++i]); |
220 |
|
cloudy = argv[i][0] == '+' ? 2 : 1; |
221 |
|
dosun = 0; |
222 |
|
break; |
223 |
+ |
case 'C': |
224 |
+ |
if (argv[i][2] == 'I' && argv[i][3] == 'E' ) { |
225 |
+ |
locus[0] = -4.607e9; |
226 |
+ |
locus[1] = 2.9678e6; |
227 |
+ |
locus[2] = 0.09911e3; |
228 |
+ |
locus[3] = 0.244063; |
229 |
+ |
locus[4] = -2.0064e9; |
230 |
+ |
locus[5] = 1.9018e6; |
231 |
+ |
locus[6] = 0.24748e3; |
232 |
+ |
locus[7] = 0.23704; |
233 |
+ |
locus[8] = -3.0; |
234 |
+ |
locus[9] = 2.87; |
235 |
+ |
locus[10] = -0.275; |
236 |
+ |
}else{ color_output = 1; |
237 |
+ |
} |
238 |
+ |
break; |
239 |
+ |
case 'l': |
240 |
+ |
locus[0] = atof(argv[++i]); |
241 |
+ |
locus[1] = atof(argv[++i]); |
242 |
+ |
locus[2] = atof(argv[++i]); |
243 |
+ |
locus[3] = atof(argv[++i]); |
244 |
+ |
locus[4] = locus[0]; |
245 |
+ |
locus[5] = locus[1]; |
246 |
+ |
locus[6] = locus[2]; |
247 |
+ |
locus[7] = locus[3]; |
248 |
+ |
locus[8] = atof(argv[++i]); |
249 |
+ |
locus[9] = atof(argv[++i]); |
250 |
+ |
locus[10] = atof(argv[++i]); |
251 |
+ |
break; |
252 |
+ |
|
253 |
|
case 't': |
254 |
|
betaturbidity = atof(argv[++i]); |
255 |
|
break; |
256 |
+ |
case 'w': |
257 |
+ |
suppress_warnings = 1; |
258 |
+ |
break; |
259 |
|
case 'b': |
260 |
|
zenithbr = atof(argv[++i]); |
261 |
|
break; |
271 |
|
case 'm': |
272 |
|
s_meridian = atof(argv[++i]) * (M_PI/180); |
273 |
|
break; |
211 |
– |
|
274 |
|
|
275 |
|
case 'O': |
276 |
< |
output = atof(argv[++i]); /*define the unit of the output of the program : |
277 |
< |
sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm) |
278 |
< |
default is set to 0*/ |
276 |
> |
output = atof(argv[++i]); /*define the unit of the output of the program: |
277 |
> |
sky and sun luminance/radiance |
278 |
> |
(0==W visible, 1==W solar radiation, 2==lm) */ |
279 |
|
break; |
280 |
|
|
281 |
|
case 'P': |
287 |
|
case 'W': /* direct normal Irradiance [W/m^2] */ |
288 |
|
input = 1; /* diffuse horizontal Irrad. [W/m^2] */ |
289 |
|
directirradiance = atof(argv[++i]); |
290 |
< |
diffusirradiance = atof(argv[++i]); |
290 |
> |
diffuseirradiance = atof(argv[++i]); |
291 |
|
break; |
292 |
|
|
293 |
|
case 'L': /* direct normal Illuminance [Lux] */ |
294 |
|
input = 2; /* diffuse horizontal Ill. [Lux] */ |
295 |
|
directilluminance = atof(argv[++i]); |
296 |
< |
diffusilluminance = atof(argv[++i]); |
296 |
> |
diffuseilluminance = atof(argv[++i]); |
297 |
|
break; |
298 |
|
|
299 |
|
case 'G': /* direct horizontal Irradiance [W/m^2] */ |
300 |
|
input = 3; /* diffuse horizontal Irrad. [W/m^2] */ |
301 |
|
directirradiance = atof(argv[++i]); |
302 |
< |
diffusirradiance = atof(argv[++i]); |
302 |
> |
diffuseirradiance = atof(argv[++i]); |
303 |
|
break; |
242 |
– |
|
304 |
|
|
305 |
+ |
case 'E': /* Erbs model based on the */ |
306 |
+ |
input = 4; /* global-horizontal irradiance [W/m^2] */ |
307 |
+ |
globalirradiance = atof(argv[++i]); |
308 |
+ |
break; |
309 |
+ |
|
310 |
+ |
case 'i': |
311 |
+ |
timeinterval = atof(argv[++i]); |
312 |
+ |
break; |
313 |
+ |
|
314 |
+ |
|
315 |
|
default: |
316 |
|
sprintf(errmsg, "unknown option: %s", argv[i]); |
317 |
< |
userror(errmsg); |
317 |
> |
usage_error(errmsg); |
318 |
|
} |
319 |
|
else |
320 |
< |
userror("bad option"); |
320 |
> |
usage_error("bad option"); |
321 |
|
|
322 |
< |
if (fabs(s_meridian-s_longitude) > 30*M_PI/180) |
323 |
< |
fprintf(stderr, |
253 |
< |
"%s: warning: %.1f hours btwn. standard meridian and longitude\n", |
322 |
> |
if (month && !tsolar && fabs(s_meridian-s_longitude) > 45*M_PI/180) |
323 |
> |
fprintf(stderr,"%s: warning: %.1f hours btwn. standard meridian and longitude\n", |
324 |
|
progname, (s_longitude-s_meridian)*12/M_PI); |
325 |
|
|
326 |
|
|
327 |
< |
/* allocation dynamique de memoire pour les pointeurs */ |
328 |
< |
if ( (c_perez = malloc(5*sizeof(double))) == NULL ) |
329 |
< |
{ |
260 |
< |
fprintf(stderr,"Out of memory error in function main !"); |
261 |
< |
exit(1); |
262 |
< |
} |
327 |
> |
/* dynamic memory allocation for the pointers */ |
328 |
> |
if ( (c_perez = calloc(5, sizeof(double))) == NULL ) |
329 |
> |
{ fprintf(stderr,"Out of memory error in function main"); return 1; } |
330 |
|
|
331 |
< |
|
331 |
> |
|
332 |
> |
atm_preci_water=exp(0.07*Td-0.075); |
333 |
|
printhead(argc, argv); |
266 |
– |
|
334 |
|
computesky(); |
335 |
|
printsky(); |
336 |
+ |
return 0; |
337 |
|
|
270 |
– |
exit(0); |
338 |
|
} |
339 |
|
|
340 |
|
|
341 |
< |
computesky() /* compute sky parameters */ |
341 |
> |
|
342 |
> |
|
343 |
> |
|
344 |
> |
void computesky() |
345 |
|
{ |
346 |
|
|
347 |
< |
/* new variables */ |
348 |
< |
int j, i; |
349 |
< |
float *lv_mod; /* 145 luminance values*/ |
350 |
< |
/* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */ |
281 |
< |
float *theta_o, *phi_o, *coeff_perez; |
347 |
> |
int j; |
348 |
> |
|
349 |
> |
float *lv_mod; /* 145 luminance values */ |
350 |
> |
float *theta_o, *phi_o; |
351 |
|
double dzeta, gamma; |
283 |
– |
double diffusion; |
352 |
|
double normfactor; |
353 |
+ |
double erbs_s0, erbs_kt; |
354 |
|
|
355 |
|
|
287 |
– |
|
356 |
|
/* compute solar direction */ |
357 |
< |
|
357 |
> |
|
358 |
|
if (month) { /* from date and time */ |
359 |
< |
int jd; |
292 |
< |
double sd, st; |
359 |
> |
double sd; |
360 |
|
|
361 |
< |
jd = jdate(month, day); /* Julian date */ |
362 |
< |
sd = sdec(jd); /* solar declination */ |
363 |
< |
if (tsolar) /* solar time */ |
364 |
< |
st = hour; |
361 |
> |
st = hour; |
362 |
> |
if (year) { /* Michalsky algorithm? */ |
363 |
> |
double mjd = mjdate(year, month, day, hour); |
364 |
> |
if (tsolar) |
365 |
> |
sd = msdec(mjd, NULL); |
366 |
> |
else |
367 |
> |
sd = msdec(mjd, &st); |
368 |
> |
} else { |
369 |
> |
int jd = jdate(month, day); /* Julian date */ |
370 |
> |
sd = sdec(jd); /* solar declination */ |
371 |
> |
if (!tsolar) /* get solar time? */ |
372 |
> |
st = hour + stadj(jd); |
373 |
> |
} |
374 |
> |
|
375 |
> |
if(timeinterval) { |
376 |
> |
|
377 |
> |
if(timeinterval<0) { |
378 |
> |
fprintf(stderr, "time interval negative\n"); |
379 |
> |
exit(1); |
380 |
> |
} |
381 |
> |
|
382 |
> |
if(fabs(solar_sunrise(month,day)-st)<=timeinterval/120) { |
383 |
> |
st= (st+timeinterval/120+solar_sunrise(month,day))/2; |
384 |
> |
if(suppress_warnings==0) |
385 |
> |
{ fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); } |
386 |
> |
} |
387 |
> |
|
388 |
> |
if(fabs(solar_sunset(month,day)-st)<timeinterval/120) { |
389 |
> |
st= (st-timeinterval/120+solar_sunset(month,day))/2; |
390 |
> |
if(suppress_warnings==0) |
391 |
> |
{ fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); } |
392 |
> |
} |
393 |
> |
|
394 |
> |
if((st<solar_sunrise(month,day)-timeinterval/120) || (st>solar_sunset(month,day)+timeinterval/120)) { |
395 |
> |
if(suppress_warnings==0) |
396 |
> |
{ fprintf(stderr, "Warning: sun position too low, printing error sky at %d %d %.3f\n",month,day,hour); } |
397 |
> |
altitude = salt(sd, st); |
398 |
> |
azimuth = sazi(sd, st); |
399 |
> |
print_error_sky(); |
400 |
> |
exit(0); |
401 |
> |
} |
402 |
> |
} |
403 |
|
else |
404 |
< |
st = hour + stadj(jd); |
404 |
> |
|
405 |
> |
if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) { |
406 |
> |
if(suppress_warnings==0) |
407 |
> |
{ fprintf(stderr, "Warning: sun altitude below zero at time step %i %i %.2f, printing error sky\n",month,day,hour); } |
408 |
> |
altitude = salt(sd, st); |
409 |
> |
azimuth = sazi(sd, st); |
410 |
> |
print_error_sky(); |
411 |
> |
exit(0); |
412 |
> |
} |
413 |
> |
|
414 |
|
altitude = salt(sd, st); |
415 |
|
azimuth = sazi(sd, st); |
416 |
|
|
417 |
|
daynumber = (double)jdate(month, day); |
418 |
< |
|
418 |
> |
|
419 |
|
} |
420 |
+ |
|
421 |
+ |
|
422 |
+ |
|
423 |
+ |
|
424 |
+ |
|
425 |
|
if (!cloudy && altitude > 87.*M_PI/180.) { |
426 |
< |
fprintf(stderr, |
426 |
> |
|
427 |
> |
if (suppress_warnings==0) { |
428 |
> |
fprintf(stderr, |
429 |
|
"%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n", |
430 |
|
progname); |
431 |
< |
printf( |
311 |
< |
"# warning - sun too close to zenith, reducing altitude to 87 degrees\n"); |
431 |
> |
} |
432 |
|
altitude = 87.*M_PI/180.; |
433 |
|
} |
434 |
+ |
|
435 |
+ |
|
436 |
+ |
|
437 |
|
sundir[0] = -sin(azimuth)*cos(altitude); |
438 |
|
sundir[1] = -cos(azimuth)*cos(altitude); |
439 |
|
sundir[2] = sin(altitude); |
441 |
|
|
442 |
|
/* calculation for the new functions */ |
443 |
|
sunzenith = 90 - altitude*180/M_PI; |
444 |
< |
|
322 |
< |
|
444 |
> |
|
445 |
|
|
446 |
< |
/* compute the inputs for the calculation of the light distribution over the sky*/ |
447 |
< |
if (input==0) |
446 |
> |
/* compute the inputs for the calculation of the light distribution over the sky*/ |
447 |
> |
if (input==0) /* P */ |
448 |
|
{ |
449 |
|
check_parametrization(); |
450 |
< |
diffusirradiance = diffus_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/ |
450 |
> |
diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/ |
451 |
|
directirradiance = direct_irradiance_from_sky_clearness(); |
452 |
|
check_irradiances(); |
453 |
|
|
454 |
|
if (output==0 || output==2) |
455 |
|
{ |
456 |
< |
diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
456 |
> |
diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
457 |
|
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
458 |
|
check_illuminances(); |
459 |
|
} |
460 |
|
} |
461 |
|
|
462 |
|
|
463 |
< |
else if (input==1) |
463 |
> |
else if (input==1) /* W */ |
464 |
|
{ |
465 |
|
check_irradiances(); |
466 |
|
skybrightness = sky_brightness(); |
467 |
|
skyclearness = sky_clearness(); |
468 |
+ |
|
469 |
|
check_parametrization(); |
470 |
< |
|
470 |
> |
|
471 |
|
if (output==0 || output==2) |
472 |
|
{ |
473 |
< |
diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
473 |
> |
diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
474 |
|
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
475 |
|
check_illuminances(); |
476 |
|
} |
478 |
|
} |
479 |
|
|
480 |
|
|
481 |
< |
else if (input==2) |
481 |
> |
else if (input==2) /* L */ |
482 |
|
{ |
483 |
|
check_illuminances(); |
484 |
|
illu_to_irra_index(); |
486 |
|
} |
487 |
|
|
488 |
|
|
489 |
< |
else if (input==3) |
489 |
> |
else if (input==3) /* G */ |
490 |
|
{ |
491 |
|
if (altitude<=0) |
492 |
|
{ |
493 |
< |
fprintf(stderr, "solar zenith angle larger than 90� \n the models used are not more valid\n"); |
494 |
< |
exit(1); |
493 |
> |
if (suppress_warnings==0) |
494 |
> |
fprintf(stderr, "Warning: sun altitude < 0, proceed with irradiance values of zero\n"); |
495 |
> |
directirradiance = 0; |
496 |
> |
diffuseirradiance = 0; |
497 |
> |
} else { |
498 |
> |
|
499 |
> |
directirradiance=directirradiance/sin(altitude); |
500 |
|
} |
501 |
< |
|
374 |
< |
directirradiance=directirradiance/sin(altitude); |
501 |
> |
|
502 |
|
check_irradiances(); |
503 |
|
skybrightness = sky_brightness(); |
504 |
|
skyclearness = sky_clearness(); |
506 |
|
|
507 |
|
if (output==0 || output==2) |
508 |
|
{ |
509 |
< |
diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
509 |
> |
diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
510 |
|
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
511 |
|
check_illuminances(); |
512 |
|
} |
513 |
|
|
514 |
|
} |
515 |
|
|
389 |
– |
|
390 |
– |
else {fprintf(stderr,"error in giving the input arguments"); exit(1);} |
516 |
|
|
517 |
+ |
else if (input==4) /* E */ /* Implementation of the Erbs model. W.Sprenger (04/13) */ |
518 |
+ |
{ |
519 |
+ |
|
520 |
+ |
if (altitude<=0) |
521 |
+ |
{ |
522 |
+ |
if (suppress_warnings==0 && globalirradiance > 50) |
523 |
+ |
fprintf(stderr, "Warning: global irradiance higher than 50 W/m^2 while the sun altitude is lower than zero\n"); |
524 |
+ |
globalirradiance = 0; diffuseirradiance = 0; directirradiance = 0; |
525 |
+ |
|
526 |
+ |
} else { |
527 |
+ |
|
528 |
+ |
erbs_s0 = solar_constant_e*get_eccentricity()*sin(altitude); |
529 |
+ |
|
530 |
+ |
if (globalirradiance>erbs_s0) |
531 |
+ |
{ |
532 |
+ |
if (suppress_warnings==0) |
533 |
+ |
fprintf(stderr, "Warning: global irradiance is higher than the time-dependent solar constant s0\n"); |
534 |
+ |
globalirradiance=erbs_s0*0.999; |
535 |
+ |
} |
536 |
+ |
|
537 |
+ |
erbs_kt=globalirradiance/erbs_s0; |
538 |
+ |
|
539 |
+ |
if (erbs_kt<=0.22) diffuseirradiance=globalirradiance*(1-0.09*erbs_kt); |
540 |
+ |
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)); |
541 |
+ |
else if (erbs_kt<1) diffuseirradiance=globalirradiance*(0.165); |
542 |
+ |
|
543 |
+ |
directirradiance=globalirradiance-diffuseirradiance; |
544 |
+ |
|
545 |
+ |
printf("# erbs_s0, erbs_kt, irr_dir_h, irr_diff: %.3f %.3f %.3f %.3f\n", erbs_s0, erbs_kt, directirradiance, diffuseirradiance); |
546 |
+ |
printf("# WARNING: the -E option is only recommended for a rough estimation!\n"); |
547 |
+ |
|
548 |
+ |
directirradiance=directirradiance/sin(altitude); |
549 |
+ |
|
550 |
+ |
} |
551 |
+ |
|
552 |
+ |
check_irradiances(); |
553 |
+ |
skybrightness = sky_brightness(); |
554 |
+ |
skyclearness = sky_clearness(); |
555 |
+ |
check_parametrization(); |
556 |
|
|
557 |
+ |
if (output==0 || output==2) |
558 |
+ |
{ |
559 |
+ |
diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
560 |
+ |
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
561 |
+ |
check_illuminances(); |
562 |
+ |
} |
563 |
+ |
|
564 |
+ |
} |
565 |
+ |
|
566 |
+ |
|
567 |
+ |
|
568 |
|
|
569 |
< |
/* normalization factor for the relative sky luminance distribution, diffuse part*/ |
569 |
> |
else { fprintf(stderr,"error at the input arguments"); exit(1); } |
570 |
|
|
396 |
– |
/* allocation dynamique de memoire pour les pointeurs */ |
397 |
– |
if ( (coeff_perez = malloc(8*20*sizeof(float))) == NULL ) |
398 |
– |
{ |
399 |
– |
fprintf(stderr,"Out of memory error in function main !"); |
400 |
– |
exit(1); |
401 |
– |
} |
571 |
|
|
572 |
< |
/* read the coefficients for the Perez sky luminance model */ |
573 |
< |
if (lect_coeff_perez(DATFILE, &coeff_perez) > 0) |
574 |
< |
{ |
406 |
< |
fprintf(stderr,"lect_coeff_perez does not work\n"); |
407 |
< |
exit(2); |
408 |
< |
} |
409 |
< |
|
572 |
> |
|
573 |
> |
/* normalization factor for the relative sky luminance distribution, diffuse part*/ |
574 |
> |
|
575 |
|
if ( (lv_mod = malloc(145*sizeof(float))) == NULL) |
576 |
|
{ |
577 |
|
fprintf(stderr,"Out of memory in function main"); |
579 |
|
} |
580 |
|
|
581 |
|
/* read the angles */ |
582 |
< |
theta_o = theta_ordered("defangle.dat"); |
583 |
< |
phi_o = phi_ordered("defangle.dat"); |
582 |
> |
theta_o = defangle_theta; |
583 |
> |
phi_o = defangle_phi; |
584 |
> |
|
585 |
|
|
586 |
< |
/* parameters for the perez model */ |
586 |
> |
/* parameters for the perez model */ |
587 |
|
coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez); |
588 |
|
|
589 |
< |
/*calculation of the modelled luminance */ |
589 |
> |
|
590 |
> |
|
591 |
> |
/*calculation of the modelled luminance */ |
592 |
|
for (j=0;j<145;j++) |
593 |
|
{ |
594 |
|
theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith)); |
595 |
+ |
|
596 |
|
*(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez); |
597 |
< |
/*printf("theta, phi, lv_mod %lf\t %lf\t %lf\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j));*/ |
597 |
> |
|
598 |
> |
/* fprintf(stderr,"theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); */ |
599 |
|
} |
600 |
< |
|
600 |
> |
|
601 |
|
/* integration of luminance for the normalization factor, diffuse part of the sky*/ |
602 |
+ |
|
603 |
|
diffnormalization = integ_lv(lv_mod, theta_o); |
433 |
– |
/*printf("perez integration %lf\n", diffnormalization);*/ |
604 |
|
|
435 |
– |
|
605 |
|
|
606 |
|
|
607 |
< |
/*normalization coefficient in lumen or in watt*/ |
607 |
> |
/*normalization coefficient in lumen or in watt*/ |
608 |
|
if (output==0) |
609 |
|
{ |
610 |
< |
diffnormalization = diffusilluminance/diffnormalization/WHTEFFICACY; |
610 |
> |
diffnormalization = diffuseilluminance/diffnormalization/WHTEFFICACY; |
611 |
|
} |
612 |
|
else if (output==1) |
613 |
|
{ |
614 |
< |
diffnormalization = diffusirradiance/diffnormalization; |
614 |
> |
diffnormalization = diffuseirradiance/diffnormalization; |
615 |
|
} |
616 |
|
else if (output==2) |
617 |
|
{ |
618 |
< |
diffnormalization = diffusilluminance/diffnormalization; |
618 |
> |
diffnormalization = diffuseilluminance/diffnormalization; |
619 |
|
} |
620 |
|
|
621 |
< |
else {fprintf(stderr,"output argument : wrong number"); exit(1);} |
621 |
> |
else {fprintf(stderr,"Wrong output specification.\n"); exit(1);} |
622 |
|
|
623 |
|
|
624 |
|
|
625 |
|
|
626 |
< |
/* calculation for the solar source */ |
626 |
> |
/* calculation for the solar source */ |
627 |
|
if (output==0) |
628 |
|
solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY; |
629 |
|
|
633 |
|
else |
634 |
|
solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180))); |
635 |
|
|
467 |
– |
|
636 |
|
|
637 |
|
|
638 |
< |
/* Compute the ground radiance */ |
639 |
< |
zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez); |
640 |
< |
zenithbr*=diffnormalization; |
473 |
< |
fprintf(stderr, "gendaylit : the actual zenith radiance(W/m^2/sr) or luminance(cd/m^2) is : %.0lf\n", zenithbr); |
638 |
> |
/* Compute the ground radiance */ |
639 |
> |
zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez); |
640 |
> |
zenithbr*=diffnormalization; |
641 |
|
|
642 |
< |
if (skyclearness==1) |
642 |
> |
if (skyclearness==1) |
643 |
|
normfactor = 0.777778; |
644 |
|
|
645 |
< |
if (skyclearness>=6) |
645 |
> |
if (skyclearness>=6) |
646 |
|
{ |
647 |
|
F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]); |
648 |
|
normfactor = normsc()/F2/M_PI; |
649 |
|
} |
650 |
|
|
651 |
< |
if ( (skyclearness>1) && (skyclearness<6) ) |
651 |
> |
if ( (skyclearness>1) && (skyclearness<6) ) |
652 |
|
{ |
653 |
|
S_INTER=1; |
654 |
|
F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude)); |
655 |
|
normfactor = normsc()/F2/M_PI; |
656 |
|
} |
657 |
|
|
658 |
< |
groundbr = zenithbr*normfactor; |
492 |
< |
printf("# Ground ambient level: %.1f\n", groundbr); |
658 |
> |
groundbr = zenithbr*normfactor; |
659 |
|
|
660 |
< |
if (dosun&&(skyclearness>1)) |
661 |
< |
groundbr += 6.8e-5/M_PI*solarradiance*sundir[2]; |
660 |
> |
if (dosun&&(skyclearness>1)) |
661 |
> |
groundbr += 6.8e-5/M_PI*solarradiance*sundir[2]; |
662 |
|
|
663 |
< |
groundbr *= gprefl; |
663 |
> |
groundbr *= gprefl; |
664 |
|
|
665 |
|
|
666 |
+ |
|
667 |
+ |
if(*(c_perez+1)>0) |
668 |
+ |
{ |
669 |
+ |
if(suppress_warnings==0) |
670 |
+ |
{ fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1));} |
671 |
+ |
print_error_sky(); |
672 |
+ |
exit(0); |
673 |
+ |
} |
674 |
|
|
675 |
+ |
|
676 |
|
return; |
677 |
|
} |
678 |
|
|
680 |
|
|
681 |
|
|
682 |
|
|
683 |
+ |
double solar_sunset(int month,int day) |
684 |
+ |
{ |
685 |
+ |
float W; |
686 |
+ |
extern double s_latitude; |
687 |
+ |
W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day)))); |
688 |
+ |
return(12+(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15)); |
689 |
+ |
} |
690 |
|
|
691 |
|
|
692 |
< |
printsky() /* print out sky */ |
692 |
> |
|
693 |
> |
|
694 |
> |
double solar_sunrise(int month,int day) |
695 |
|
{ |
696 |
+ |
float W; |
697 |
+ |
extern double s_latitude; |
698 |
+ |
W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day)))); |
699 |
+ |
return(12-(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15)); |
700 |
+ |
} |
701 |
+ |
|
702 |
+ |
|
703 |
+ |
|
704 |
+ |
|
705 |
+ |
void printsky() |
706 |
+ |
{ |
707 |
+ |
|
708 |
+ |
printf("# Local solar time: %.2f\n", st); |
709 |
+ |
printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI); |
710 |
+ |
printf("# epsilon, delta, atmospheric precipitable water content : %.4f %.4f %.4f \n", skyclearness, skybrightness,atm_preci_water ); |
711 |
+ |
|
712 |
+ |
|
713 |
|
if (dosun&&(skyclearness>1)) |
714 |
< |
{ |
714 |
> |
{ |
715 |
|
printf("\nvoid light solar\n"); |
716 |
|
printf("0\n0\n"); |
717 |
|
printf("3 %.3e %.3e %.3e\n", solarradiance, solarradiance, solarradiance); |
718 |
|
printf("\nsolar source sun\n"); |
719 |
|
printf("0\n0\n"); |
720 |
|
printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle); |
721 |
< |
} |
521 |
< |
|
522 |
< |
if (dosun&&(skyclearness==1)) |
523 |
< |
{ |
721 |
> |
} else if (dosun) { |
722 |
|
printf("\nvoid light solar\n"); |
723 |
|
printf("0\n0\n"); |
724 |
|
printf("3 0.0 0.0 0.0\n"); |
725 |
|
printf("\nsolar source sun\n"); |
726 |
|
printf("0\n0\n"); |
727 |
|
printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle); |
728 |
< |
} |
728 |
> |
} |
729 |
> |
/* print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018 */ |
730 |
> |
if (color_output==1 && skyclearness < 4.5 && skyclearness >1.065 ) |
731 |
> |
{ |
732 |
> |
fprintf(stderr, " warning: sky clearness(epsilon)= %f \n",skyclearness); |
733 |
> |
fprintf(stderr, " warning: intermediate sky!! \n"); |
734 |
> |
fprintf(stderr, " warning: color model for intermediate sky pending \n"); |
735 |
> |
fprintf(stderr, " warning: no color output ! \n"); |
736 |
> |
color_output=0; |
737 |
> |
} |
738 |
> |
if (color_output==1) |
739 |
> |
{ |
740 |
> |
printf("\nvoid colorfunc skyfunc\n"); |
741 |
> |
printf("4 skybright_r skybright_g skybright_b perezlum_c.cal\n"); |
742 |
> |
printf("0\n"); |
743 |
> |
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, |
744 |
> |
*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4), |
745 |
> |
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]); |
746 |
> |
}else{ |
747 |
> |
printf("\nvoid brightfunc skyfunc\n"); |
748 |
> |
printf("2 skybright perezlum.cal\n"); |
749 |
> |
printf("0\n"); |
750 |
> |
printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr, |
751 |
> |
*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4), |
752 |
> |
sundir[0], sundir[1], sundir[2]); |
753 |
> |
} |
754 |
|
|
755 |
+ |
} |
756 |
|
|
757 |
+ |
|
758 |
+ |
|
759 |
+ |
void print_error_sky() |
760 |
+ |
{ |
761 |
+ |
|
762 |
+ |
|
763 |
+ |
sundir[0] = -sin(azimuth)*cos(altitude); |
764 |
+ |
sundir[1] = -cos(azimuth)*cos(altitude); |
765 |
+ |
sundir[2] = sin(altitude); |
766 |
+ |
|
767 |
+ |
printf("# Local solar time: %.2f\n", st); |
768 |
+ |
printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI); |
769 |
+ |
|
770 |
|
printf("\nvoid brightfunc skyfunc\n"); |
771 |
|
printf("2 skybright perezlum.cal\n"); |
772 |
|
printf("0\n"); |
773 |
< |
printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr, |
537 |
< |
*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4), |
538 |
< |
sundir[0], sundir[1], sundir[2]); |
773 |
> |
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]); |
774 |
|
} |
775 |
+ |
|
776 |
|
|
777 |
|
|
778 |
< |
printdefaults() /* print default values */ |
778 |
> |
|
779 |
> |
|
780 |
> |
void printdefaults() /* print default values */ |
781 |
|
{ |
782 |
|
printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl); |
783 |
|
if (zenithbr > 0.0) |
790 |
|
} |
791 |
|
|
792 |
|
|
793 |
< |
userror(msg) /* print usage error and quit */ |
794 |
< |
char *msg; |
793 |
> |
|
794 |
> |
|
795 |
> |
void usage_error(char* msg) /* print usage error and quit */ |
796 |
|
{ |
797 |
|
if (msg != NULL) |
798 |
< |
fprintf(stderr, "%s: Use error - %s\n", progname, msg); |
799 |
< |
fprintf(stderr, "Usage: %s month day hour [-P|-W|-L] direct_value diffus_value [options]\n", progname); |
800 |
< |
fprintf(stderr, "or : %s -ang altitude azimuth [-P|-W|-L] direct_value diffus_value [options]\n", progname); |
798 |
> |
fprintf(stderr, "%s: Use error - %s\n\n", progname, msg); |
799 |
> |
fprintf(stderr, "Usage: %s month day hour [-y year] [...]\n", progname); |
800 |
> |
fprintf(stderr, " or: %s -ang altitude azimuth [...]\n", progname); |
801 |
> |
fprintf(stderr, " followed by: -P epsilon delta [options]\n"); |
802 |
> |
fprintf(stderr, " or: [-W|-L|-G] direct_value diffuse_value [options]\n"); |
803 |
> |
fprintf(stderr, " or: -E global_irradiance [options]\n\n"); |
804 |
> |
fprintf(stderr, " Description:\n"); |
805 |
|
fprintf(stderr, " -P epsilon delta (these are the Perez parameters) \n"); |
806 |
|
fprintf(stderr, " -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n"); |
807 |
|
fprintf(stderr, " -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n"); |
808 |
|
fprintf(stderr, " -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n"); |
809 |
+ |
fprintf(stderr, " -E global-horizontal-irradiance (W/m^2)\n\n"); |
810 |
+ |
fprintf(stderr, " Output specification with option:\n"); |
811 |
|
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"); |
812 |
+ |
fprintf(stderr, " gendaylit version 2.5 (2018/04/18) \n\n"); |
813 |
|
exit(1); |
814 |
|
} |
815 |
|
|
816 |
|
|
817 |
|
|
818 |
< |
double |
819 |
< |
normsc() /* compute normalization factor (E0*F2/L0) */ |
818 |
> |
|
819 |
> |
double normsc() /* compute normalization factor (E0*F2/L0) */ |
820 |
|
{ |
821 |
|
static double nfc[2][5] = { |
822 |
|
/* clear sky approx. */ |
839 |
|
|
840 |
|
|
841 |
|
|
842 |
< |
printhead(ac, av) /* print command header */ |
843 |
< |
register int ac; |
844 |
< |
register char **av; |
842 |
> |
|
843 |
> |
|
844 |
> |
void printhead(int ac, char** av) /* print command header */ |
845 |
|
{ |
846 |
|
putchar('#'); |
847 |
|
while (ac--) { |
856 |
|
|
857 |
|
|
858 |
|
|
613 |
– |
|
614 |
– |
|
615 |
– |
|
616 |
– |
|
617 |
– |
|
618 |
– |
|
859 |
|
/* Perez models */ |
860 |
|
|
861 |
|
/* Perez global horizontal luminous efficacy model */ |
865 |
|
double value; |
866 |
|
double category_bounds[10], a[10], b[10], c[10], d[10]; |
867 |
|
int category_total_number, category_number, i; |
868 |
< |
|
869 |
< |
|
870 |
< |
if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup) |
871 |
< |
fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n"); |
872 |
< |
|
868 |
> |
|
869 |
> |
check_parametrization(); |
870 |
> |
|
871 |
> |
|
872 |
> |
/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
873 |
> |
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n"); */ |
874 |
> |
|
875 |
> |
|
876 |
|
/* initialize category bounds (clearness index bounds) */ |
877 |
|
|
878 |
|
category_total_number = 8; |
927 |
|
|
928 |
|
|
929 |
|
|
687 |
– |
|
930 |
|
for (i=1; i<=category_total_number; i++) |
931 |
|
{ |
932 |
|
if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) ) |
940 |
|
} |
941 |
|
|
942 |
|
|
943 |
+ |
|
944 |
+ |
|
945 |
|
/* global horizontal diffuse efficacy model, according to PEREZ */ |
946 |
|
double glob_h_diffuse_effi_PEREZ() |
947 |
|
{ |
949 |
|
double category_bounds[10], a[10], b[10], c[10], d[10]; |
950 |
|
int category_total_number, category_number, i; |
951 |
|
|
952 |
+ |
check_parametrization(); |
953 |
|
|
954 |
< |
if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup) |
955 |
< |
fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n"); |
956 |
< |
|
954 |
> |
|
955 |
> |
/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
956 |
> |
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */ |
957 |
> |
|
958 |
|
/* initialize category bounds (clearness index bounds) */ |
959 |
|
|
960 |
|
category_total_number = 8; |
961 |
|
|
962 |
+ |
//XXX: category_bounds > 0.1 |
963 |
|
category_bounds[1] = 1; |
964 |
|
category_bounds[2] = 1.065; |
965 |
|
category_bounds[3] = 1.230; |
1010 |
|
|
1011 |
|
|
1012 |
|
|
1013 |
< |
|
1013 |
> |
category_number = -1; |
1014 |
|
for (i=1; i<=category_total_number; i++) |
1015 |
|
{ |
1016 |
|
if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) ) |
1017 |
|
category_number = i; |
1018 |
|
} |
1019 |
|
|
1020 |
+ |
if (category_number == -1) { |
1021 |
+ |
if (suppress_warnings==0) |
1022 |
+ |
fprintf(stderr, "Warning: sky clearness (= %.3f) too high, printing error sky\n", skyclearness); |
1023 |
+ |
print_error_sky(); |
1024 |
+ |
exit(0); |
1025 |
+ |
} |
1026 |
+ |
|
1027 |
+ |
|
1028 |
|
value = a[category_number] + b[category_number]*atm_preci_water + c[category_number]*cos(sunzenith*M_PI/180) + |
1029 |
|
d[category_number]*log(skybrightness); |
1030 |
|
|
1031 |
|
return(value); |
1032 |
+ |
|
1033 |
|
} |
1034 |
|
|
1035 |
|
|
1036 |
+ |
|
1037 |
+ |
|
1038 |
+ |
|
1039 |
+ |
|
1040 |
|
/* direct normal efficacy model, according to PEREZ */ |
1041 |
|
|
1042 |
|
double direct_n_effi_PEREZ() |
1047 |
|
int category_total_number, category_number, i; |
1048 |
|
|
1049 |
|
|
1050 |
< |
if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup) |
1051 |
< |
fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n"); |
1050 |
> |
/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
1051 |
> |
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");*/ |
1052 |
|
|
1053 |
|
|
1054 |
|
/* initialize category bounds (clearness index bounds) */ |
1122 |
|
/*check the range of epsilon and delta indexes of the perez parametrization*/ |
1123 |
|
void check_parametrization() |
1124 |
|
{ |
1125 |
< |
if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup) |
1125 |
> |
|
1126 |
> |
if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) |
1127 |
|
{ |
1128 |
< |
fprintf(stderr,"sky clearness or sky brightness out of range %lf\t %lf\n", skyclearness, skybrightness); |
1129 |
< |
exit(1); |
1128 |
> |
|
1129 |
> |
/* limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */ |
1130 |
> |
|
1131 |
> |
if (skyclearness<skyclearinf){ |
1132 |
> |
/* if (suppress_warnings==0) |
1133 |
> |
fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */ |
1134 |
> |
skyclearness=skyclearinf; |
1135 |
|
} |
1136 |
+ |
if (skyclearness>skyclearsup){ |
1137 |
+ |
/* if (suppress_warnings==0) |
1138 |
+ |
fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */ |
1139 |
+ |
skyclearness=skyclearsup-0.001; |
1140 |
+ |
} |
1141 |
+ |
if (skybrightness<skybriginf){ |
1142 |
+ |
/* if (suppress_warnings==0) |
1143 |
+ |
fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */ |
1144 |
+ |
skybrightness=skybriginf; |
1145 |
+ |
} |
1146 |
+ |
if (skybrightness>skybrigsup){ |
1147 |
+ |
/* if (suppress_warnings==0) |
1148 |
+ |
fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */ |
1149 |
+ |
skybrightness=skybrigsup; |
1150 |
+ |
} |
1151 |
+ |
|
1152 |
+ |
return; } |
1153 |
|
else return; |
1154 |
|
} |
1155 |
|
|
1156 |
|
|
1157 |
< |
/* likelihood of the direct and diffuse components */ |
1157 |
> |
|
1158 |
> |
|
1159 |
> |
|
1160 |
> |
/* validity of the direct and diffuse components */ |
1161 |
|
void check_illuminances() |
1162 |
|
{ |
1163 |
< |
if (!( (directilluminance>=0) && (directilluminance<=solar_constant_l*1000) && (diffusilluminance>0) )) |
1164 |
< |
{ |
1165 |
< |
fprintf(stderr,"direct or diffuse illuminances out of range\n"); |
1166 |
< |
exit(1); |
1163 |
> |
if (directilluminance < 0) { |
1164 |
> |
if(suppress_warnings==0) |
1165 |
> |
{ fprintf(stderr,"Warning: direct illuminance < 0. Using 0.0\n"); } |
1166 |
> |
directilluminance = 0.0; |
1167 |
|
} |
1168 |
< |
return; |
1168 |
> |
if (diffuseilluminance < 0) { |
1169 |
> |
if(suppress_warnings==0) |
1170 |
> |
{ fprintf(stderr,"Warning: diffuse illuminance < 0. Using 0.0\n"); } |
1171 |
> |
diffuseilluminance = 0.0; |
1172 |
> |
} |
1173 |
> |
|
1174 |
> |
if (directilluminance+diffuseilluminance==0 && altitude > 0) { |
1175 |
> |
if(suppress_warnings==0) |
1176 |
> |
{ fprintf(stderr,"Warning: zero illuminance at sun altitude > 0, printing error sky\n"); } |
1177 |
> |
print_error_sky(); |
1178 |
> |
exit(0); |
1179 |
> |
} |
1180 |
> |
|
1181 |
> |
if (directilluminance > solar_constant_l) { |
1182 |
> |
if(suppress_warnings==0) |
1183 |
> |
{ fprintf(stderr,"Warning: direct illuminance exceeds solar constant\n"); } |
1184 |
> |
print_error_sky(); |
1185 |
> |
exit(0); |
1186 |
> |
} |
1187 |
|
} |
1188 |
|
|
1189 |
|
|
1190 |
|
void check_irradiances() |
1191 |
|
{ |
1192 |
< |
if (!( (directirradiance>=0) && (directirradiance<=solar_constant_e) && (diffusirradiance>0) )) |
1193 |
< |
{ |
1194 |
< |
fprintf(stderr,"direct or diffuse irradiances out of range\n"); |
1195 |
< |
exit(1); |
1196 |
< |
} |
1197 |
< |
return; |
1192 |
> |
if (directirradiance < 0) { |
1193 |
> |
if(suppress_warnings==0) |
1194 |
> |
{ fprintf(stderr,"Warning: direct irradiance < 0. Using 0.0\n"); } |
1195 |
> |
directirradiance = 0.0; |
1196 |
> |
} |
1197 |
> |
if (diffuseirradiance < 0) { |
1198 |
> |
if(suppress_warnings==0) |
1199 |
> |
{ fprintf(stderr,"Warning: diffuse irradiance < 0. Using 0.0\n"); } |
1200 |
> |
diffuseirradiance = 0.0; |
1201 |
> |
} |
1202 |
> |
|
1203 |
> |
if (directirradiance+diffuseirradiance==0 && altitude > 0) { |
1204 |
> |
if(suppress_warnings==0) |
1205 |
> |
{ fprintf(stderr,"Warning: zero irradiance at sun altitude > 0, printing error sky\n"); } |
1206 |
> |
print_error_sky(); |
1207 |
> |
exit(0); |
1208 |
> |
} |
1209 |
> |
|
1210 |
> |
if (directirradiance > solar_constant_e) { |
1211 |
> |
if(suppress_warnings==0) |
1212 |
> |
{ fprintf(stderr,"Warning: direct irradiance exceeds solar constant\n"); } |
1213 |
> |
print_error_sky(); |
1214 |
> |
exit(0); |
1215 |
> |
} |
1216 |
|
} |
1217 |
|
|
1218 |
|
|
1222 |
|
{ |
1223 |
|
double value; |
1224 |
|
|
1225 |
< |
value = diffusirradiance * air_mass() / ( solar_constant_e*get_eccentricity()); |
1225 |
> |
value = diffuseirradiance * air_mass() / ( solar_constant_e*get_eccentricity()); |
1226 |
|
|
1227 |
|
return(value); |
1228 |
|
} |
1231 |
|
/* Perez sky's clearness */ |
1232 |
|
double sky_clearness() |
1233 |
|
{ |
1234 |
< |
double value; |
1234 |
> |
double value; |
1235 |
|
|
1236 |
< |
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) ; |
1236 |
> |
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) ; |
1237 |
|
|
1238 |
< |
return(value); |
1238 |
> |
return(value); |
1239 |
|
} |
1240 |
|
|
1241 |
|
|
1242 |
|
|
1243 |
|
/* diffus horizontal irradiance from Perez sky's brightness */ |
1244 |
< |
double diffus_irradiance_from_sky_brightness() |
1244 |
> |
double diffuse_irradiance_from_sky_brightness() |
1245 |
|
{ |
1246 |
|
double value; |
1247 |
|
|
1256 |
|
{ |
1257 |
|
double value; |
1258 |
|
|
1259 |
< |
value = diffus_irradiance_from_sky_brightness(); |
1259 |
> |
value = diffuse_irradiance_from_sky_brightness(); |
1260 |
|
value = value * ( (skyclearness-1) * (1+1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) ); |
1261 |
|
|
1262 |
|
return(value); |
1263 |
|
} |
1264 |
|
|
1265 |
|
|
1266 |
< |
void illu_to_irra_index(void) |
1266 |
> |
|
1267 |
> |
|
1268 |
> |
void illu_to_irra_index() |
1269 |
|
{ |
1270 |
< |
double test1=0.1, test2=0.1; |
1270 |
> |
double test1=0.1, test2=0.1, d_eff; |
1271 |
|
int counter=0; |
1272 |
|
|
1273 |
< |
diffusirradiance = diffusilluminance*solar_constant_e/(solar_constant_l*1000); |
1274 |
< |
directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000); |
1273 |
> |
diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l); |
1274 |
> |
directirradiance = directilluminance*solar_constant_e/(solar_constant_l); |
1275 |
|
skyclearness = sky_clearness(); |
1276 |
|
skybrightness = sky_brightness(); |
1277 |
< |
if (skyclearness>12) skyclearness=12; |
954 |
< |
if (skybrightness<0.05) skybrightness=0.01; |
955 |
< |
|
956 |
< |
|
957 |
< |
while ( ((fabs(diffusirradiance-test1)>10) || (fabs(directirradiance-test2)>10) |
958 |
< |
|| skyclearness>skyclearinf || skyclearness<skyclearsup |
959 |
< |
|| skybrightness>skybriginf || skybrightness<skybrigsup ) |
960 |
< |
&& !(counter==5) ) |
961 |
< |
{ |
962 |
< |
/*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/ |
1277 |
> |
check_parametrization(); |
1278 |
|
|
1279 |
< |
test1=diffusirradiance; |
1279 |
> |
|
1280 |
> |
while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10) |
1281 |
> |
|| (!(skyclearness<skyclearinf || skyclearness>skyclearsup)) |
1282 |
> |
|| (!(skybrightness<skybriginf || skybrightness>skybrigsup)) ) |
1283 |
> |
&& !(counter==9) ) |
1284 |
> |
{ |
1285 |
> |
|
1286 |
> |
test1=diffuseirradiance; |
1287 |
|
test2=directirradiance; |
1288 |
|
counter++; |
1289 |
|
|
1290 |
< |
diffusirradiance = diffusilluminance/glob_h_diffuse_effi_PEREZ(); |
1291 |
< |
directirradiance = directilluminance/direct_n_effi_PEREZ(); |
970 |
< |
/*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/ |
1290 |
> |
diffuseirradiance = diffuseilluminance/glob_h_diffuse_effi_PEREZ(); |
1291 |
> |
d_eff = direct_n_effi_PEREZ(); |
1292 |
|
|
1293 |
+ |
|
1294 |
+ |
if (d_eff < 0.1) |
1295 |
+ |
directirradiance = 0; |
1296 |
+ |
else |
1297 |
+ |
directirradiance = directilluminance/d_eff; |
1298 |
+ |
|
1299 |
|
skybrightness = sky_brightness(); |
1300 |
|
skyclearness = sky_clearness(); |
1301 |
< |
if (skyclearness>12) skyclearness=12; |
1302 |
< |
if (skybrightness<0.05) skybrightness=0.01; |
976 |
< |
|
977 |
< |
/*fprintf(stderr, "%lf\t %lf\n", skybrightness, skyclearness);*/ |
978 |
< |
|
1301 |
> |
check_parametrization(); |
1302 |
> |
|
1303 |
|
} |
1304 |
|
|
1305 |
|
|
1306 |
|
return; |
1307 |
|
} |
1308 |
|
|
1309 |
< |
|
986 |
< |
int lect_coeff_perez(char *filename,float **coeff_perez) |
1309 |
> |
static int get_numlin(float epsilon) |
1310 |
|
{ |
1311 |
< |
FILE *fcoeff_perez; |
1312 |
< |
float temp; |
1313 |
< |
int i,j; |
1314 |
< |
|
1315 |
< |
if ((fcoeff_perez = frlibopen(filename)) == NULL) |
1316 |
< |
{ |
1317 |
< |
fprintf(stderr,"file %s cannot be opened\n", filename); |
1318 |
< |
return 1; /* il y a un probleme de fichier */ |
1319 |
< |
} |
1320 |
< |
else |
1321 |
< |
{ |
1322 |
< |
/*printf("file %s open\n", filename);*/ |
1323 |
< |
} |
1324 |
< |
|
1325 |
< |
rewind(fcoeff_perez); /* on se place en debut de fichier */ |
1003 |
< |
|
1004 |
< |
for (i=0;i<8;i++) |
1005 |
< |
for (j=0;j<20;j++) |
1006 |
< |
{ |
1007 |
< |
fscanf(fcoeff_perez,"%f",&temp); |
1008 |
< |
*(*coeff_perez+i*20+j) = temp; |
1009 |
< |
} |
1010 |
< |
fclose(fcoeff_perez); |
1011 |
< |
|
1012 |
< |
return 0; /* tout est OK */ |
1311 |
> |
if (epsilon < 1.065) |
1312 |
> |
return 0; |
1313 |
> |
else if (epsilon < 1.230) |
1314 |
> |
return 1; |
1315 |
> |
else if (epsilon < 1.500) |
1316 |
> |
return 2; |
1317 |
> |
else if (epsilon < 1.950) |
1318 |
> |
return 3; |
1319 |
> |
else if (epsilon < 2.800) |
1320 |
> |
return 4; |
1321 |
> |
else if (epsilon < 4.500) |
1322 |
> |
return 5; |
1323 |
> |
else if (epsilon < 6.200) |
1324 |
> |
return 6; |
1325 |
> |
return 7; |
1326 |
|
} |
1327 |
|
|
1015 |
– |
|
1016 |
– |
|
1328 |
|
/* sky luminance perez model */ |
1329 |
< |
double calc_rel_lum_perez(double dzeta,double gamma,double Z, |
1019 |
< |
double epsilon,double Delta,float *coeff_perez) |
1329 |
> |
double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]) |
1330 |
|
{ |
1331 |
+ |
|
1332 |
|
float x[5][4]; |
1333 |
|
int i,j,num_lin; |
1334 |
|
double c_perez[5]; |
1335 |
|
|
1336 |
|
if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) ) |
1337 |
|
{ |
1338 |
< |
fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n"); |
1338 |
> |
fprintf(stderr,"Error: epsilon out of range in function calc_rel_lum_perez!\n"); |
1339 |
|
exit(1); |
1340 |
|
} |
1341 |
|
|
1344 |
|
{ |
1345 |
|
if ( Delta < 0.2 ) Delta = 0.2; |
1346 |
|
} |
1347 |
< |
|
1348 |
< |
if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0; |
1349 |
< |
if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1; |
1350 |
< |
if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2; |
1040 |
< |
if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3; |
1041 |
< |
if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4; |
1042 |
< |
if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5; |
1043 |
< |
if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6; |
1044 |
< |
if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7; |
1045 |
< |
|
1347 |
> |
|
1348 |
> |
|
1349 |
> |
num_lin = get_numlin(epsilon); |
1350 |
> |
|
1351 |
|
for (i=0;i<5;i++) |
1352 |
|
for (j=0;j<4;j++) |
1353 |
|
{ |
1354 |
|
x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j); |
1355 |
< |
/* printf("x %d %d vaut %f\n",i,j,x[i][j]); */ |
1355 |
> |
/* fprintf(stderr,"x %d %d vaut %f\n",i,j,x[i][j]); */ |
1356 |
|
} |
1357 |
|
|
1358 |
|
|
1379 |
|
|
1380 |
|
|
1381 |
|
/* coefficients for the sky luminance perez model */ |
1382 |
< |
void coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez) |
1382 |
> |
void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]) |
1383 |
|
{ |
1384 |
|
float x[5][4]; |
1385 |
|
int i,j,num_lin; |
1386 |
|
|
1387 |
|
if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) ) |
1388 |
|
{ |
1389 |
< |
fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n"); |
1389 |
> |
fprintf(stderr,"Error: epsilon out of range in function coeff_lum_perez!\n"); |
1390 |
|
exit(1); |
1391 |
|
} |
1392 |
|
|
1395 |
|
{ |
1396 |
|
if ( Delta < 0.2 ) Delta = 0.2; |
1397 |
|
} |
1398 |
+ |
|
1399 |
+ |
|
1400 |
+ |
num_lin = get_numlin(epsilon); |
1401 |
|
|
1402 |
< |
if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0; |
1095 |
< |
if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1; |
1096 |
< |
if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2; |
1097 |
< |
if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3; |
1098 |
< |
if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4; |
1099 |
< |
if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5; |
1100 |
< |
if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6; |
1101 |
< |
if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7; |
1402 |
> |
/*fprintf(stderr,"numlin %d\n", num_lin);*/ |
1403 |
|
|
1404 |
|
for (i=0;i<5;i++) |
1405 |
|
for (j=0;j<4;j++) |
1431 |
|
} |
1432 |
|
|
1433 |
|
|
1434 |
+ |
|
1435 |
|
/* degrees into radians */ |
1436 |
|
double radians(double degres) |
1437 |
|
{ |
1438 |
< |
return degres*M_PI/180.0; |
1438 |
> |
return degres*(M_PI/180.); |
1439 |
|
} |
1440 |
|
|
1441 |
+ |
|
1442 |
|
/* radian into degrees */ |
1443 |
|
double degres(double radians) |
1444 |
|
{ |
1445 |
< |
return radians/M_PI*180.0; |
1445 |
> |
return radians*(180./M_PI); |
1446 |
|
} |
1447 |
|
|
1448 |
+ |
|
1449 |
|
/* calculation of the angles dzeta and gamma */ |
1450 |
|
void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z) |
1451 |
|
{ |
1455 |
|
else if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1.1 ) |
1456 |
|
{ |
1457 |
|
printf("error in calculation of gamma (angle between point and sun"); |
1458 |
< |
exit(3); |
1458 |
> |
exit(1); |
1459 |
|
} |
1460 |
|
else |
1461 |
|
*gamma = acos(cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)); |
1462 |
|
} |
1463 |
|
|
1464 |
|
|
1161 |
– |
/********************************************************************************/ |
1162 |
– |
/* Fonction: theta_ordered */ |
1163 |
– |
/* */ |
1164 |
– |
/* In: char *filename */ |
1165 |
– |
/* */ |
1166 |
– |
/* Out: float * */ |
1167 |
– |
/* */ |
1168 |
– |
/* Update: 29/08/93 */ |
1169 |
– |
/* */ |
1170 |
– |
/* Rem: theta en degres */ |
1171 |
– |
/* */ |
1172 |
– |
/* But: fournit les valeurs de theta du fichier d'entree a la memoire */ |
1173 |
– |
/* */ |
1174 |
– |
/********************************************************************************/ |
1175 |
– |
float *theta_ordered(char *filename) |
1176 |
– |
{ |
1177 |
– |
int i; |
1178 |
– |
float buffer,*ptr; |
1179 |
– |
FILE *file_in; |
1465 |
|
|
1181 |
– |
if ( (file_in = frlibopen(filename)) == NULL ) |
1182 |
– |
{ |
1183 |
– |
fprintf(stderr,"Cannot open file %s in function theta_ordered\n",filename); |
1184 |
– |
exit(1); |
1185 |
– |
} |
1186 |
– |
|
1187 |
– |
rewind(file_in); |
1188 |
– |
|
1189 |
– |
if ( (ptr = malloc(145*sizeof(float))) == NULL ) |
1190 |
– |
{ |
1191 |
– |
fprintf(stderr,"Out of memory in function theta_ordered\n"); |
1192 |
– |
exit(1); |
1193 |
– |
} |
1194 |
– |
|
1195 |
– |
for (i=0;i<145;i++) |
1196 |
– |
{ |
1197 |
– |
fscanf(file_in,"%f",&buffer); |
1198 |
– |
*(ptr+i) = buffer; |
1199 |
– |
fscanf(file_in,"%f",&buffer); |
1200 |
– |
} |
1201 |
– |
|
1202 |
– |
fclose(file_in); |
1203 |
– |
return ptr; |
1204 |
– |
} |
1205 |
– |
|
1206 |
– |
|
1207 |
– |
/********************************************************************************/ |
1208 |
– |
/* Fonction: phi_ordered */ |
1209 |
– |
/* */ |
1210 |
– |
/* In: char *filename */ |
1211 |
– |
/* */ |
1212 |
– |
/* Out: float * */ |
1213 |
– |
/* */ |
1214 |
– |
/* Update: 29/08/93 */ |
1215 |
– |
/* */ |
1216 |
– |
/* Rem: valeurs de Phi en DEGRES */ |
1217 |
– |
/* */ |
1218 |
– |
/* But: mettre les angles contenus dans le fichier d'entree dans la memoire */ |
1219 |
– |
/* */ |
1220 |
– |
/********************************************************************************/ |
1221 |
– |
float *phi_ordered(char *filename) |
1222 |
– |
{ |
1223 |
– |
int i; |
1224 |
– |
float buffer,*ptr; |
1225 |
– |
FILE *file_in; |
1226 |
– |
|
1227 |
– |
if ( (file_in = frlibopen(filename)) == NULL ) |
1228 |
– |
{ |
1229 |
– |
fprintf(stderr,"Cannot open file %s in function phi_ordered\n",filename); |
1230 |
– |
exit(1); |
1231 |
– |
} |
1232 |
– |
|
1233 |
– |
rewind(file_in); |
1234 |
– |
|
1235 |
– |
if ( (ptr = malloc(145*sizeof(float))) == NULL ) |
1236 |
– |
{ |
1237 |
– |
fprintf(stderr,"Out of memory in function phi_ordered"); |
1238 |
– |
exit(1); |
1239 |
– |
} |
1240 |
– |
|
1241 |
– |
for (i=0;i<145;i++) |
1242 |
– |
{ |
1243 |
– |
fscanf(file_in,"%f",&buffer); |
1244 |
– |
fscanf(file_in,"%f",&buffer); |
1245 |
– |
*(ptr+i) = buffer; |
1246 |
– |
} |
1247 |
– |
|
1248 |
– |
fclose(file_in); |
1249 |
– |
return ptr; |
1250 |
– |
} |
1251 |
– |
|
1252 |
– |
|
1253 |
– |
/********************************************************************************/ |
1254 |
– |
/* Fonction: integ_lv */ |
1255 |
– |
/* */ |
1256 |
– |
/* In: float *lv,*theta */ |
1257 |
– |
/* int sun_pos */ |
1258 |
– |
/* */ |
1259 |
– |
/* Out: double */ |
1260 |
– |
/* */ |
1261 |
– |
/* Update: 29/08/93 */ |
1262 |
– |
/* */ |
1263 |
– |
/* Rem: */ |
1264 |
– |
/* */ |
1265 |
– |
/* But: calcul l'integrale de luminance relative sans la dir. du soleil */ |
1266 |
– |
/* */ |
1267 |
– |
/********************************************************************************/ |
1466 |
|
double integ_lv(float *lv,float *theta) |
1467 |
|
{ |
1468 |
|
int i; |
1469 |
|
double buffer=0.0; |
1470 |
< |
|
1470 |
> |
|
1471 |
|
for (i=0;i<145;i++) |
1472 |
+ |
{ |
1473 |
|
buffer += (*(lv+i))*cos(radians(*(theta+i))); |
1474 |
< |
|
1475 |
< |
return buffer*2*M_PI/144; |
1476 |
< |
|
1474 |
> |
} |
1475 |
> |
|
1476 |
> |
return buffer*(2.*M_PI/145.); |
1477 |
|
} |
1478 |
|
|
1479 |
|
|
1480 |
|
|
1282 |
– |
|
1283 |
– |
|
1284 |
– |
|
1481 |
|
/* enter day number(double), return E0 = square(R0/R): eccentricity correction factor */ |
1482 |
|
|
1483 |
|
double get_eccentricity() |
1490 |
|
0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle); |
1491 |
|
|
1492 |
|
return (E0); |
1297 |
– |
|
1493 |
|
} |
1494 |
|
|
1495 |
|
|
1497 |
|
double air_mass() |
1498 |
|
{ |
1499 |
|
double m; |
1305 |
– |
|
1500 |
|
if (sunzenith>90) |
1501 |
|
{ |
1502 |
< |
fprintf(stderr, "solar zenith angle larger than 90� in fuction air_mass():\n the models used are not more valid\n"); |
1503 |
< |
exit(1); |
1502 |
> |
if(suppress_warnings==0) |
1503 |
> |
{ fprintf(stderr, "Warning: air mass has reached the maximal value\n"); } |
1504 |
> |
sunzenith=90; |
1505 |
|
} |
1311 |
– |
|
1506 |
|
m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) ); |
1507 |
|
return(m); |
1508 |
|
} |
1315 |
– |
|
1316 |
– |
|
1317 |
– |
double get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut) |
1318 |
– |
|
1319 |
– |
{ |
1320 |
– |
|
1321 |
– |
double angle; |
1322 |
– |
|
1323 |
– |
|
1324 |
– |
if (sun_zenith == 0) |
1325 |
– |
puts("WARNING: zenith_angle = 0 in function get_angle_sun_vert_plan"); |
1326 |
– |
|
1327 |
– |
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) ); |
1328 |
– |
angle = angle*180/M_PI; |
1329 |
– |
return(angle); |
1330 |
– |
} |
1331 |
– |
|
1332 |
– |
|
1333 |
– |
|
1334 |
– |
|
1335 |
– |
|
1509 |
|
|
1510 |
|
|
1511 |
|
|