ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gendaylit.c
(Generate patch)

Comparing ray/src/gen/gendaylit.c (file contents):
Revision 2.11 by greg, Tue Apr 30 17:05:27 2013 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
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 < *
11 < *      24.1.2006                some adjustments for cygwin compilation, inclusion of RADIANCE3.7 libraries, by J. Wienold
12 < *               2011/10/08      [email protected]:
13 < *                                               - integrated coeff_perez.dat and defangles.dat
14 < *                                               - avoid some segfaults caused by out of range parameters and
15 < *                                               - numerically dangerous range checks
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 < /*
19 < *  gendaylit.c         program to generate the angular distribution of the daylight.
20 < *                      Our zenith is along the Z-axis, the X-axis
21 < *                      points east, and the Y-axis points north.
22 < */
23 <
24 < #define _USE_MATH_DEFINES
25 <
10 > #define  _USE_MATH_DEFINES
11   #include  <stdio.h>
12   #include  <string.h>
13   #include  <math.h>
14   #include  <stdlib.h>
15  
16   #include  "color.h"
17 + #include  "sun.h"
18   #include  "paths.h"
19  
20   #define  DOT(v1,v2)     (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
# Line 38 | Line 24 | double  normsc();
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,0.9738,0.2809,0.0356,-0.1246,-0.5718,0.9938,
28 <        -1.2219,-0.7730,1.4148,1.1016,-0.2054,0.0367,-3.9128,0.9156,6.9750,0.1774,6.4477,-0.1239,-1.5798,-0.5081,-1.7812,0.1080,0.2624,0.0672,-0.2190,-0.4285,
29 <        -1.1000,-0.2515,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,-2.6204,-0.0156,0.1597,0.4199,-0.5562,
30 <        -0.5484,-0.6654,-0.2672,0.7117,0.7234,-0.6219,-5.6812,2.6297,33.3389,-18.3000,-62.2500,52.0781,-3.5000,0.0016,1.1477,0.1062,0.4659,-0.3296,-0.0876,-0.0329,
31 <        -0.6000,-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,1.4062,0.3988,0.0032,0.0766,-0.0656,-0.1294,
32 <        -1.0156,-0.3670,1.0078,1.4051,0.2875,-0.5328,-3.8500,3.3750,14.0000,-0.9999,-7.1406,7.5469,-3.4000,-0.1078,-1.0750,1.5702,-0.0672,0.4016,0.3017,-0.4844,
33 <        -1.0000,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,0.2656,1.0468,-0.3788,-2.4517,1.4656,
34 <        -1.0500,0.0289,0.4260,0.3590,-0.3250,0.1156,0.7781,0.0025,31.0625,-14.5000,-46.1148,55.3750,-7.2312,0.4050,13.3500,0.6234,1.5000,-0.6426,1.8564,0.5636};
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 < float defangle_theta[] = {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, 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, 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, 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, 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, 24, 24, 24, 24, 24, 24, 24, 24, 12, 12, 12, 12, 12, 12, 0};
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[] = {0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 0, 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 0, 60, 120, 180, 240, 300, 0};
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 */
63 > /* Perez sky parametrization: epsilon and delta calculations from the direct and diffuse irradiances */
64   double sky_brightness();
65   double sky_clearness();
66  
# Line 62 | Line 68 | double sky_clearness();
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  
66 /* 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();
# Line 74 | Line 82 | void   check_illuminances();
82   void    illu_to_irra_index();
83   void    print_error_sky();
84  
85 <
78 < /* Perez sky luminance model */
79 < double  calc_rel_lum_perez(double dzeta,double gamma,double Z,
80 <                double epsilon,double Delta,float coeff_perez[]);
81 < /* coefficients for the sky luminance perez model */
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);
# Line 86 | Line 90 | void   theta_phi_to_dzeta_gamma(double theta,double phi,
90   double  integ_lv(float *lv,float *theta);
91  
92   void printdefaults();
93 + void check_sun_position();
94   void computesky();
95   void printhead(int ac, char** av);
96 < void userror(char* msg);
96 > void usage_error(char* msg);
97   void printsky();
98  
99   FILE * frlibopen(char* fname);
# Line 97 | Line 102 | FILE * frlibopen(char* fname);
102   double  get_eccentricity();
103   double  air_mass();
104  
105 < extern int jdate(int month, int day);
106 < extern double stadj(int  jd);
107 < extern double sdec(int  jd);
108 < extern double salt(double sd, double st);
104 < extern double sazi(double sd, double st);
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  
106
107 /* sun calculation constants */
108 extern double  s_latitude;
109 extern double  s_longitude;
110 extern double  s_meridian;
111
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 134 | Line 132 | double  altitude, azimuth;                     /* or solar angles */
132   /* definition of the sky conditions through the Perez parametrization */
133   double  skyclearness = 0;
134   double  skybrightness = 0;
135 < double  solarradiance;  /*radiance of the sun disk and of the circumsolar area*/
136 < double  diffuseilluminance, directilluminance, diffuseirradiance, directirradiance;
137 < double  sunzenith, daynumber=150, atm_preci_water=2;
135 > double  solarradiance;
136 > double  diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance;
137 > double  sunzenith, daynumber, atm_preci_water=2;
138  
139 < double sunaltitude_border = 0;
139 > /*double  sunaltitude_border = 0;*/
140   double  diffnormalization = 0;
141 < double dirnormalization = 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*/
146 <
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 = 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 +
173   int main(int argc, char** argv)
174   {
175          int  i;
# Line 177 | Line 180 | int main(int argc, char** argv)
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 185 | Line 188 | int main(int argc, 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 201 | Line 204 | int main(int argc, char** argv)
204                                  cloudy = 0;
205                                  dosun = argv[i][0] == '+';
206                                  break;
204                        case 'r':
207                          case 'R':
208                                  u_solar = argv[i][1] == 'R' ? -1 : 1;
209                                  solarbr = atof(argv[++i]);
# Line 210 | Line 212 | int main(int argc, 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;
# Line 233 | Line 265 | int main(int argc, char** argv)
265                                  break;
266                          
267                          case 'O':
268 <                                output = atoi(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) */
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 261 | Line 294 | int main(int argc, char** argv)
294                                  diffuseirradiance = atof(argv[++i]);
295                                  break;
296                          
297 <                        case 'l':
298 <                                sunaltitude_border = atof(argv[++i]);
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,
278 <                    "%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 */
319 >        /* dynamic memory allocation for the pointers */
320          if ( (c_perez = calloc(5, sizeof(double))) == NULL )
321 <        {
285 <                fprintf(stderr,"Out of memory error in function main !");
286 <                return 1;
287 <        }
321 >        { fprintf(stderr,"Out of memory error in function main"); return 1; }
322  
323 +        
324          printhead(argc, argv);
290
325          computesky();
326          printsky();
293
327          return 0;
328 +
329   }
330  
331  
332 < void computesky()                       /* compute sky parameters */
332 >
333 >
334 >
335 > void computesky()
336   {
337  
301        /* new variables */
338          int     j;
339 <        float   *lv_mod;  /* 145 luminance values*/
340 <          /* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */
339 >        
340 >        float   *lv_mod;  /* 145 luminance values */
341          float   *theta_o, *phi_o;
342          double  dzeta, gamma;
343          double  normfactor;
344 +        double  erbs_s0, erbs_kt;
345  
346  
310
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 320 | Line 356 | void 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 <        
333 <        
334 <        /* if loop for the -l option. 01/2013 Sprenger  */
335 <        
336 <        if (altitude*180/M_PI < sunaltitude_border) {
337 <        
338 <        if (suppress_warnings==0)
339 <            fprintf(stderr, "Warning: sun altitude (%.3f degrees) below the border (%.3f degrees)\n",altitude*180/M_PI,sunaltitude_border);
340 <        print_error_sky();
341 <        exit(0);
342 <        }
343 <        
344 <        
408 >
409 >
410                          
346        
347                        
411          if (!cloudy && altitude > 87.*M_PI/180.) {
412                  
413                  if (suppress_warnings==0) {
# Line 355 | Line 418 | void computesky()                      /* compute sky parameters */
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 362 | Line 427 | void computesky()                      /* compute sky parameters */
427                  
428          /* calculation for the new functions */
429          sunzenith = 90 - altitude*180/M_PI;
430 <        
366 <        
430 >                        
431  
432          /* compute the inputs for the calculation of the light distribution over the sky*/
433 <        if (input==0)
433 >        if (input==0)           /* P */
434                  {
435                  check_parametrization();
436                  diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
# Line 382 | Line 446 | void computesky()                      /* compute sky parameters */
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                          diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
# Line 399 | Line 464 | void 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 407 | Line 472 | void 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                                  if (suppress_warnings==0)
480 <                                     fprintf(stderr, "Warning: solar zenith angle larger than 90 degrees; using zero irradiance to proceed\n");
480 >                                     fprintf(stderr, "Warning: sun altitude < 0, proceed with irradiance values of zero\n");
481                                  directirradiance = 0;
482                                  diffuseirradiance = 0;
483                          } else {
484 <                                directirradiance=directirradiance/sin(altitude);
484 >                        
485 >                                directirradiance=directirradiance/sin(altitude);
486                          }
487 +                                
488                  check_irradiances();
489                  skybrightness = sky_brightness();
490                  skyclearness =  sky_clearness();
# Line 432 | Line 499 | void computesky()                      /* compute sky parameters */
499  
500                  }
501  
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 <        else    {fprintf(stderr,"error in giving the input arguments"); exit(1);}
555 >        else    { fprintf(stderr,"error at the input arguments"); exit(1); }
556  
557  
558          
559          /* normalization factor for the relative sky luminance distribution, diffuse part*/
560 <
560 >        
561          if ( (lv_mod = malloc(145*sizeof(float))) == NULL)
562          {
563                  fprintf(stderr,"Out of memory in function main");
# Line 448 | Line 567 | void computesky()                      /* compute sky parameters */
567          /* read the angles */
568          theta_o = defangle_theta;
569          phi_o = defangle_phi;
570 +        
571  
572          /* parameters for the perez model */
573          coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez);
574  
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 %f\t %f\t %f\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);
465        /*printf("perez integration %lf\n", diffnormalization);*/
590          
467        
591  
592  
593          /*normalization coefficient in lumen or in watt*/
# Line 496 | Line 619 | void computesky()                      /* compute sky parameters */
619          else
620                  solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
621          
499                        
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;
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)
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;
644 >        groundbr = zenithbr*normfactor;
645  
646 < if (dosun&&(skyclearness>1))
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  
665  
666  
667  
668 < void print_error_sky()
668 >
669 > double solar_sunset(int month,int day)
670   {
671 <        sundir[0] = -sin(azimuth)*cos(altitude);
672 <        sundir[1] = -cos(azimuth)*cos(altitude);
673 <        sundir[2] = sin(altitude);
674 <        
543 <        printf("\nvoid brightfunc skyfunc\n");
544 <        printf("2 skybright perezlum.cal\n");
545 <        printf("0\n");
546 <        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]);
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   }
548        
676  
677  
678 < void 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          {              
700                  printf("\nvoid light solar\n");
# Line 566 | Line 711 | void printsky()                        /* print out sky */
711                  printf("0\n0\n");
712                  printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
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,
574 <                *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
575 <                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 +
764 +
765   void printdefaults()                    /* print default values */
766   {
767          printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
# Line 589 | Line 775 | void printdefaults()                   /* print default values */
775   }
776  
777  
778 < void userror(char* msg)                 /* print usage error and quit */
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|-G] direct_value diffuse_value [options]\n", progname);
785 <        fprintf(stderr, "or:    %s -ang altitude azimuth [-P|-W|-L|-G] direct_value diffuse_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.00 (2013/01/28)  \n");
797 >        fprintf(stderr, "       gendaylit version 2.5 (2018/04/18)  \n\n");
798          exit(1);
799   }
800  
801  
802  
803 +
804   double normsc()           /* compute normalization factor (E0*F2/L0) */
805   {
806          static double  nfc[2][5] = {
# Line 629 | Line 824 | double normsc()                  /* compute normalization fac
824  
825  
826  
827 +
828 +
829   void printhead(int ac, char** av)               /* print command header */
830   {
831          putchar('#');
# Line 653 | 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) && suppress_warnings==0)
856 <     fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \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 712 | Line 912 | if ((skyclearness<skyclearinf || skyclearness>skyclear
912  
913  
914  
715
915          for (i=1; i<=category_total_number; i++)
916          {
917                  if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
# Line 726 | Line 925 | if ((skyclearness<skyclearinf || skyclearness>skyclear
925   }
926  
927  
928 +
929 +
930   /* global horizontal diffuse efficacy model, according to PEREZ */
931   double glob_h_diffuse_effi_PEREZ()
932   {
# Line 733 | 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) && suppress_warnings==0)
940 <     fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_effi_PEREZ \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
947 > //XXX:  category_bounds > 0.1
948          category_bounds[1] = 1;
949          category_bounds[2] = 1.065;
950          category_bounds[3] = 1.230;
# Line 792 | Line 995 | if ((skyclearness<skyclearinf || skyclearness>skyclear
995  
996  
997  
795
998          category_number = -1;
999          for (i=1; i<=category_total_number; i++)
1000          {
# Line 802 | Line 1004 | if ((skyclearness<skyclearinf || skyclearness>skyclear
1004  
1005          if (category_number == -1) {
1006                  if (suppress_warnings==0)
1007 <                    fprintf(stderr, "ERROR: Model parameters out of range\n");
1007 >                    fprintf(stderr, "Warning: sky clearness (= %.3f) too high, printing error sky\n", skyclearness);
1008                  print_error_sky();
1009 <                exit(1);
1009 >                exit(0);
1010          }
1011                  
1012  
# Line 812 | Line 1014 | if ((skyclearness<skyclearinf || skyclearness>skyclear
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 825 | 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) && suppress_warnings==0)
1036 <   fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \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 900 | Line 1107 | return(value);
1107   /*check the range of epsilon and delta indexes of the perez parametrization*/
1108   void check_parametrization()
1109   {
1110 +
1111   if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup)
1112                  {
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;
909                        if (suppress_warnings==0)
910                            fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness);
1120                  }
1121                  if (skyclearness>skyclearsup){
1122 <                        skyclearness=skyclearsup-0.1;
1123 <                        if (suppress_warnings==0)
1124 <                            fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness);
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;
919                        if (suppress_warnings==0)
920                            fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness);
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;
924                        if (suppress_warnings==0)
925                            fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness);
1135                  }
1136  
1137          return; }
# Line 930 | Line 1139 | if (skyclearness<skyclearinf || skyclearness>skyclears
1139   }
1140  
1141  
1142 +
1143 +
1144 +
1145   /* validity of the direct and diffuse components */
1146   void    check_illuminances()
1147   {
1148          if (directilluminance < 0) {
1149 <                fprintf(stderr,"WARNING: direct illuminance < 0. Using 0.0\n");
1149 >                if(suppress_warnings==0)
1150 >                { fprintf(stderr,"Warning: direct illuminance < 0. Using 0.0\n"); }
1151                  directilluminance = 0.0;
1152          }
1153          if (diffuseilluminance < 0) {
1154 <                fprintf(stderr,"WARNING: diffuse illuminance < 0. Using 0.0\n");
1154 >                if(suppress_warnings==0)
1155 >                { fprintf(stderr,"Warning: diffuse illuminance < 0. Using 0.0\n"); }
1156                  diffuseilluminance = 0.0;
1157          }
1158 <        if (directilluminance > solar_constant_l*1000.0) {
1159 <                fprintf(stderr,"ERROR: direct illuminance exceeds solar constant\n");
1160 <                exit(1);
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) {
1178 <                fprintf(stderr,"WARNING: direct irradiance < 0. Using 0.0\n");
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 <                fprintf(stderr,"WARNING: diffuse irradiance < 0. Using 0.0\n");
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 <                fprintf(stderr,"ERROR: direct irradiance exceeds solar constant\n");
1197 <                exit(1);
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          
# Line 1012 | Line 1248 | double direct_irradiance_from_sky_clearness()
1248   }
1249  
1250  
1251 +
1252 +
1253   void illu_to_irra_index()
1254   {
1255   double  test1=0.1, test2=0.1, d_eff;
1256   int     counter=0;      
1257  
1258 < diffuseirradiance = diffuseilluminance*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   check_parametrization();
1263 <        
1263 >
1264 >
1265   while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
1266 <                || skyclearness>skyclearinf || skyclearness<skyclearsup
1267 <                || skybrightness>skybriginf || skybrightness<skybrigsup )
1268 <                 && !(counter==5) )
1266 >                || (!(skyclearness<skyclearinf || skyclearness>skyclearsup))
1267 >                || (!(skybrightness<skybriginf || skybrightness>skybrigsup)) )
1268 >                 && !(counter==9) )
1269          {
1270 <
1270 >        
1271          test1=diffuseirradiance;
1272          test2=directirradiance;
1273          counter++;
1274          
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
1281 >        else    
1282                  directirradiance = directilluminance/d_eff;
1283          
1284          skybrightness = sky_brightness();
1285          skyclearness =  sky_clearness();
1286          check_parametrization();
1287 <        
1287 >                
1288          }
1289  
1290  
# Line 1072 | Line 1313 | static int get_numlin(float epsilon)
1313   /* sky luminance perez model */
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 1087 | Line 1329 | double calc_rel_lum_perez(double dzeta,double gamma,do
1329          {
1330                  if ( Delta < 0.2 ) Delta = 0.2;
1331          }
1332 <
1332 >        
1333 >        
1334          num_lin = get_numlin(epsilon);
1335 <
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 1128 | Line 1371 | void coeff_lum_perez(double Z, double epsilon, double
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 1137 | Line 1380 | void coeff_lum_perez(double Z, double epsilon, double
1380          {
1381                  if ( Delta < 0.2 ) Delta = 0.2;
1382          }
1383 <
1383 >        
1384 >        
1385          num_lin = get_numlin(epsilon);
1386  
1387 <        //fprintf(stderr,"numlin %d\n", num_lin);
1387 >        /*fprintf(stderr,"numlin %d\n", num_lin);*/
1388  
1389          for (i=0;i<5;i++)
1390                  for (j=0;j<4;j++)
# Line 1172 | 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 1193 | 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));
# Line 1201 | Line 1448 | void theta_phi_to_dzeta_gamma(double theta,double phi,
1448  
1449  
1450  
1204 /********************************************************************************/
1205 /*      Fonction: integ_lv                                                      */
1206 /*                                                                              */
1207 /*      In: float *lv,*theta                                                    */
1208 /*          int sun_pos                                                         */
1209 /*                                                                              */
1210 /*      Out: double                                                             */
1211 /*                                                                              */
1212 /*      Update: 29/08/93                                                        */
1213 /*                                                                              */
1214 /*      Rem:                                                                    */
1215 /*                                                                              */
1216 /*      But: calcul l'integrale de luminance relative sans la dir. du soleil    */
1217 /*                                                                              */
1218 /********************************************************************************/
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;
1228
1462   }
1463  
1464  
1465  
1233
1234
1235
1466   /* enter day number(double), return E0 = square(R0/R):  eccentricity correction factor  */
1467  
1468   double get_eccentricity()
# Line 1245 | Line 1475 | double get_eccentricity()
1475              0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle);
1476  
1477          return (E0);
1248
1478   }
1479  
1480  
# Line 1253 | Line 1482 | double get_eccentricity()
1482   double  air_mass()
1483   {
1484   double  m;
1256
1485   if (sunzenith>90)
1486          {
1487 <        fprintf(stderr, "Solar zenith angle larger than 90 degrees in function air_mass()\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          }
1262        
1491   m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1492   return(m);
1493   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines