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.8 by greg, Wed Oct 5 17:33:36 2011 UTC vs.
Revision 2.9 by greg, Wed Jan 30 01:02:42 2013 UTC

# Line 1 | Line 1
1   #ifndef lint
2   static const char RCSid[] = "$Id$";
3   #endif
4 < /*        Copyright (c) 1994    *Fraunhofer Institut for Solar Energy Systems
5 < *                              Oltmannstr 5, D-79100 Freiburg, Germany
4 > /*      Copyright (c) 1994,2006 *Fraunhofer Institut for Solar Energy Systems
5 > *                              Heidenhofstr. 2, D-79110 Freiburg, Germany
6   *                              *Agence de l'Environnement et de la Maitrise de l'Energie
7   *                              Centre de Valbonne, 500 route des Lucioles, 06565 Sophia Antipolis Cedex, France
8   *                              *BOUYGUES
9   *                              1 Avenue Eugene Freyssinet, Saint-Quentin-Yvelines, France
10 + *
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
16   */
17  
12 –
13 –
18   /*
19   *  gendaylit.c         program to generate the angular distribution of the daylight.
20   *                      Our zenith is along the Z-axis, the X-axis
# Line 21 | Line 25 | static const char RCSid[] = "$Id$";
25   #include  <string.h>
26   #include  <math.h>
27   #include  <stdlib.h>
24 – #include  <ctype.h>
28  
26 – #include  "rtmath.h"
27 – #include  "rtio.h"
29   #include  "color.h"
30   #include  "paths.h"
31  
32 < 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);
32 > #define  DOT(v1,v2)     (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
33  
34   double  normsc();
35  
36 < #define DATFILE         "coeff_perez.dat"
36 > /*static        char *rcsid="$Header$";*/
37  
38 + float coeff_perez[] = {
39 +        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,
40 +        -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,
41 +        -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,
42 +        -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,
43 +        -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,
44 +        -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,
45 +        -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,
46 +        -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};
47  
48  
49 + 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};
50 +
51 + 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};
52 +
53 +
54 +
55   /* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */
56   double sky_brightness();
57   double sky_clearness();
46 – void computesky();
58  
59   /* calculation of the direct and diffuse components from the Perez parametrization */
60 < double  diffus_irradiance_from_sky_brightness();
60 > double  diffuse_irradiance_from_sky_brightness();
61   double  direct_irradiance_from_sky_clearness();
62  
63  
# Line 59 | Line 70 | void   check_parametrization();
70   void    check_irradiances();
71   void    check_illuminances();
72   void    illu_to_irra_index();
73 + void    print_error_sky();
74  
75  
76   /* Perez sky luminance model */
65 – int     lect_coeff_perez(char *filename,float **coeff_perez);
77   double  calc_rel_lum_perez(double dzeta,double gamma,double Z,
78 <                double epsilon,double Delta,float *coeff_perez);
78 >                double epsilon,double Delta,float coeff_perez[]);
79   /* coefficients for the sky luminance perez model */
80 < void    coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez);
80 > void    coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]);
81   double  radians(double degres);
82   double  degres(double radians);
83   void    theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z);
84   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);
85  
86 + void printdefaults();
87 + void computesky();
88 + void printhead(int ac, char** av);
89 + void userror(char* msg);
90 + void printsky();
91  
92 + FILE * frlibopen(char* fname);
93  
94   /* astronomy and geometry*/
95   double  get_eccentricity();
96   double  air_mass();
83 – double  get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut);
97  
98 + extern int jdate(int month, int day);
99 + extern double stadj(int  jd);
100 + extern double sdec(int  jd);
101 + extern double salt(double sd, double st);
102 + extern double sazi(double sd, double st);
103  
86 – /* date*/
87 – int     jdate(int month, int day);
104  
89 –
90 –
91 –
92 –
105   /* sun calculation constants */
106   extern double  s_latitude;
107   extern double  s_longitude;
# 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;
133 > double  skyclearness = 0;
134 > double  skybrightness = 0;
135   double  solarradiance;  /*radiance of the sun disk and of the circumsolar area*/
136 < double  diffusilluminance, directilluminance, diffusirradiance, directirradiance;
136 > double  diffuseilluminance, directilluminance, diffuseirradiance, directirradiance;
137   double  sunzenith, daynumber=150, atm_preci_water=2;
138  
139 < double  diffnormalization, dirnormalization;
139 > double sunaltitude_border = 0;
140 > double  diffnormalization = 0;
141 > double dirnormalization = 0;
142   double  *c_perez;
143  
144   int     output=0;       /*define the unit of the output (sky luminance or radiance): visible watt=0, solar watt=1, lumen=2*/
145   int     input=0;        /*define the input for the calulation*/
146  
147 + int     suppress_warnings=0;
148 +
149          /* default values */
150   int  cloudy = 0;                                /* 1=standard, 2=uniform */
151   int  dosun = 1;
# Line 139 | Line 156 | int    S_INTER=0;
156  
157          /* computed values */
158   double  sundir[3];
159 < double  groundbr;
159 > double  groundbr = 0;
160   double  F2;
161   double  solarbr = 0.0;
162   int     u_solar = 0;                            /* -1=irradiance, 1=radiance */
# Line 148 | Line 165 | char  *progname;
165   char  errmsg[128];
166  
167  
168 < main(argc, argv)
152 < int  argc;
153 < char  *argv[];
168 > int main(int argc, char** argv)
169   {
170          int  i;
171  
172          progname = argv[0];
173          if (argc == 2 && !strcmp(argv[1], "-defaults")) {
174                  printdefaults();
175 <                exit(0);
175 >                return 0;
176          }
177          if (argc < 4)
178                  userror("arg count");
179          if (!strcmp(argv[1], "-ang")) {
180 <                altitude = atof(argv[2]) * (PI/180);
181 <                azimuth = atof(argv[3]) * (PI/180);
180 >                altitude = atof(argv[2]) * (M_PI/180);
181 >                azimuth = atof(argv[3]) * (M_PI/180);
182                  month = 0;
183          } else {
184                  month = atoi(argv[1]);
# Line 196 | Line 211 | char  *argv[];
211                          case 't':
212                                  betaturbidity = atof(argv[++i]);
213                                  break;
214 +                        case 'w':
215 +                                suppress_warnings = 1;
216 +                                break;                  
217                          case 'b':
218                                  zenithbr = atof(argv[++i]);
219                                  break;
# Line 203 | Line 221 | char  *argv[];
221                                  gprefl = atof(argv[++i]);
222                                  break;
223                          case 'a':
224 <                                s_latitude = atof(argv[++i]) * (PI/180);
224 >                                s_latitude = atof(argv[++i]) * (M_PI/180);
225                                  break;
226                          case 'o':
227 <                                s_longitude = atof(argv[++i]) * (PI/180);
227 >                                s_longitude = atof(argv[++i]) * (M_PI/180);
228                                  break;
229                          case 'm':
230 <                                s_meridian = atof(argv[++i]) * (PI/180);
230 >                                s_meridian = atof(argv[++i]) * (M_PI/180);
231                                  break;
214 –
232                          
233                          case 'O':
234                                  output = atof(argv[++i]);       /*define the unit of the output of the program :
235 <                                                                sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm)
219 <                                                                default is set to 0*/
235 >                                                                sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm) */
236                                  break;
237                                  
238                          case 'P':
# Line 228 | Line 244 | char  *argv[];
244                          case 'W':                                       /* direct normal Irradiance [W/m^2] */
245                                  input = 1;                              /* diffuse horizontal Irrad. [W/m^2] */
246                                  directirradiance = atof(argv[++i]);
247 <                                diffusirradiance = atof(argv[++i]);
247 >                                diffuseirradiance = atof(argv[++i]);
248                                  break;
249                                  
250                          case 'L':                                       /* direct normal Illuminance [Lux] */
251                                  input = 2;                              /* diffuse horizontal Ill. [Lux] */
252                                  directilluminance = atof(argv[++i]);
253 <                                diffusilluminance = atof(argv[++i]);
253 >                                diffuseilluminance = atof(argv[++i]);
254                                  break;
255                          
256                          case 'G':                                       /* direct horizontal Irradiance [W/m^2] */
257                                  input = 3;                              /* diffuse horizontal Irrad. [W/m^2] */
258                                  directirradiance = atof(argv[++i]);
259 <                                diffusirradiance = atof(argv[++i]);
259 >                                diffuseirradiance = atof(argv[++i]);
260                                  break;
245 –                                
261                          
262 +                        case 'l':
263 +                                sunaltitude_border = atof(argv[++i]);
264 +                                break;
265 +                        
266 +                        
267                          default:
268                                  sprintf(errmsg, "unknown option: %s", argv[i]);
269                                  userror(errmsg);
# Line 251 | Line 271 | char  *argv[];
271                  else
272                          userror("bad option");
273  
274 <        if (fabs(s_meridian-s_longitude) > 30*PI/180)
274 >        if (fabs(s_meridian-s_longitude) > 30*M_PI/180)
275                  fprintf(stderr,
276                      "%s: warning: %.1f hours btwn. standard meridian and longitude\n",
277 <                    progname, (s_longitude-s_meridian)*12/PI);
277 >                    progname, (s_longitude-s_meridian)*12/M_PI);
278  
279  
280          /* allocation dynamique de memoire pour les pointeurs */
281 <        if ( (c_perez = malloc(5*sizeof(double))) == NULL )
281 >        if ( (c_perez = calloc(5, sizeof(double))) == NULL )
282          {
283                  fprintf(stderr,"Out of memory error in function main !");
284 <                exit(1);
284 >                return 1;
285          }
286  
267 –
287          printhead(argc, argv);
288  
289          computesky();
290          printsky();
291  
292 <        exit(0);
292 >        return 0;
293   }
294  
295  
296 < void
278 < computesky()                    /* compute sky parameters */
296 > void computesky()                       /* compute sky parameters */
297   {
298  
299          /* new variables */
300          int     j;
301          float   *lv_mod;  /* 145 luminance values*/
302            /* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */
303 <        float   *theta_o, *phi_o, *coeff_perez;
303 >        float   *theta_o, *phi_o;
304          double  dzeta, gamma;
305          double  normfactor;
306  
# Line 306 | Line 324 | computesky()                   /* compute sky parameters */
324                  daynumber = (double)jdate(month, day);
325  
326          }
327 <        if (!cloudy && altitude > 87.*PI/180.) {
328 <                fprintf(stderr,
327 >        
328 >        
329 >        
330 >        
331 >        
332 >        /* if loop for the -l option. 01/2013 Sprenger  */
333 >        
334 >        if (altitude*180/M_PI < sunaltitude_border) {
335 >        
336 >        if (suppress_warnings==0)
337 >            fprintf(stderr, "Warning: sun altitude (%.3f degrees) below the border (%.3f degrees)\n",altitude*180/M_PI,sunaltitude_border);
338 >        print_error_sky();
339 >        exit(0);
340 >        }
341 >        
342 >        
343 >                        
344 >        
345 >                        
346 >        if (!cloudy && altitude > 87.*M_PI/180.) {
347 >                
348 >                if (suppress_warnings==0) {
349 >                    fprintf(stderr,
350                      "%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n",
351                      progname);
352 <                printf(
353 <                    "# warning - sun too close to zenith, reducing altitude to 87 degrees\n");
315 <                altitude = 87.*PI/180.;
352 >                }
353 >                altitude = 87.*M_PI/180.;
354          }
355 +        
356          sundir[0] = -sin(azimuth)*cos(altitude);
357          sundir[1] = -cos(azimuth)*cos(altitude);
358          sundir[2] = sin(altitude);
359          
360                  
361          /* calculation for the new functions */
362 <        sunzenith = 90 - altitude*180/PI;
362 >        sunzenith = 90 - altitude*180/M_PI;
363          
364          
365  
366 < /* compute the inputs for the calculation of the light distribution over the sky*/
366 >        /* compute the inputs for the calculation of the light distribution over the sky*/
367          if (input==0)
368                  {
369                  check_parametrization();
370 <                diffusirradiance = diffus_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
370 >                diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
371                  directirradiance = direct_irradiance_from_sky_clearness();
372                  check_irradiances();
373                  
374                  if (output==0 || output==2)
375                          {
376 <                        diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
376 >                        diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
377                          directilluminance = directirradiance*direct_n_effi_PEREZ();
378                          check_illuminances();
379                          }
# Line 350 | Line 389 | computesky()                   /* compute sky parameters */
389  
390                  if (output==0 || output==2)
391                          {
392 <                        diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
392 >                        diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
393                          directilluminance = directirradiance*direct_n_effi_PEREZ();
394                          check_illuminances();
395                          }
# Line 370 | Line 409 | computesky()                   /* compute sky parameters */
409                  {
410                          if (altitude<=0)
411                          {
412 <                        fprintf(stderr, "solar zenith angle larger than 90� \n the models used are not more valid\n");
413 <                        exit(1);
412 >                                if (suppress_warnings==0)
413 >                                     fprintf(stderr, "Warning: solar zenith angle larger than 90 degrees; using zero irradiance to proceed\n");
414 >                                directirradiance = 0;
415 >                                diffuseirradiance = 0;
416 >                        } else {
417 >                                directirradiance=directirradiance/sin(altitude);
418                          }
376 –
377 –                directirradiance=directirradiance/sin(altitude);
419                  check_irradiances();
420                  skybrightness = sky_brightness();
421                  skyclearness =  sky_clearness();
# Line 382 | Line 423 | computesky()                   /* compute sky parameters */
423  
424                  if (output==0 || output==2)
425                          {
426 <                        diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
426 >                        diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
427                          directilluminance = directirradiance*direct_n_effi_PEREZ();
428                          check_illuminances();
429                          }
# Line 394 | Line 435 | computesky()                   /* compute sky parameters */
435  
436  
437          
438 < /* normalization factor for the relative sky luminance distribution, diffuse part*/
438 >        /* normalization factor for the relative sky luminance distribution, diffuse part*/
439  
399 –        /* allocation dynamique de memoire pour les pointeurs */
400 –        if ( (coeff_perez = malloc(8*20*sizeof(float))) == NULL )
401 –        {
402 –                fprintf(stderr,"Out of memory error in function main !");
403 –                exit(1);
404 –        }
405 –
406 – /* read the coefficients for the Perez sky luminance model */
407 –        if (lect_coeff_perez(DATFILE, &coeff_perez) > 0)
408 –        {
409 –                fprintf(stderr,"lect_coeff_perez does not work\n");
410 –                exit(2);
411 –        }
412 –
440          if ( (lv_mod = malloc(145*sizeof(float))) == NULL)
441          {
442                  fprintf(stderr,"Out of memory in function main");
# Line 417 | Line 444 | computesky()                   /* compute sky parameters */
444          }
445  
446          /* read the angles */
447 <        theta_o = theta_ordered("defangle.dat");
448 <        phi_o = phi_ordered("defangle.dat");
447 >        theta_o = defangle_theta;
448 >        phi_o = defangle_phi;
449  
450 < /* parameters for the perez model */
450 >        /* parameters for the perez model */
451          coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez);
452  
453 < /*calculation of the modelled luminance */
453 >        /*calculation of the modelled luminance */
454          for (j=0;j<145;j++)
455          {
456                  theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith));
457                  *(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez);
458 <                /*printf("theta, phi, lv_mod %lf\t %lf\t %lf\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j));*/
458 >                // printf("theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j));
459          }
460  
461          /* integration of luminance for the normalization factor, diffuse part of the sky*/
# Line 438 | Line 465 | computesky()                   /* compute sky parameters */
465          
466  
467  
468 < /*normalization coefficient in lumen or in watt*/
468 >        /*normalization coefficient in lumen or in watt*/
469          if (output==0)
470                  {
471 <                diffnormalization = diffusilluminance/diffnormalization/WHTEFFICACY;
471 >                diffnormalization = diffuseilluminance/diffnormalization/WHTEFFICACY;
472                  }
473          else if (output==1)
474                  {
475 <                diffnormalization = diffusirradiance/diffnormalization;
475 >                diffnormalization = diffuseirradiance/diffnormalization;
476                  }
477          else if (output==2)
478                  {
479 <                diffnormalization = diffusilluminance/diffnormalization;
479 >                diffnormalization = diffuseilluminance/diffnormalization;
480                  }
481  
482 <        else    {fprintf(stderr,"output argument : wrong number"); exit(1);}
482 >        else    {fprintf(stderr,"Wrong output specification.\n"); exit(1);}
483  
484  
485  
486  
487 < /* calculation for the solar source */  
487 >        /* calculation for the solar source */  
488          if (output==0)
489 <                solarradiance = directilluminance/(2*PI*(1-cos(half_sun_angle*PI/180)))/WHTEFFICACY;
489 >                solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY;
490                  
491          else if (output==1)
492 <                solarradiance = directirradiance/(2*PI*(1-cos(half_sun_angle*PI/180)));
492 >                solarradiance = directirradiance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
493          
494          else
495 <                solarradiance = directilluminance/(2*PI*(1-cos(half_sun_angle*PI/180)));
495 >                solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
496          
497                          
498  
# Line 473 | Line 500 | computesky()                   /* compute sky parameters */
500   /* Compute the ground radiance */
501   zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
502   zenithbr*=diffnormalization;
503 < /*
477 < fprintf(stderr, "gendaylit : the actual zenith radiance(W/m^2/sr) or luminance(cd/m^2) is : %.0lf\n", zenithbr);
478 < */
479 <
503 >        
504   if (skyclearness==1)
505          normfactor = 0.777778;
506                  
507   if (skyclearness>=6)
508          {              
509 <        F2 = 0.274*(0.91 + 10.0*exp(-3.0*(PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]);
510 <        normfactor = normsc()/F2/PI;
509 >        F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]);
510 >        normfactor = normsc()/F2/M_PI;
511          }
512  
513   if ( (skyclearness>1) && (skyclearness<6) )
514          {
515          S_INTER=1;
516 <        F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(PI/2.0-altitude)*(.4441+1.48*altitude));
517 <        normfactor = normsc()/F2/PI;
516 >        F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude));
517 >        normfactor = normsc()/F2/M_PI;
518          }
519  
520   groundbr = zenithbr*normfactor;
497 – printf("# Ground ambient level: %.1f\n", groundbr);
521  
522   if (dosun&&(skyclearness>1))
523 < groundbr += 6.8e-5/PI*solarradiance*sundir[2];          
523 >        groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];                
524  
525   groundbr *= gprefl;
526  
# Line 509 | Line 532 | return;
532  
533  
534  
535 + void print_error_sky()
536 + {
537 +        sundir[0] = -sin(azimuth)*cos(altitude);
538 +        sundir[1] = -cos(azimuth)*cos(altitude);
539 +        sundir[2] = sin(altitude);
540 +        
541 +        printf("\nvoid brightfunc skyfunc\n");
542 +        printf("2 skybright perezlum.cal\n");
543 +        printf("0\n");
544 +        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]);
545 + }
546 +        
547  
548  
549 <
515 < printsky()                      /* print out sky */
549 > void printsky()                 /* print out sky */
550   {
551          if (dosun&&(skyclearness>1))
552 <                {              
552 >        {              
553                  printf("\nvoid light solar\n");
554                  printf("0\n0\n");
555                  printf("3 %.3e %.3e %.3e\n", solarradiance, solarradiance, solarradiance);
556                  printf("\nsolar source sun\n");
557                  printf("0\n0\n");
558                  printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
559 <                }
526 <                
527 <        if (dosun&&(skyclearness==1))
528 <                {              
559 >        } else if (dosun) {
560                  printf("\nvoid light solar\n");
561                  printf("0\n0\n");
562                  printf("3 0.0 0.0 0.0\n");
563                  printf("\nsolar source sun\n");
564                  printf("0\n0\n");
565                  printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
566 <                }
566 >        }
567          
537 –
568          printf("\nvoid brightfunc skyfunc\n");
569          printf("2 skybright perezlum.cal\n");
570          printf("0\n");
571          printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr,
572 <            *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
573 <            sundir[0], sundir[1], sundir[2]);
572 >                *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
573 >                sundir[0], sundir[1], sundir[2]);
574   }
575  
576  
577 < printdefaults()                 /* print default values */
577 > void printdefaults()                    /* print default values */
578   {
579          printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
580          if (zenithbr > 0.0)
581                  printf("-b %f\t\t\t# Zenith radiance (watts/ster/m^2\n", zenithbr);
582          else
583                  printf("-t %f\t\t\t# Atmospheric betaturbidity\n", betaturbidity);
584 <        printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/PI));
585 <        printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/PI));
586 <        printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/PI));
584 >        printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/M_PI));
585 >        printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/M_PI));
586 >        printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/M_PI));
587   }
588  
589  
590 < userror(msg)                    /* print usage error and quit */
561 < char  *msg;
590 > void userror(char* msg)                 /* print usage error and quit */
591   {
592          if (msg != NULL)
593                  fprintf(stderr, "%s: Use error - %s\n", progname, msg);
594 <        fprintf(stderr, "Usage: %s month day hour [-P|-W|-L] direct_value diffus_value [options]\n", progname);
595 <        fprintf(stderr, "or   : %s -ang altitude azimuth [-P|-W|-L] direct_value diffus_value [options]\n", progname);
594 >        fprintf(stderr, "Usage: %s month day hour [-P|-W|-L|-G] direct_value diffuse_value [options]\n", progname);
595 >        fprintf(stderr, "or:    %s -ang altitude azimuth [-P|-W|-L|-G] direct_value diffuse_value [options]\n", progname);
596          fprintf(stderr, "       -P epsilon delta  (these are the Perez parameters) \n");
597          fprintf(stderr, "       -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
598          fprintf(stderr, "       -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n");
599          fprintf(stderr, "       -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
600          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");
601 +        fprintf(stderr, "       gendaylit version 2.00 (2013/01/28)  \n");
602          exit(1);
603   }
604  
605  
606  
607 < double
578 < normsc()                        /* compute normalization factor (E0*F2/L0) */
607 > double normsc()           /* compute normalization factor (E0*F2/L0) */
608   {
609          static double  nfc[2][5] = {
610                                  /* clear sky approx. */
# Line 588 | Line 617 | normsc()                       /* compute normalization factor (E0*F2/L0)
617          register int  i;
618                                          /* polynomial approximation */
619          nf = nfc[S_INTER];
620 <        x = (altitude - PI/4.0)/(PI/4.0);
620 >        x = (altitude - M_PI/4.0)/(M_PI/4.0);
621          nsc = nf[i=4];
622          while (i--)
623                  nsc = nsc*x + nf[i];
# Line 598 | Line 627 | normsc()                       /* compute normalization factor (E0*F2/L0)
627  
628  
629  
630 < printhead(ac, av)               /* print command header */
602 < register int  ac;
603 < register char  **av;
630 > void printhead(int ac, char** av)               /* print command header */
631   {
632          putchar('#');
633          while (ac--) {
# Line 613 | Line 640 | register char  **av;
640  
641  
642  
616 – void
617 – skip_comments(FILE *fp)         /* skip comments in file */
618 – {
619 –        int     c;
620 –        
621 –        while ((c = getc(fp)) != EOF)
622 –                if (c == '#') {
623 –                        while ((c = getc(fp)) != EOF)
624 –                                if (c == '\n')
625 –                                        break;
626 –                } else if (!isspace(c)) {
627 –                        ungetc(c, fp);
628 –                        break;
629 –                }
630 – }
643  
644  
633 –
645   /* Perez models */
646  
647   /* Perez global horizontal luminous efficacy model */
# Line 642 | Line 653 | double glob_h_effi_PEREZ()
653          int     category_total_number, category_number, i;
654  
655  
656 < if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
657 < fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
656 > if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
657 >     fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n");
658  
659          /* initialize category bounds (clearness index bounds) */
660  
# Line 707 | Line 718 | fprintf(stderr, "Warning : skyclearness or skybrightne
718          }
719  
720          value = a[category_number] + b[category_number]*atm_preci_water  +
721 <            c[category_number]*cos(sunzenith*PI/180) +  d[category_number]*log(skybrightness);
721 >            c[category_number]*cos(sunzenith*M_PI/180) +  d[category_number]*log(skybrightness);
722  
723          return(value);
724   }
# Line 721 | Line 732 | double glob_h_diffuse_effi_PEREZ()
732          int     category_total_number, category_number, i;
733  
734          
735 < if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
736 < fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
735 > if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
736 >     fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_effi_PEREZ \n");
737  
738   /* initialize category bounds (clearness index bounds) */
739  
740          category_total_number = 8;
741  
742 + //XXX: category_bounds > 0.1
743          category_bounds[1] = 1;
744          category_bounds[2] = 1.065;
745          category_bounds[3] = 1.230;
# Line 779 | Line 791 | fprintf(stderr, "Warning : skyclearness or skybrightne
791  
792  
793  
794 +        category_number = -1;
795          for (i=1; i<=category_total_number; i++)
796          {
797                  if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
798                          category_number = i;
799          }
800  
801 <        value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*cos(sunzenith*PI/180) +
801 >        if (category_number == -1) {
802 >                if (suppress_warnings==0)
803 >                    fprintf(stderr, "ERROR: Model parameters out of range\n");
804 >                print_error_sky();
805 >                exit(1);
806 >        }
807 >                
808 >
809 >        value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*cos(sunzenith*M_PI/180) +
810              d[category_number]*log(skybrightness);
811  
812          return(value);
# Line 802 | Line 823 | double         category_bounds[10], a[10], b[10], c[10], d[10
823   int     category_total_number, category_number, i;
824  
825  
826 < if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
827 < fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
826 > if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
827 >   fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");
828  
829  
830   /* initialize category bounds (clearness index bounds) */
# Line 866 | Line 887 | if ( (skyclearness >= category_bounds[i]) && (skyclear
887   category_number = i;
888   }
889  
890 < value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*exp(5.73*sunzenith*PI/180 - 5) +  d[category_number]*skybrightness;
890 > value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*exp(5.73*sunzenith*M_PI/180 - 5) +  d[category_number]*skybrightness;
891  
892   if (value < 0) value = 0;
893  
# Line 877 | Line 898 | return(value);
898   /*check the range of epsilon and delta indexes of the perez parametrization*/
899   void check_parametrization()
900   {
901 < if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
901 > if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup)
902                  {
903 <                fprintf(stderr,"sky clearness or sky brightness out of range %lf\t %lf\n", skyclearness, skybrightness);
904 <                exit(1);        
903 >
904 > /*  limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */
905 >                if (skyclearness<skyclearinf){
906 >                        skyclearness=skyclearinf;
907 >                        if (suppress_warnings==0)
908 >                            fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness);
909                  }
910 +                if (skyclearness>skyclearsup){
911 +                        skyclearness=skyclearsup-0.1;
912 +                        if (suppress_warnings==0)
913 +                            fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness);
914 +                }
915 +                if (skybrightness<skybriginf){
916 +                        skybrightness=skybriginf;
917 +                        if (suppress_warnings==0)
918 +                            fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness);
919 +                }
920 +                if (skybrightness>skybrigsup){
921 +                        skybrightness=skybrigsup;
922 +                        if (suppress_warnings==0)
923 +                            fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness);
924 +                }
925 +
926 +        return; }
927          else return;
928   }
929  
930  
931 < /* likelihood of the direct and diffuse components */
931 > /* validity of the direct and diffuse components */
932   void    check_illuminances()
933   {
934 <        if (!( (directilluminance>=0) && (directilluminance<=solar_constant_l*1000) && (diffusilluminance>0) ))
935 <        {
936 <        fprintf(stderr,"direct or diffuse illuminances out of range\n");
895 <        exit(1);
934 >        if (directilluminance < 0) {
935 >                fprintf(stderr,"WARNING: direct illuminance < 0. Using 0.0\n");
936 >                directilluminance = 0.0;
937          }
938 < return;
938 >        if (diffuseilluminance < 0) {
939 >                fprintf(stderr,"WARNING: diffuse illuminance < 0. Using 0.0\n");
940 >                diffuseilluminance = 0.0;
941 >        }
942 >        if (directilluminance > solar_constant_l*1000.0) {
943 >                fprintf(stderr,"ERROR: direct illuminance exceeds solar constant\n");
944 >                exit(1);
945 >        }
946   }
947  
948  
949   void    check_irradiances()
950   {
951 <        if (!( (directirradiance>=0) && (directirradiance<=solar_constant_e) && (diffusirradiance>0) ))
952 <        {
953 <        fprintf(stderr,"direct or diffuse irradiances out of range\n");
954 <        exit(1);
955 <        }      
956 < return;
951 >        if (directirradiance < 0) {
952 >                fprintf(stderr,"WARNING: direct irradiance < 0. Using 0.0\n");
953 >                directirradiance = 0.0;
954 >        }
955 >        if (diffuseirradiance < 0) {
956 >                fprintf(stderr,"WARNING: diffuse irradiance < 0. Using 0.0\n");
957 >                diffuseirradiance = 0.0;
958 >        }
959 >        if (directirradiance > solar_constant_e) {
960 >                fprintf(stderr,"ERROR: direct irradiance exceeds solar constant\n");
961 >                exit(1);
962 >        }
963   }
964          
965  
# Line 915 | Line 969 | double sky_brightness()
969   {
970   double value;
971  
972 < value = diffusirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
972 > value = diffuseirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
973  
974   return(value);
975   }
# Line 924 | Line 978 | return(value);
978   /* Perez sky's clearness */
979   double sky_clearness()
980   {
981 < double value;
981 >        double value;
982  
983 < value = ( (diffusirradiance + directirradiance)/(diffusirradiance) + 1.041*sunzenith*PI/180*sunzenith*PI/180*sunzenith*PI/180 ) / (1 + 1.041*sunzenith*PI/180*sunzenith*PI/180*sunzenith*PI/180) ;
983 >        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) ;
984  
985 < return(value);
985 >        return(value);
986   }
987  
988  
989  
990   /* diffus horizontal irradiance from Perez sky's brightness */
991 < double diffus_irradiance_from_sky_brightness()
991 > double diffuse_irradiance_from_sky_brightness()
992   {
993          double value;
994  
# Line 949 | Line 1003 | double direct_irradiance_from_sky_clearness()
1003   {
1004          double value;
1005  
1006 <        value = diffus_irradiance_from_sky_brightness();
1007 <        value = value * ( (skyclearness-1) * (1+1.041*sunzenith*PI/180*sunzenith*PI/180*sunzenith*PI/180) );
1006 >        value = diffuse_irradiance_from_sky_brightness();
1007 >        value = value * ( (skyclearness-1) * (1+1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) );
1008  
1009          return(value);
1010   }
1011  
1012  
1013 < void illu_to_irra_index(void)
1013 > void illu_to_irra_index()
1014   {
1015 < double  test1=0.1, test2=0.1;
1015 > double  test1=0.1, test2=0.1, d_eff;
1016   int     counter=0;      
1017  
1018 < diffusirradiance = diffusilluminance*solar_constant_e/(solar_constant_l*1000);
1018 > diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l*1000);
1019   directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000);
1020   skyclearness =  sky_clearness();
1021   skybrightness = sky_brightness();
1022 < if (skyclearness>12) skyclearness=12;
969 < if (skybrightness<0.05) skybrightness=0.01;
1022 > check_parametrization();
1023          
1024 <        
972 < while ( ((fabs(diffusirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
1024 > while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
1025                  || skyclearness>skyclearinf || skyclearness<skyclearsup
1026                  || skybrightness>skybriginf || skybrightness<skybrigsup )
1027                   && !(counter==5) )
1028          {
977 –                /*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
1029  
1030 <        test1=diffusirradiance;
1030 >        test1=diffuseirradiance;
1031          test2=directirradiance;
1032          counter++;
1033          
1034 <        diffusirradiance = diffusilluminance/glob_h_diffuse_effi_PEREZ();
1035 <        directirradiance = directilluminance/direct_n_effi_PEREZ();
1036 <        /*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
1034 >        diffuseirradiance = diffuseilluminance/glob_h_diffuse_effi_PEREZ();
1035 >        d_eff = direct_n_effi_PEREZ();
1036 >        if (d_eff < 0.1)
1037 >                directirradiance = 0;
1038 >        else
1039 >                directirradiance = directilluminance/d_eff;
1040          
1041          skybrightness = sky_brightness();
1042          skyclearness =  sky_clearness();
1043 <        if (skyclearness>12) skyclearness=12;
1044 <        if (skybrightness<0.05) skybrightness=0.01;
991 <
992 <                /*fprintf(stderr, "%lf\t %lf\n", skybrightness, skyclearness);*/
993 <
1043 >        check_parametrization();
1044 >        
1045          }
1046  
1047  
1048   return;
1049   }              
1050  
1051 <
1001 < int lect_coeff_perez(char *filename,float **coeff_perez)
1051 > static int get_numlin(float epsilon)
1052   {
1053 <        FILE *fcoeff_perez;
1054 <        float temp;
1055 <        int i,j;
1056 <
1057 <        if ((fcoeff_perez = frlibopen(filename)) == NULL)
1058 <        {
1059 <                fprintf(stderr,"file %s cannot be opened\n", filename);
1060 <                return 1; /* il y a un probleme de fichier */
1061 <        }
1062 <        else
1063 <        {
1064 <                /*printf("file %s  open\n", filename);*/
1065 <        }
1066 <        
1067 <        skip_comments(fcoeff_perez);
1018 <
1019 <        for (i=0;i<8;i++)
1020 <                for (j=0;j<20;j++)
1021 <                {
1022 <                        fscanf(fcoeff_perez,"%f",&temp);
1023 <                        *(*coeff_perez+i*20+j) = temp;
1024 <                }
1025 <        fclose(fcoeff_perez);
1026 <
1027 <        return 0; /* tout est OK */
1053 >        if (epsilon < 1.065)
1054 >                return 0;
1055 >        else if (epsilon < 1.230)
1056 >                return 1;
1057 >        else if (epsilon < 1.500)
1058 >                return 2;
1059 >        else if (epsilon < 1.950)
1060 >                return 3;
1061 >        else if (epsilon < 2.800)
1062 >                return 4;
1063 >        else if (epsilon < 4.500)
1064 >                return 5;
1065 >        else if (epsilon < 6.200)
1066 >                return 6;
1067 >        return 7;
1068   }
1069  
1030 –
1031 –
1070   /* sky luminance perez model */
1071 < double calc_rel_lum_perez(double dzeta,double gamma,double Z,
1034 < double epsilon,double Delta,float *coeff_perez)
1071 > double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[])
1072   {
1073          float x[5][4];
1074          int i,j,num_lin;
# Line 1039 | Line 1076 | double epsilon,double Delta,float *coeff_perez)
1076  
1077          if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
1078          {
1079 <                fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n");
1079 >                fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez!\n");
1080                  exit(1);
1081          }
1082  
# Line 1049 | Line 1086 | double epsilon,double Delta,float *coeff_perez)
1086                  if ( Delta < 0.2 ) Delta = 0.2;
1087          }
1088  
1089 <        if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0;
1053 <        if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
1054 <        if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2;
1055 <        if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3;
1056 <        if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4;
1057 <        if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5;
1058 <        if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6;
1059 <        if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7;
1089 >        num_lin = get_numlin(epsilon);
1090  
1091          for (i=0;i<5;i++)
1092                  for (j=0;j<4;j++)
# Line 1089 | Line 1119 | double epsilon,double Delta,float *coeff_perez)
1119  
1120  
1121   /* coefficients for the sky luminance perez model */
1122 < void coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez)
1122 > void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[])
1123   {
1124          float x[5][4];
1125          int i,j,num_lin;
# Line 1106 | Line 1136 | void coeff_lum_perez(double Z, double epsilon, double
1136                  if ( Delta < 0.2 ) Delta = 0.2;
1137          }
1138  
1139 <        if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0;
1110 <        if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
1111 <        if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2;
1112 <        if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3;
1113 <        if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4;
1114 <        if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5;
1115 <        if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6;
1116 <        if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7;
1139 >        num_lin = get_numlin(epsilon);
1140  
1141 +        //fprintf(stderr,"numlin %d\n", num_lin);
1142 +
1143          for (i=0;i<5;i++)
1144                  for (j=0;j<4;j++)
1145                  {
# Line 1148 | Line 1173 | void coeff_lum_perez(double Z, double epsilon, double
1173   /* degrees into radians */
1174   double radians(double degres)
1175   {
1176 <        return degres*PI/180.0;
1176 >        return degres*M_PI/180.0;
1177   }
1178  
1179   /* radian into degrees */
1180   double degres(double radians)
1181   {
1182 <        return radians/PI*180.0;
1182 >        return radians/M_PI*180.0;
1183   }
1184  
1185   /* calculation of the angles dzeta and gamma */
# Line 1173 | Line 1198 | void theta_phi_to_dzeta_gamma(double theta,double phi,
1198   }
1199  
1200  
1176 – /********************************************************************************/
1177 – /*      Fonction: theta_ordered                                                 */
1178 – /*                                                                              */
1179 – /*      In: char *filename                                                      */
1180 – /*                                                                              */
1181 – /*      Out: float *                                                            */
1182 – /*                                                                              */
1183 – /*      Update: 29/08/93                                                        */
1184 – /*                                                                              */
1185 – /*      Rem: theta en degres                                                    */
1186 – /*                                                                              */
1187 – /*      But: fournit les valeurs de theta du fichier d'entree a la memoire      */
1188 – /*                                                                              */
1189 – /********************************************************************************/
1190 – float *theta_ordered(char *filename)
1191 – {
1192 –        int i;
1193 –        float buffer,*ptr;
1194 –        FILE *file_in;
1201  
1196 –        if ( (file_in = frlibopen(filename)) == NULL )
1197 –        {
1198 –                fprintf(stderr,"Cannot open file %s in function theta_ordered\n",filename);
1199 –                exit(1);
1200 –        }
1201 –        
1202 –        skip_comments(file_in);
1203 –
1204 –        if ( (ptr = malloc(145*sizeof(float))) == NULL )
1205 –        {
1206 –                fprintf(stderr,"Out of memory in function theta_ordered\n");
1207 –                exit(1);
1208 –        }
1209 –
1210 –        for (i=0;i<145;i++)
1211 –        {
1212 –                fscanf(file_in,"%f",&buffer);
1213 –                *(ptr+i) = buffer;
1214 –                fscanf(file_in,"%f",&buffer);
1215 –        }
1216 –
1217 –        fclose(file_in);
1218 –        return ptr;
1219 – }
1220 –
1221 –
1202   /********************************************************************************/
1223 – /*      Fonction: phi_ordered                                                   */
1224 – /*                                                                              */
1225 – /*      In: char *filename                                                      */
1226 – /*                                                                              */
1227 – /*      Out: float *                                                            */
1228 – /*                                                                              */
1229 – /*      Update: 29/08/93                                                        */
1230 – /*                                                                              */
1231 – /*      Rem: valeurs de Phi en DEGRES                                           */
1232 – /*                                                                              */
1233 – /*      But: mettre les angles contenus dans le fichier d'entree dans la memoire */
1234 – /*                                                                              */
1235 – /********************************************************************************/
1236 – float *phi_ordered(char *filename)
1237 – {
1238 –        int i;
1239 –        float buffer,*ptr;
1240 –        FILE *file_in;
1241 –
1242 –        if ( (file_in = frlibopen(filename)) == NULL )
1243 –        {
1244 –                fprintf(stderr,"Cannot open file %s in function phi_ordered\n",filename);
1245 –                exit(1);
1246 –        }
1247 –        
1248 –        skip_comments(file_in);
1249 –
1250 –        if ( (ptr = malloc(145*sizeof(float))) == NULL )
1251 –        {
1252 –                fprintf(stderr,"Out of memory in function phi_ordered");
1253 –                exit(1);
1254 –        }
1255 –
1256 –        for (i=0;i<145;i++)
1257 –        {
1258 –                fscanf(file_in,"%f",&buffer);
1259 –                fscanf(file_in,"%f",&buffer);
1260 –                *(ptr+i) = buffer;
1261 –        }
1262 –
1263 –        fclose(file_in);
1264 –        return ptr;
1265 – }
1266 –
1267 –
1268 – /********************************************************************************/
1203   /*      Fonction: integ_lv                                                      */
1204   /*                                                                              */
1205   /*      In: float *lv,*theta                                                    */
# Line 1288 | Line 1222 | double integ_lv(float *lv,float *theta)
1222          for (i=0;i<145;i++)
1223                  buffer += (*(lv+i))*cos(radians(*(theta+i)));
1224  
1225 <        return buffer*2*PI/144;
1225 >        return buffer*2*M_PI/144;
1226  
1227   }
1228  
# Line 1304 | Line 1238 | double get_eccentricity()
1238          double day_angle;
1239          double E0;
1240  
1241 <        day_angle  = 2*PI*(daynumber - 1)/365;
1241 >        day_angle  = 2*M_PI*(daynumber - 1)/365;
1242          E0         = 1.00011+0.034221*cos(day_angle)+0.00128*sin(day_angle)+
1243              0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle);
1244  
# Line 1320 | Line 1254 | double m;
1254  
1255   if (sunzenith>90)
1256          {
1257 <        fprintf(stderr, "solar zenith angle larger than 90� in fuction air_mass():\n the models used are not more valid\n");
1257 >        fprintf(stderr, "Solar zenith angle larger than 90 degrees in function air_mass()\n");
1258          exit(1);
1259          }
1260          
1261 < m = 1/( cos(sunzenith*PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1261 > m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1262   return(m);
1263   }
1330 –
1331 –
1332 – double get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut)
1333 –
1334 – {
1335 –
1336 – double angle;
1337 –
1338 –
1339 – if (sun_zenith == 0)
1340 –        puts("WARNING: zenith_angle = 0 in function get_angle_sun_vert_plan");
1341 –
1342 – angle = acos( cos(sun_zenith*PI/180)*cos(direction_zenith*PI/180) + sin(sun_zenith*PI/180)*sin(direction_zenith*PI/180)*cos((sun_azimut-direction_azimut)*PI/180) );
1343 – angle = angle*180/PI;
1344 – return(angle);
1345 – }
1346 –
1347 –
1348 –
1349 –
1350 –
1264  
1265  
1266  

Diff Legend

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