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

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

Diff Legend

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