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

Comparing ray/src/gen/gendaylit.c (file contents):
Revision 2.1 by greg, Sat Jun 6 20:18:32 2009 UTC vs.
Revision 2.7 by greg, Wed Oct 5 17:20:55 2011 UTC

# Line 21 | Line 21 | static const char RCSid[] = "$Id$";
21   #include  <string.h>
22   #include  <math.h>
23   #include  <stdlib.h>
24 + #include  <ctype.h>
25  
26   #include  "rtio.h"
27   #include  "fvect.h"
# Line 42 | Line 43 | double  normsc();
43   /* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */
44   double sky_brightness();
45   double sky_clearness();
46 + void computesky();
47  
48   /* calculation of the direct and diffuse components from the Perez parametrization */
49   double  diffus_irradiance_from_sky_brightness();
# Line 71 | Line 73 | void   theta_phi_to_dzeta_gamma(double theta,double phi,
73   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);
77  
78  
79  
# Line 159 | Line 162 | char  *argv[];
162          if (argc < 4)
163                  userror("arg count");
164          if (!strcmp(argv[1], "-ang")) {
165 <                altitude = atof(argv[2]) * (M_PI/180);
166 <                azimuth = atof(argv[3]) * (M_PI/180);
165 >                altitude = atof(argv[2]) * (PI/180);
166 >                azimuth = atof(argv[3]) * (PI/180);
167                  month = 0;
168          } else {
169                  month = atoi(argv[1]);
# Line 200 | Line 203 | char  *argv[];
203                                  gprefl = atof(argv[++i]);
204                                  break;
205                          case 'a':
206 <                                s_latitude = atof(argv[++i]) * (M_PI/180);
206 >                                s_latitude = atof(argv[++i]) * (PI/180);
207                                  break;
208                          case 'o':
209 <                                s_longitude = atof(argv[++i]) * (M_PI/180);
209 >                                s_longitude = atof(argv[++i]) * (PI/180);
210                                  break;
211                          case 'm':
212 <                                s_meridian = atof(argv[++i]) * (M_PI/180);
212 >                                s_meridian = atof(argv[++i]) * (PI/180);
213                                  break;
214  
215                          
# Line 248 | Line 251 | char  *argv[];
251                  else
252                          userror("bad option");
253  
254 <        if (fabs(s_meridian-s_longitude) > 30*M_PI/180)
254 >        if (fabs(s_meridian-s_longitude) > 30*PI/180)
255                  fprintf(stderr,
256                      "%s: warning: %.1f hours btwn. standard meridian and longitude\n",
257 <                    progname, (s_longitude-s_meridian)*12/M_PI);
257 >                    progname, (s_longitude-s_meridian)*12/PI);
258  
259  
260          /* allocation dynamique de memoire pour les pointeurs */
# Line 271 | Line 274 | char  *argv[];
274   }
275  
276  
277 + void
278   computesky()                    /* compute sky parameters */
279   {
280  
281          /* new variables */
282 <        int     j, i;
282 >        int     j;
283          float   *lv_mod;  /* 145 luminance values*/
284            /* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */
285          float   *theta_o, *phi_o, *coeff_perez;
286          double  dzeta, gamma;
283 –        double  diffusion;
287          double  normfactor;
288  
289  
# Line 303 | Line 306 | computesky()                   /* compute sky parameters */
306                  daynumber = (double)jdate(month, day);
307  
308          }
309 <        if (!cloudy && altitude > 87.*M_PI/180.) {
309 >        if (!cloudy && altitude > 87.*PI/180.) {
310                  fprintf(stderr,
311                      "%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n",
312                      progname);
313                  printf(
314                      "# warning - sun too close to zenith, reducing altitude to 87 degrees\n");
315 <                altitude = 87.*M_PI/180.;
315 >                altitude = 87.*PI/180.;
316          }
317          sundir[0] = -sin(azimuth)*cos(altitude);
318          sundir[1] = -cos(azimuth)*cos(altitude);
# Line 317 | Line 320 | computesky()                   /* compute sky parameters */
320          
321                  
322          /* calculation for the new functions */
323 <        sunzenith = 90 - altitude*180/M_PI;
323 >        sunzenith = 90 - altitude*180/PI;
324          
325          
326  
# Line 456 | Line 459 | computesky()                   /* compute sky parameters */
459  
460   /* calculation for the solar source */  
461          if (output==0)
462 <                solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY;
462 >                solarradiance = directilluminance/(2*PI*(1-cos(half_sun_angle*PI/180)))/WHTEFFICACY;
463                  
464          else if (output==1)
465 <                solarradiance = directirradiance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
465 >                solarradiance = directirradiance/(2*PI*(1-cos(half_sun_angle*PI/180)));
466          
467          else
468 <                solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
468 >                solarradiance = directilluminance/(2*PI*(1-cos(half_sun_angle*PI/180)));
469          
470                          
471  
# Line 470 | Line 473 | computesky()                   /* compute sky parameters */
473   /* Compute the ground radiance */
474   zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
475   zenithbr*=diffnormalization;
476 + /*
477   fprintf(stderr, "gendaylit : the actual zenith radiance(W/m^2/sr) or luminance(cd/m^2) is : %.0lf\n", zenithbr);
478 <        
478 > */
479 >
480   if (skyclearness==1)
481          normfactor = 0.777778;
482                  
483   if (skyclearness>=6)
484          {              
485 <        F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]);
486 <        normfactor = normsc()/F2/M_PI;
485 >        F2 = 0.274*(0.91 + 10.0*exp(-3.0*(PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]);
486 >        normfactor = normsc()/F2/PI;
487          }
488  
489   if ( (skyclearness>1) && (skyclearness<6) )
490          {
491          S_INTER=1;
492 <        F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude));
493 <        normfactor = normsc()/F2/M_PI;
492 >        F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(PI/2.0-altitude)*(.4441+1.48*altitude));
493 >        normfactor = normsc()/F2/PI;
494          }
495  
496   groundbr = zenithbr*normfactor;
497   printf("# Ground ambient level: %.1f\n", groundbr);
498  
499   if (dosun&&(skyclearness>1))
500 < groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];                
500 > groundbr += 6.8e-5/PI*solarradiance*sundir[2];          
501  
502   groundbr *= gprefl;
503  
# Line 546 | Line 551 | printdefaults()                        /* print default values */
551                  printf("-b %f\t\t\t# Zenith radiance (watts/ster/m^2\n", zenithbr);
552          else
553                  printf("-t %f\t\t\t# Atmospheric betaturbidity\n", betaturbidity);
554 <        printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/M_PI));
555 <        printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/M_PI));
556 <        printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/M_PI));
554 >        printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/PI));
555 >        printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/PI));
556 >        printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/PI));
557   }
558  
559  
# Line 583 | Line 588 | normsc()                       /* compute normalization factor (E0*F2/L0)
588          register int  i;
589                                          /* polynomial approximation */
590          nf = nfc[S_INTER];
591 <        x = (altitude - M_PI/4.0)/(M_PI/4.0);
591 >        x = (altitude - PI/4.0)/(PI/4.0);
592          nsc = nf[i=4];
593          while (i--)
594                  nsc = nsc*x + nf[i];
# Line 608 | Line 613 | register char  **av;
613  
614  
615  
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 + }
631  
632  
633  
614 –
615 –
616 –
617 –
618 –
634   /* Perez models */
635  
636   /* Perez global horizontal luminous efficacy model */
# Line 692 | Line 707 | fprintf(stderr, "Warning : skyclearness or skybrightne
707          }
708  
709          value = a[category_number] + b[category_number]*atm_preci_water  +
710 <            c[category_number]*cos(sunzenith*M_PI/180) +  d[category_number]*log(skybrightness);
710 >            c[category_number]*cos(sunzenith*PI/180) +  d[category_number]*log(skybrightness);
711  
712          return(value);
713   }
# Line 770 | Line 785 | fprintf(stderr, "Warning : skyclearness or skybrightne
785                          category_number = i;
786          }
787  
788 <        value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*cos(sunzenith*M_PI/180) +
788 >        value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*cos(sunzenith*PI/180) +
789              d[category_number]*log(skybrightness);
790  
791          return(value);
# Line 851 | Line 866 | if ( (skyclearness >= category_bounds[i]) && (skyclear
866   category_number = i;
867   }
868  
869 < 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;
869 > value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*exp(5.73*sunzenith*PI/180 - 5) +  d[category_number]*skybrightness;
870  
871   if (value < 0) value = 0;
872  
# Line 911 | Line 926 | double sky_clearness()
926   {
927   double value;
928  
929 < 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) ;
929 > 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) ;
930  
931   return(value);
932   }
# Line 935 | Line 950 | double direct_irradiance_from_sky_clearness()
950          double value;
951  
952          value = diffus_irradiance_from_sky_brightness();
953 <        value = value * ( (skyclearness-1) * (1+1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) );
953 >        value = value * ( (skyclearness-1) * (1+1.041*sunzenith*PI/180*sunzenith*PI/180*sunzenith*PI/180) );
954  
955          return(value);
956   }
# Line 998 | Line 1013 | int lect_coeff_perez(char *filename,float **coeff_pere
1013          {
1014                  /*printf("file %s  open\n", filename);*/
1015          }
1016 +        
1017 +        skip_comments(fcoeff_perez);
1018  
1002 –        rewind(fcoeff_perez); /* on se place en debut de fichier */
1003 –
1019          for (i=0;i<8;i++)
1020                  for (j=0;j<20;j++)
1021                  {
# Line 1133 | Line 1148 | void coeff_lum_perez(double Z, double epsilon, double
1148   /* degrees into radians */
1149   double radians(double degres)
1150   {
1151 <        return degres*M_PI/180.0;
1151 >        return degres*PI/180.0;
1152   }
1153  
1154   /* radian into degrees */
1155   double degres(double radians)
1156   {
1157 <        return radians/M_PI*180.0;
1157 >        return radians/PI*180.0;
1158   }
1159  
1160   /* calculation of the angles dzeta and gamma */
# Line 1183 | Line 1198 | float *theta_ordered(char *filename)
1198                  fprintf(stderr,"Cannot open file %s in function theta_ordered\n",filename);
1199                  exit(1);
1200          }
1201 +        
1202 +        skip_comments(file_in);
1203  
1187 –        rewind(file_in);
1188 –
1204          if ( (ptr = malloc(145*sizeof(float))) == NULL )
1205          {
1206                  fprintf(stderr,"Out of memory in function theta_ordered\n");
# Line 1229 | Line 1244 | float *phi_ordered(char *filename)
1244                  fprintf(stderr,"Cannot open file %s in function phi_ordered\n",filename);
1245                  exit(1);
1246          }
1247 +        
1248 +        skip_comments(file_in);
1249  
1233 –        rewind(file_in);
1234 –
1250          if ( (ptr = malloc(145*sizeof(float))) == NULL )
1251          {
1252                  fprintf(stderr,"Out of memory in function phi_ordered");
# Line 1273 | Line 1288 | double integ_lv(float *lv,float *theta)
1288          for (i=0;i<145;i++)
1289                  buffer += (*(lv+i))*cos(radians(*(theta+i)));
1290  
1291 <        return buffer*2*M_PI/144;
1291 >        return buffer*2*PI/144;
1292  
1293   }
1294  
# Line 1289 | Line 1304 | double get_eccentricity()
1304          double day_angle;
1305          double E0;
1306  
1307 <        day_angle  = 2*M_PI*(daynumber - 1)/365;
1307 >        day_angle  = 2*PI*(daynumber - 1)/365;
1308          E0         = 1.00011+0.034221*cos(day_angle)+0.00128*sin(day_angle)+
1309              0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle);
1310  
# Line 1309 | Line 1324 | if (sunzenith>90)
1324          exit(1);
1325          }
1326          
1327 < m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1327 > m = 1/( cos(sunzenith*PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1328   return(m);
1329   }
1330  
# Line 1324 | Line 1339 | double angle;
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*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) );
1343 < angle = angle*180/M_PI;
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  

Diff Legend

– Removed lines
+ Added lines
< Changed lines (old)
> Changed lines (new)