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.5 by greg, Wed Aug 10 22:30:31 2011 UTC vs.
Revision 2.21 by greg, Thu Jan 28 19:03:15 2021 UTC

# Line 1 | Line 1
1 < #ifndef lint
2 < static const char RCSid[] = "$Id$";
3 < #endif
4 < /*        Copyright (c) 1994    *Fraunhofer Institut for Solar Energy Systems
5 < *                              Oltmannstr 5, D-79100 Freiburg, Germany
1 > /*      Copyright (c) 1994,2006 *Fraunhofer Institut for Solar Energy Systems
2 > *                              Heidenhofstr. 2, D-79110 Freiburg, Germany
3   *                              *Agence de l'Environnement et de la Maitrise de l'Energie
4   *                              Centre de Valbonne, 500 route des Lucioles, 06565 Sophia Antipolis Cedex, France
5   *                              *BOUYGUES
6   *                              1 Avenue Eugene Freyssinet, Saint-Quentin-Yvelines, France
7 + *  print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018
8 + *  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>
24 – #include  <ctype.h>
16  
26 – #include  "rtio.h"
27 – #include  "fvect.h"
17   #include  "color.h"
18 + #include  "sun.h"
19   #include  "paths.h"
20  
21 < extern int jdate(int month, int day);
32 < extern double stadj(int  jd);
33 < extern double sdec(int  jd);
34 < extern double salt(double sd, double st);
35 < extern double sazi(double sd, double st);
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();
46 – void computesky();
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  
53 – /* 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 */
65 < int     lect_coeff_perez(char *filename,float **coeff_perez);
66 < double  calc_rel_lum_perez(double dzeta,double gamma,double Z,
67 <                double epsilon,double Delta,float *coeff_perez);
68 < /* coefficients for the sky luminance perez model */
69 < void    coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez);
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);
74 – float   *theta_ordered(char *filename);
75 – float   *phi_ordered(char *filename);
76 – void    skip_comments(FILE *fp);
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();
83 – 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  
86 – /* date*/
87 – int     jdate(int month, int day);
88 –
89 –
90 –
91 –
92 –
93 – /* sun calculation constants */
94 – extern double  s_latitude;
95 – extern double  s_longitude;
96 – extern double  s_meridian;
97 –
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 */
# Line 118 | Line 130 | double  altitude, azimuth;                     /* or solar angles */
130  
131  
132   /* definition of the sky conditions through the Perez parametrization */
133 < double  skyclearness, skybrightness;
134 < double  solarradiance;  /*radiance of the sun disk and of the circumsolar area*/
135 < double  diffusilluminance, directilluminance, diffusirradiance, directirradiance;
136 < double  sunzenith, daynumber=150, atm_preci_water=2;
133 > double  skyclearness = 0;
134 > double  skybrightness = 0;
135 > double  solarradiance;
136 > double  diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance;
137 > double  sunzenith, daynumber, atm_preci_water, 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;
153 < char  *argv[];
172 >
173 > int main(int argc, char** argv)
174   {
175          int  i;
176  
177          progname = argv[0];
178          if (argc == 2 && !strcmp(argv[1], "-defaults")) {
179                  printdefaults();
180 <                exit(0);
180 >                return 0;
181          }
182          if (argc < 4)
183 <                userror("arg count");
183 >                usage_error("arg count");
184          if (!strcmp(argv[1], "-ang")) {
185                  altitude = atof(argv[2]) * (M_PI/180);
186                  azimuth = atof(argv[3]) * (M_PI/180);
# Line 168 | Line 188 | char  *argv[];
188          } else {
189                  month = atoi(argv[1]);
190                  if (month < 1 || month > 12)
191 <                        userror("bad month");
191 >                        usage_error("bad month");
192                  day = atoi(argv[2]);
193                  if (day < 1 || day > 31)
194 <                        userror("bad day");
194 >                        usage_error("bad day");
195                  hour = atof(argv[3]);
196                  if (hour < 0 || hour >= 24)
197 <                        userror("bad hour");
197 >                        usage_error("bad hour");
198                  tsolar = argv[3][0] == '+';
199          }
200          for (i = 4; i < argc; i++)
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]);
# Line 193 | Line 220 | char  *argv[];
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;
# Line 211 | Line 271 | char  *argv[];
271                          case 'm':
272                                  s_meridian = atof(argv[++i]) * (M_PI/180);
273                                  break;
214 –
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':
# Line 228 | Line 287 | char  *argv[];
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;
245 –                                
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,
256 <                    "%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 <        {
263 <                fprintf(stderr,"Out of memory error in function main !");
264 <                exit(1);
265 <        }
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);
269 –
334          computesky();
335          printsky();
336 +        return 0;
337  
273 –        exit(0);
338   }
339  
340  
341 < void
342 < 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 */
285 <        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;
287 –        double  diffusion;
352          double  normfactor;
353 +        double  erbs_s0, erbs_kt;
354  
355  
291 –
356          /* compute solar direction */
357 <
357 >                
358          if (month) {                    /* from date and time */
359 <                int  jd;
296 <                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(
315 <                    "# 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);
# Line 322 | Line 441 | computesky()                   /* compute sky parameters */
441                  
442          /* calculation for the new functions */
443          sunzenith = 90 - altitude*180/M_PI;
444 <        
326 <        
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                          }
# Line 359 | Line 478 | computesky()                   /* compute sky parameters */
478                  }
479                          
480          
481 <        else if (input==2)
481 >        else if (input==2)      /* L */
482                  {              
483                  check_illuminances();
484                  illu_to_irra_index();
# Line 367 | Line 486 | computesky()                   /* compute sky parameters */
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 <
378 <                directirradiance=directirradiance/sin(altitude);
501 >                                
502                  check_irradiances();
503                  skybrightness = sky_brightness();
504                  skyclearness =  sky_clearness();
# Line 383 | Line 506 | computesky()                   /* compute sky parameters */
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  
393 –        
394 –        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  
400 –        /* allocation dynamique de memoire pour les pointeurs */
401 –        if ( (coeff_perez = malloc(8*20*sizeof(float))) == NULL )
402 –        {
403 –                fprintf(stderr,"Out of memory error in function main !");
404 –                exit(1);
405 –        }
571  
572 < /* read the coefficients for the Perez sky luminance model */
573 <        if (lect_coeff_perez(DATFILE, &coeff_perez) > 0)
574 <        {
410 <                fprintf(stderr,"lect_coeff_perez does not work\n");
411 <                exit(2);
412 <        }
413 <
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");
# Line 418 | Line 579 | computesky()                   /* compute sky parameters */
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);
437 –        /*printf("perez integration %lf\n", diffnormalization);*/
604          
439 –        
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                  
# Line 468 | Line 633 | computesky()                   /* compute sky parameters */
633          else
634                  solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
635          
471 –                        
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;
641 < /*
642 < fprintf(stderr, "gendaylit : the actual zenith radiance(W/m^2/sr) or luminance(cd/m^2) is : %.0lf\n", zenithbr);
479 < */
480 <
481 < if (skyclearness==1)
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)
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;
498 < 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  
# Line 511 | Line 680 | return;
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 <                }
527 <                
528 <        if (dosun&&(skyclearness==1))
529 <                {              
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,
543 <            *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
544 <            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)
# Line 558 | Line 790 | printdefaults()                        /* print default values */
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. */
# Line 599 | Line 839 | normsc()                       /* compute normalization factor (E0*F2/L0)
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--) {
# Line 614 | Line 854 | register char  **av;
854  
855  
856  
617 – void
618 – skip_comments(FILE *fp)         /* skip comments in file */
619 – {
620 –        int     c;
621 –        
622 –        while ((c = getc(fp)) != EOF)
623 –                if (c == '#') {
624 –                        while ((c = getc(fp)) != EOF)
625 –                                if (c == '\n')
626 –                                        break;
627 –                } else if (!isspace(c)) {
628 –                        ungetc(c, fp);
629 –                        break;
630 –                }
631 – }
857  
858  
634 –
859   /* Perez models */
860  
861   /* Perez global horizontal luminous efficacy model */
# Line 641 | Line 865 | double glob_h_effi_PEREZ()
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;
# Line 700 | Line 927 | fprintf(stderr, "Warning : skyclearness or skybrightne
927  
928  
929  
703 –
930          for (i=1; i<=category_total_number; i++)
931          {
932                  if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
# Line 714 | Line 940 | fprintf(stderr, "Warning : skyclearness or skybrightne
940   }
941  
942  
943 +
944 +
945   /* global horizontal diffuse efficacy model, according to PEREZ */
946   double glob_h_diffuse_effi_PEREZ()
947   {
# Line 721 | Line 949 | double glob_h_diffuse_effi_PEREZ()
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;
# Line 779 | Line 1010 | fprintf(stderr, "Warning : skyclearness or skybrightne
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()
# Line 803 | Line 1047 | double         category_bounds[10], a[10], b[10], c[10], d[10
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) */
# Line 878 | Line 1122 | return(value);
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  
# Line 916 | Line 1222 | double sky_brightness()
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   }
# Line 925 | Line 1231 | return(value);
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  
# Line 950 | Line 1256 | double direct_irradiance_from_sky_clearness()
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;
970 < if (skybrightness<0.05) skybrightness=0.01;
971 <        
972 <        
973 < while ( ((fabs(diffusirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
974 <                || skyclearness>skyclearinf || skyclearness<skyclearsup
975 <                || skybrightness>skybriginf || skybrightness<skybrigsup )
976 <                 && !(counter==5) )
977 <        {
978 <                /*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
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();
986 <        /*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;
992 <
993 <                /*fprintf(stderr, "%lf\t %lf\n", skybrightness, skyclearness);*/
994 <
1301 >        check_parametrization();
1302 >                
1303          }
1304  
1305  
1306   return;
1307   }              
1308  
1309 <
1002 < 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 <        skip_comments(fcoeff_perez);
1019 <
1020 <        for (i=0;i<8;i++)
1021 <                for (j=0;j<20;j++)
1022 <                {
1023 <                        fscanf(fcoeff_perez,"%f",&temp);
1024 <                        *(*coeff_perez+i*20+j) = temp;
1025 <                }
1026 <        fclose(fcoeff_perez);
1027 <
1028 <        return 0; /* tout est OK */
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  
1031 –
1032 –
1328   /* sky luminance perez model */
1329 < double calc_rel_lum_perez(double dzeta,double gamma,double Z,
1035 < 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  
# Line 1049 | Line 1344 | double epsilon,double Delta,float *coeff_perez)
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;
1056 <        if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3;
1057 <        if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4;
1058 <        if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5;
1059 <        if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6;
1060 <        if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7;
1061 <
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  
# Line 1090 | Line 1379 | double epsilon,double Delta,float *coeff_perez)
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  
# Line 1106 | Line 1395 | void coeff_lum_perez(double Z, double epsilon, double
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;
1111 <        if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
1112 <        if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2;
1113 <        if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3;
1114 <        if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4;
1115 <        if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5;
1116 <        if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6;
1117 <        if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7;
1402 >        /*fprintf(stderr,"numlin %d\n", num_lin);*/
1403  
1404          for (i=0;i<5;i++)
1405                  for (j=0;j<4;j++)
# Line 1146 | Line 1431 | void coeff_lum_perez(double Z, double epsilon, double
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   {
# Line 1167 | Line 1455 | void theta_phi_to_dzeta_gamma(double theta,double phi,
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  
1177 – /********************************************************************************/
1178 – /*      Fonction: theta_ordered                                                 */
1179 – /*                                                                              */
1180 – /*      In: char *filename                                                      */
1181 – /*                                                                              */
1182 – /*      Out: float *                                                            */
1183 – /*                                                                              */
1184 – /*      Update: 29/08/93                                                        */
1185 – /*                                                                              */
1186 – /*      Rem: theta en degres                                                    */
1187 – /*                                                                              */
1188 – /*      But: fournit les valeurs de theta du fichier d'entree a la memoire      */
1189 – /*                                                                              */
1190 – /********************************************************************************/
1191 – float *theta_ordered(char *filename)
1192 – {
1193 –        int i;
1194 –        float buffer,*ptr;
1195 –        FILE *file_in;
1465  
1466 <        if ( (file_in = frlibopen(filename)) == NULL )
1198 <        {
1199 <                fprintf(stderr,"Cannot open file %s in function theta_ordered\n",filename);
1200 <                exit(1);
1201 <        }
1202 <        
1203 <        skip_comments(file_in);
1204 <
1205 <        if ( (ptr = malloc(145*sizeof(float))) == NULL )
1206 <        {
1207 <                fprintf(stderr,"Out of memory in function theta_ordered\n");
1208 <                exit(1);
1209 <        }
1210 <
1211 <        for (i=0;i<145;i++)
1212 <        {
1213 <                fscanf(file_in,"%f",&buffer);
1214 <                *(ptr+i) = buffer;
1215 <                fscanf(file_in,"%f",&buffer);
1216 <        }
1217 <
1218 <        fclose(file_in);
1219 <        return ptr;
1220 < }
1221 <
1222 <
1223 < /********************************************************************************/
1224 < /*      Fonction: phi_ordered                                                   */
1225 < /*                                                                              */
1226 < /*      In: char *filename                                                      */
1227 < /*                                                                              */
1228 < /*      Out: float *                                                            */
1229 < /*                                                                              */
1230 < /*      Update: 29/08/93                                                        */
1231 < /*                                                                              */
1232 < /*      Rem: valeurs de Phi en DEGRES                                           */
1233 < /*                                                                              */
1234 < /*      But: mettre les angles contenus dans le fichier d'entree dans la memoire */
1235 < /*                                                                              */
1236 < /********************************************************************************/
1237 < float *phi_ordered(char *filename)
1466 > double integ_lv(float *lv,float *theta)
1467   {
1468          int i;
1469 <        float buffer,*ptr;
1241 <        FILE *file_in;
1242 <
1243 <        if ( (file_in = frlibopen(filename)) == NULL )
1244 <        {
1245 <                fprintf(stderr,"Cannot open file %s in function phi_ordered\n",filename);
1246 <                exit(1);
1247 <        }
1469 >        double buffer=0.0;
1470          
1249 –        skip_comments(file_in);
1250 –
1251 –        if ( (ptr = malloc(145*sizeof(float))) == NULL )
1252 –        {
1253 –                fprintf(stderr,"Out of memory in function phi_ordered");
1254 –                exit(1);
1255 –        }
1256 –
1471          for (i=0;i<145;i++)
1472          {
1473 <                fscanf(file_in,"%f",&buffer);
1260 <                fscanf(file_in,"%f",&buffer);
1261 <                *(ptr+i) = buffer;
1473 >                buffer += (*(lv+i))*cos(radians(*(theta+i)));
1474          }
1475 <
1476 <        fclose(file_in);
1265 <        return ptr;
1475 >                        
1476 >        return buffer*(2.*M_PI/145.);
1477   }
1478  
1479  
1269 – /********************************************************************************/
1270 – /*      Fonction: integ_lv                                                      */
1271 – /*                                                                              */
1272 – /*      In: float *lv,*theta                                                    */
1273 – /*          int sun_pos                                                         */
1274 – /*                                                                              */
1275 – /*      Out: double                                                             */
1276 – /*                                                                              */
1277 – /*      Update: 29/08/93                                                        */
1278 – /*                                                                              */
1279 – /*      Rem:                                                                    */
1280 – /*                                                                              */
1281 – /*      But: calcul l'integrale de luminance relative sans la dir. du soleil    */
1282 – /*                                                                              */
1283 – /********************************************************************************/
1284 – double integ_lv(float *lv,float *theta)
1285 – {
1286 –        int i;
1287 –        double buffer=0.0;
1480  
1289 –        for (i=0;i<145;i++)
1290 –                buffer += (*(lv+i))*cos(radians(*(theta+i)));
1291 –
1292 –        return buffer*2*M_PI/144;
1293 –
1294 – }
1295 –
1296 –
1297 –
1298 –
1299 –
1300 –
1481   /* enter day number(double), return E0 = square(R0/R):  eccentricity correction factor  */
1482  
1483   double get_eccentricity()
# Line 1310 | Line 1490 | double get_eccentricity()
1490              0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle);
1491  
1492          return (E0);
1313 –
1493   }
1494  
1495  
# Line 1318 | Line 1497 | double get_eccentricity()
1497   double  air_mass()
1498   {
1499   double  m;
1321 –
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          }
1327 –        
1506   m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1507   return(m);
1508   }
1331 –
1332 –
1333 – double get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut)
1334 –
1335 – {
1336 –
1337 – double angle;
1338 –
1339 –
1340 – if (sun_zenith == 0)
1341 –        puts("WARNING: zenith_angle = 0 in function get_angle_sun_vert_plan");
1342 –
1343 – angle = acos( cos(sun_zenith*M_PI/180)*cos(direction_zenith*M_PI/180) + sin(sun_zenith*M_PI/180)*sin(direction_zenith*M_PI/180)*cos((sun_azimut-direction_azimut)*M_PI/180) );
1344 – angle = angle*180/M_PI;
1345 – return(angle);
1346 – }
1347 –
1348 –
1349 –
1350 –
1351 –
1509  
1510  
1511  

Diff Legend

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