--- ray/src/gen/gendaylit.c 2013/01/30 01:02:42 2.9 +++ ray/src/gen/gendaylit.c 2013/09/06 16:54:06 2.14 @@ -1,58 +1,62 @@ -#ifndef lint -static const char RCSid[] = "$Id: gendaylit.c,v 2.9 2013/01/30 01:02:42 greg Exp $"; -#endif /* Copyright (c) 1994,2006 *Fraunhofer Institut for Solar Energy Systems * Heidenhofstr. 2, D-79110 Freiburg, Germany * *Agence de l'Environnement et de la Maitrise de l'Energie * Centre de Valbonne, 500 route des Lucioles, 06565 Sophia Antipolis Cedex, France * *BOUYGUES * 1 Avenue Eugene Freyssinet, Saint-Quentin-Yvelines, France -* -* 24.1.2006 some adjustments for cygwin compilation, inclusion of RADIANCE3.7 libraries, by J. Wienold -* 2011/10/08 chr@riap.de: -* - integrated coeff_perez.dat and defangles.dat -* - avoid some segfaults caused by out of range parameters and -* - numerically dangerous range checks */ -/* - * gendaylit.c program to generate the angular distribution of the daylight. - * Our zenith is along the Z-axis, the X-axis - * points east, and the Y-axis points north. -*/ - #include #include #include #include #include "color.h" +#include "sun.h" #include "paths.h" #define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]) +#define _USE_MATH_DEFINES double normsc(); -/*static char *rcsid="$Header: /usr/local/cvs/radiance/ray/src/gen/gendaylit.c,v 2.9 2013/01/30 01:02:42 greg Exp $";*/ +/*static char *rcsid="$Header: /usr/local/cvs/radiance/ray/src/gen/gendaylit.c,v 2.14 2013/09/06 16:54:06 greg Exp $";*/ float coeff_perez[] = { - 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, - -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, - -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, - -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, - -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, - -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, - -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, - -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}; + 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,-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,-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,-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,-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,-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,-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,-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}; -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}; +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}; -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}; +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}; -/* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */ +/* Perez sky parametrization: epsilon and delta calculations from the direct and diffuse irradiances */ double sky_brightness(); double sky_clearness(); @@ -60,11 +64,13 @@ double sky_clearness(); double diffuse_irradiance_from_sky_brightness(); double direct_irradiance_from_sky_clearness(); +/* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : */ +/* input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */ -/* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */ double glob_h_effi_PEREZ(); double glob_h_diffuse_effi_PEREZ(); double direct_n_effi_PEREZ(); + /*likelihood check of the epsilon, delta, direct and diffuse components*/ void check_parametrization(); void check_irradiances(); @@ -72,11 +78,7 @@ void check_illuminances(); void illu_to_irra_index(); void print_error_sky(); - -/* Perez sky luminance model */ -double calc_rel_lum_perez(double dzeta,double gamma,double Z, - double epsilon,double Delta,float coeff_perez[]); -/* coefficients for the sky luminance perez model */ +double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]); void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]); double radians(double degres); double degres(double radians); @@ -84,9 +86,10 @@ void theta_phi_to_dzeta_gamma(double theta,double phi, double integ_lv(float *lv,float *theta); void printdefaults(); +void check_sun_position(); void computesky(); void printhead(int ac, char** av); -void userror(char* msg); +void usage_error(char* msg); void printsky(); FILE * frlibopen(char* fname); @@ -95,27 +98,26 @@ FILE * frlibopen(char* fname); double get_eccentricity(); double air_mass(); -extern int jdate(int month, int day); -extern double stadj(int jd); -extern double sdec(int jd); -extern double salt(double sd, double st); -extern double sazi(double sd, double st); +double solar_sunset(int month, int day); +double solar_sunrise(int month, int day); +double stadj(); +int jdate(int month, int day); /* sun calculation constants */ -extern double s_latitude; -extern double s_longitude; -extern double s_meridian; +extern double s_latitude; +extern double s_longitude; +extern double s_meridian; const double AU = 149597890E3; const double solar_constant_e = 1367; /* solar constant W/m^2 */ -const double solar_constant_l = 127.5; /* solar constant klux */ +const double solar_constant_l = 127500; /* solar constant lux */ const double half_sun_angle = 0.2665; const double half_direct_angle = 2.85; -const double skyclearinf = 1.000; /* limitations for the variation of the Perez parameters */ -const double skyclearsup = 12.1; +const double skyclearinf = 1.0; /* limitations for the variation of the Perez parameters */ +const double skyclearsup = 12.01; const double skybriginf = 0.01; const double skybrigsup = 0.6; @@ -132,23 +134,24 @@ double altitude, azimuth; /* or solar angles */ /* definition of the sky conditions through the Perez parametrization */ double skyclearness = 0; double skybrightness = 0; -double solarradiance; /*radiance of the sun disk and of the circumsolar area*/ -double diffuseilluminance, directilluminance, diffuseirradiance, directirradiance; -double sunzenith, daynumber=150, atm_preci_water=2; +double solarradiance; +double diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance; +double sunzenith, daynumber, atm_preci_water=2; -double sunaltitude_border = 0; +/*double sunaltitude_border = 0;*/ double diffnormalization = 0; -double dirnormalization = 0; +double dirnormalization = 0; double *c_perez; -int output=0; /*define the unit of the output (sky luminance or radiance): visible watt=0, solar watt=1, lumen=2*/ -int input=0; /*define the input for the calulation*/ +int output=0; /* define the unit of the output (sky luminance or radiance): */ + /* visible watt=0, solar watt=1, lumen=2 */ +int input=0; /* define the input for the calulation */ int suppress_warnings=0; /* default values */ -int cloudy = 0; /* 1=standard, 2=uniform */ -int dosun = 1; +int cloudy = 0; /* 1=standard, 2=uniform */ +int dosun = 1; double zenithbr = -1.0; double betaturbidity = 0.1; double gprefl = 0.2; @@ -160,11 +163,14 @@ double groundbr = 0; double F2; double solarbr = 0.0; int u_solar = 0; /* -1=irradiance, 1=radiance */ +float timeinterval = 0; -char *progname; -char errmsg[128]; +char *progname; +char errmsg[128]; +double st; + int main(int argc, char** argv) { int i; @@ -175,7 +181,7 @@ int main(int argc, char** argv) return 0; } if (argc < 4) - userror("arg count"); + usage_error("arg count"); if (!strcmp(argv[1], "-ang")) { altitude = atof(argv[2]) * (M_PI/180); azimuth = atof(argv[3]) * (M_PI/180); @@ -183,13 +189,13 @@ int main(int argc, char** argv) } else { month = atoi(argv[1]); if (month < 1 || month > 12) - userror("bad month"); + usage_error("bad month"); day = atoi(argv[2]); if (day < 1 || day > 31) - userror("bad day"); + usage_error("bad day"); hour = atof(argv[3]); if (hour < 0 || hour >= 24) - userror("bad hour"); + usage_error("bad hour"); tsolar = argv[3][0] == '+'; } for (i = 4; i < argc; i++) @@ -199,7 +205,6 @@ int main(int argc, char** argv) cloudy = 0; dosun = argv[i][0] == '+'; break; - case 'r': case 'R': u_solar = argv[i][1] == 'R' ? -1 : 1; solarbr = atof(argv[++i]); @@ -231,8 +236,9 @@ int main(int argc, char** argv) break; case 'O': - output = atof(argv[++i]); /*define the unit of the output of the program : - sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm) */ + output = atof(argv[++i]); /*define the unit of the output of the program: + sky and sun luminance/radiance + (0==W visible, 1==W solar radiation, 2==lm) */ break; case 'P': @@ -259,58 +265,61 @@ int main(int argc, char** argv) diffuseirradiance = atof(argv[++i]); break; - case 'l': - sunaltitude_border = atof(argv[++i]); + case 'E': /* Erbs model based on the */ + input = 4; /* global-horizontal irradiance [W/m^2] */ + globalirradiance = atof(argv[++i]); break; + case 'i': + timeinterval = atof(argv[++i]); + break; + default: sprintf(errmsg, "unknown option: %s", argv[i]); - userror(errmsg); + usage_error(errmsg); } else - userror("bad option"); + usage_error("bad option"); - if (fabs(s_meridian-s_longitude) > 30*M_PI/180) - fprintf(stderr, - "%s: warning: %.1f hours btwn. standard meridian and longitude\n", + if (month && !tsolar && fabs(s_meridian-s_longitude) > 45*M_PI/180) + fprintf(stderr,"%s: warning: %.1f hours btwn. standard meridian and longitude\n", progname, (s_longitude-s_meridian)*12/M_PI); - /* allocation dynamique de memoire pour les pointeurs */ + /* dynamic memory allocation for the pointers */ if ( (c_perez = calloc(5, sizeof(double))) == NULL ) - { - fprintf(stderr,"Out of memory error in function main !"); - return 1; - } + { fprintf(stderr,"Out of memory error in function main"); return 1; } + printhead(argc, argv); - computesky(); printsky(); - return 0; + } -void computesky() /* compute sky parameters */ + + + +void computesky() { - /* new variables */ int j; - float *lv_mod; /* 145 luminance values*/ - /* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */ + + float *lv_mod; /* 145 luminance values */ float *theta_o, *phi_o; double dzeta, gamma; double normfactor; + double erbs_s0, erbs_kt; - /* compute solar direction */ - + if (month) { /* from date and time */ int jd; - double sd, st; + double sd; jd = jdate(month, day); /* Julian date */ sd = sdec(jd); /* solar declination */ @@ -318,31 +327,58 @@ void computesky() /* compute sky parameters */ st = hour; else st = hour + stadj(jd); + + + if(timeinterval) { + + if(timeinterval<0) { + fprintf(stderr, "time interval negative\n"); + exit(1); + } + + if(fabs(solar_sunrise(month,day)-st)<=timeinterval/120) { + st= (st+timeinterval/120+solar_sunrise(month,day))/2; + if(suppress_warnings==0) + { fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); } + } + + if(fabs(solar_sunset(month,day)-st)solar_sunset(month,day)+timeinterval/120)) { + if(suppress_warnings==0) + { fprintf(stderr, "Warning: sun position too low, printing error sky at %d %d %.3f\n",month,day,hour); } + altitude = salt(sd, st); + azimuth = sazi(sd, st); + print_error_sky(); + exit(0); + } + } + else + + if(stsolar_sunset(month,day)) { + if(suppress_warnings==0) + { fprintf(stderr, "Warning: sun altitude below zero at time step %i %i %.2f, printing error sky\n",month,day,hour); } + altitude = salt(sd, st); + azimuth = sazi(sd, st); + print_error_sky(); + exit(0); + } + altitude = salt(sd, st); azimuth = sazi(sd, st); daynumber = (double)jdate(month, day); - + } - - - - /* if loop for the -l option. 01/2013 Sprenger */ - - if (altitude*180/M_PI < sunaltitude_border) { - - if (suppress_warnings==0) - fprintf(stderr, "Warning: sun altitude (%.3f degrees) below the border (%.3f degrees)\n",altitude*180/M_PI,sunaltitude_border); - print_error_sky(); - exit(0); - } - - + + - - if (!cloudy && altitude > 87.*M_PI/180.) { if (suppress_warnings==0) { @@ -353,6 +389,8 @@ void computesky() /* compute sky parameters */ altitude = 87.*M_PI/180.; } + + sundir[0] = -sin(azimuth)*cos(altitude); sundir[1] = -cos(azimuth)*cos(altitude); sundir[2] = sin(altitude); @@ -360,11 +398,10 @@ void computesky() /* compute sky parameters */ /* calculation for the new functions */ sunzenith = 90 - altitude*180/M_PI; - - + /* compute the inputs for the calculation of the light distribution over the sky*/ - if (input==0) + if (input==0) /* P */ { check_parametrization(); diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/ @@ -380,13 +417,14 @@ void computesky() /* compute sky parameters */ } - else if (input==1) + else if (input==1) /* W */ { check_irradiances(); skybrightness = sky_brightness(); skyclearness = sky_clearness(); + check_parametrization(); - + if (output==0 || output==2) { diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ @@ -397,7 +435,7 @@ void computesky() /* compute sky parameters */ } - else if (input==2) + else if (input==2) /* L */ { check_illuminances(); illu_to_irra_index(); @@ -405,17 +443,19 @@ void computesky() /* compute sky parameters */ } - else if (input==3) + else if (input==3) /* G */ { if (altitude<=0) { if (suppress_warnings==0) - fprintf(stderr, "Warning: solar zenith angle larger than 90 degrees; using zero irradiance to proceed\n"); + fprintf(stderr, "Warning: sun altitude < 0, proceed with irradiance values of zero\n"); directirradiance = 0; diffuseirradiance = 0; } else { - directirradiance=directirradiance/sin(altitude); + + directirradiance=directirradiance/sin(altitude); } + check_irradiances(); skybrightness = sky_brightness(); skyclearness = sky_clearness(); @@ -430,13 +470,65 @@ void computesky() /* compute sky parameters */ } + + else if (input==4) /* E */ /* Implementation of the Erbs model. W.Sprenger (04/13) */ + { + + if (altitude<=0) + { + if (suppress_warnings==0 && globalirradiance > 50) + fprintf(stderr, "Warning: global irradiance higher than 50 W/m^2 while the sun altitude is lower than zero\n"); + globalirradiance = 0; diffuseirradiance = 0; directirradiance = 0; + + } else { + + erbs_s0 = solar_constant_e*get_eccentricity()*sin(altitude); + + if (globalirradiance>erbs_s0) + { + if (suppress_warnings==0) + fprintf(stderr, "Warning: global irradiance is higher than the time-dependent solar constant s0\n"); + globalirradiance=erbs_s0*0.999; + } + + erbs_kt=globalirradiance/erbs_s0; + + if (erbs_kt<=0.22) diffuseirradiance=globalirradiance*(1-0.09*erbs_kt); + else if (erbs_kt<=0.8) diffuseirradiance=globalirradiance*(0.9511-0.1604*erbs_kt+4.388*pow(erbs_kt,2)-16.638*pow(erbs_kt,3)+12.336*pow(erbs_kt,4)); + else if (erbs_kt<1) diffuseirradiance=globalirradiance*(0.165); + + directirradiance=globalirradiance-diffuseirradiance; + + printf("# erbs_s0, erbs_kt, irr_dir_h, irr_diff: %.3f %.3f %.3f %.3f\n", erbs_s0, erbs_kt, directirradiance, diffuseirradiance); + printf("# WARNING: the -E option is only recommended for a rough estimation!\n"); + + directirradiance=directirradiance/sin(altitude); + + } + + check_irradiances(); + skybrightness = sky_brightness(); + skyclearness = sky_clearness(); + check_parametrization(); + + if (output==0 || output==2) + { + diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ + directilluminance = directirradiance*direct_n_effi_PEREZ(); + check_illuminances(); + } + + } + + + - else {fprintf(stderr,"error in giving the input arguments"); exit(1);} + else { fprintf(stderr,"error at the input arguments"); exit(1); } /* normalization factor for the relative sky luminance distribution, diffuse part*/ - + if ( (lv_mod = malloc(145*sizeof(float))) == NULL) { fprintf(stderr,"Out of memory in function main"); @@ -446,23 +538,27 @@ void computesky() /* compute sky parameters */ /* read the angles */ theta_o = defangle_theta; phi_o = defangle_phi; + /* parameters for the perez model */ coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez); + + /*calculation of the modelled luminance */ for (j=0;j<145;j++) { theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith)); + *(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez); - // printf("theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); + + /* fprintf(stderr,"theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); */ } - + /* integration of luminance for the normalization factor, diffuse part of the sky*/ + diffnormalization = integ_lv(lv_mod, theta_o); - /*printf("perez integration %lf\n", diffnormalization);*/ - /*normalization coefficient in lumen or in watt*/ @@ -494,60 +590,82 @@ void computesky() /* compute sky parameters */ else solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180))); - -/* Compute the ground radiance */ -zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez); -zenithbr*=diffnormalization; + /* Compute the ground radiance */ + zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez); + zenithbr*=diffnormalization; -if (skyclearness==1) + if (skyclearness==1) normfactor = 0.777778; -if (skyclearness>=6) + if (skyclearness>=6) { F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]); normfactor = normsc()/F2/M_PI; } -if ( (skyclearness>1) && (skyclearness<6) ) + if ( (skyclearness>1) && (skyclearness<6) ) { S_INTER=1; F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude)); normfactor = normsc()/F2/M_PI; } -groundbr = zenithbr*normfactor; + groundbr = zenithbr*normfactor; -if (dosun&&(skyclearness>1)) + if (dosun&&(skyclearness>1)) groundbr += 6.8e-5/M_PI*solarradiance*sundir[2]; -groundbr *= gprefl; + groundbr *= gprefl; + + if(*(c_perez+1)>0) + { + if(suppress_warnings==0) + { fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1));} + print_error_sky(); + exit(0); + } + return; } -void print_error_sky() + +double solar_sunset(int month,int day) { - sundir[0] = -sin(azimuth)*cos(altitude); - sundir[1] = -cos(azimuth)*cos(altitude); - sundir[2] = sin(altitude); - - printf("\nvoid brightfunc skyfunc\n"); - printf("2 skybright perezlum.cal\n"); - printf("0\n"); - 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]); + float W; + extern double s_latitude; + W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day)))); + return(12+(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15)); } - -void printsky() /* print out sky */ + + +double solar_sunrise(int month,int day) { + float W; + extern double s_latitude; + W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day)))); + return(12-(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15)); +} + + + + +void printsky() +{ + + printf("# Local solar time: %.2f\n", st); + printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI); + + if (dosun&&(skyclearness>1)) { printf("\nvoid light solar\n"); @@ -565,15 +683,39 @@ void printsky() /* print out sky */ printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle); } + printf("\nvoid brightfunc skyfunc\n"); printf("2 skybright perezlum.cal\n"); printf("0\n"); printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr, *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4), sundir[0], sundir[1], sundir[2]); + } + +void print_error_sky() +{ + + + sundir[0] = -sin(azimuth)*cos(altitude); + sundir[1] = -cos(azimuth)*cos(altitude); + sundir[2] = sin(altitude); + + printf("# Local solar time: %.2f\n", st); + printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI); + + printf("\nvoid brightfunc skyfunc\n"); + printf("2 skybright perezlum.cal\n"); + printf("0\n"); + 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]); +} + + + + + void printdefaults() /* print default values */ { printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl); @@ -587,23 +729,32 @@ void printdefaults() /* print default values */ } -void userror(char* msg) /* print usage error and quit */ + + +void usage_error(char* msg) /* print usage error and quit */ { if (msg != NULL) - fprintf(stderr, "%s: Use error - %s\n", progname, msg); - fprintf(stderr, "Usage: %s month day hour [-P|-W|-L|-G] direct_value diffuse_value [options]\n", progname); - fprintf(stderr, "or: %s -ang altitude azimuth [-P|-W|-L|-G] direct_value diffuse_value [options]\n", progname); + fprintf(stderr, "%s: Use error - %s\n\n", progname, msg); + fprintf(stderr, "Usage: %s month day hour [...]\n", progname); + fprintf(stderr, " or: %s -ang altitude azimuth [...]\n", progname); + fprintf(stderr, " followed by: -P epsilon delta [options]\n"); + fprintf(stderr, " or: [-W|-L|-G] direct_value diffuse_value [options]\n"); + fprintf(stderr, " or: -E global_irradiance [options]\n\n"); + fprintf(stderr, " Description:\n"); fprintf(stderr, " -P epsilon delta (these are the Perez parameters) \n"); fprintf(stderr, " -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n"); fprintf(stderr, " -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n"); fprintf(stderr, " -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n"); + fprintf(stderr, " -E global-horizontal-irradiance (W/m^2)\n\n"); + fprintf(stderr, " Output specification with option:\n"); 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"); - fprintf(stderr, " gendaylit version 2.00 (2013/01/28) \n"); + fprintf(stderr, " gendaylit version 2.4 (2013/09/04) \n\n"); exit(1); } + double normsc() /* compute normalization factor (E0*F2/L0) */ { static double nfc[2][5] = { @@ -627,6 +778,8 @@ double normsc() /* compute normalization fac + + void printhead(int ac, char** av) /* print command header */ { putchar('#'); @@ -651,11 +804,14 @@ double glob_h_effi_PEREZ() double value; double category_bounds[10], a[10], b[10], c[10], d[10]; int category_total_number, category_number, i; - - -if ((skyclearnessskyclearsup || skybrightnessskybrigsup) && suppress_warnings==0) - fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n"); - + + check_parametrization(); + + +/*if ((skyclearnessskyclearsup || skybrightnessskybrigsup) && suppress_warnings==0) + fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n"); */ + + /* initialize category bounds (clearness index bounds) */ category_total_number = 8; @@ -710,7 +866,6 @@ if ((skyclearnessskyclear - for (i=1; i<=category_total_number; i++) { if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) ) @@ -724,6 +879,8 @@ if ((skyclearnessskyclear } + + /* global horizontal diffuse efficacy model, according to PEREZ */ double glob_h_diffuse_effi_PEREZ() { @@ -731,15 +888,17 @@ double glob_h_diffuse_effi_PEREZ() double category_bounds[10], a[10], b[10], c[10], d[10]; int category_total_number, category_number, i; + check_parametrization(); -if ((skyclearnessskyclearsup || skybrightnessskybrigsup) && suppress_warnings==0) - fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_effi_PEREZ \n"); - + +/*if ((skyclearnessskyclearsup || skybrightnessskybrigsup) && suppress_warnings==0) + fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */ + /* initialize category bounds (clearness index bounds) */ category_total_number = 8; -//XXX: category_bounds > 0.1 +//XXX: category_bounds > 0.1 category_bounds[1] = 1; category_bounds[2] = 1.065; category_bounds[3] = 1.230; @@ -790,7 +949,6 @@ if ((skyclearnessskyclear - category_number = -1; for (i=1; i<=category_total_number; i++) { @@ -800,9 +958,9 @@ if ((skyclearnessskyclear if (category_number == -1) { if (suppress_warnings==0) - fprintf(stderr, "ERROR: Model parameters out of range\n"); + fprintf(stderr, "Warning: sky clearness (= %.3f) too high, printing error sky\n", skyclearness); print_error_sky(); - exit(1); + exit(0); } @@ -810,9 +968,14 @@ if ((skyclearnessskyclear d[category_number]*log(skybrightness); return(value); + } + + + + /* direct normal efficacy model, according to PEREZ */ double direct_n_effi_PEREZ() @@ -823,8 +986,8 @@ double category_bounds[10], a[10], b[10], c[10], d[10 int category_total_number, category_number, i; -if ((skyclearnessskyclearsup || skybrightnessskybrigsup) && suppress_warnings==0) - fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n"); +/*if ((skyclearnessskyclearsup || skybrightnessskybrigsup) && suppress_warnings==0) + fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");*/ /* initialize category bounds (clearness index bounds) */ @@ -898,29 +1061,31 @@ return(value); /*check the range of epsilon and delta indexes of the perez parametrization*/ void check_parametrization() { + if (skyclearnessskyclearsup || skybrightnessskybrigsup) { /* limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */ + if (skyclearnessskyclearsup){ - skyclearness=skyclearsup-0.1; - if (suppress_warnings==0) - fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); + /* if (suppress_warnings==0) + fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */ + skyclearness=skyclearsup-0.001; } if (skybrightnessskybrigsup){ + /* if (suppress_warnings==0) + fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */ skybrightness=skybrigsup; - if (suppress_warnings==0) - fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); } return; } @@ -928,37 +1093,64 @@ if (skyclearnessskyclears } + + + /* validity of the direct and diffuse components */ void check_illuminances() { if (directilluminance < 0) { - fprintf(stderr,"WARNING: direct illuminance < 0. Using 0.0\n"); + if(suppress_warnings==0) + { fprintf(stderr,"Warning: direct illuminance < 0. Using 0.0\n"); } directilluminance = 0.0; } if (diffuseilluminance < 0) { - fprintf(stderr,"WARNING: diffuse illuminance < 0. Using 0.0\n"); + if(suppress_warnings==0) + { fprintf(stderr,"Warning: diffuse illuminance < 0. Using 0.0\n"); } diffuseilluminance = 0.0; } - if (directilluminance > solar_constant_l*1000.0) { - fprintf(stderr,"ERROR: direct illuminance exceeds solar constant\n"); - exit(1); + + if (directilluminance+diffuseilluminance==0 && altitude > 0) { + if(suppress_warnings==0) + { fprintf(stderr,"Warning: zero illuminance at sun altitude > 0, printing error sky\n"); } + print_error_sky(); + exit(0); } + + if (directilluminance > solar_constant_l) { + if(suppress_warnings==0) + { fprintf(stderr,"Warning: direct illuminance exceeds solar constant\n"); } + print_error_sky(); + exit(0); + } } void check_irradiances() { if (directirradiance < 0) { - fprintf(stderr,"WARNING: direct irradiance < 0. Using 0.0\n"); + if(suppress_warnings==0) + { fprintf(stderr,"Warning: direct irradiance < 0. Using 0.0\n"); } directirradiance = 0.0; } if (diffuseirradiance < 0) { - fprintf(stderr,"WARNING: diffuse irradiance < 0. Using 0.0\n"); + if(suppress_warnings==0) + { fprintf(stderr,"Warning: diffuse irradiance < 0. Using 0.0\n"); } diffuseirradiance = 0.0; } + + if (directirradiance+diffuseirradiance==0 && altitude > 0) { + if(suppress_warnings==0) + { fprintf(stderr,"Warning: zero irradiance at sun altitude > 0, printing error sky\n"); } + print_error_sky(); + exit(0); + } + if (directirradiance > solar_constant_e) { - fprintf(stderr,"ERROR: direct irradiance exceeds solar constant\n"); - exit(1); + if(suppress_warnings==0) + { fprintf(stderr,"Warning: direct irradiance exceeds solar constant\n"); } + print_error_sky(); + exit(0); } } @@ -1010,38 +1202,43 @@ double direct_irradiance_from_sky_clearness() } + + void illu_to_irra_index() { double test1=0.1, test2=0.1, d_eff; int counter=0; -diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l*1000); -directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000); +diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l); +directirradiance = directilluminance*solar_constant_e/(solar_constant_l); skyclearness = sky_clearness(); skybrightness = sky_brightness(); check_parametrization(); - + + while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10) - || skyclearness>skyclearinf || skyclearnessskybriginf || skybrightnessskyclearsup)) + || (!(skybrightnessskybrigsup)) ) + && !(counter==9) ) { - + test1=diffuseirradiance; test2=directirradiance; counter++; diffuseirradiance = diffuseilluminance/glob_h_diffuse_effi_PEREZ(); d_eff = direct_n_effi_PEREZ(); + + if (d_eff < 0.1) directirradiance = 0; - else + else directirradiance = directilluminance/d_eff; skybrightness = sky_brightness(); skyclearness = sky_clearness(); check_parametrization(); - + } @@ -1070,13 +1267,14 @@ static int get_numlin(float epsilon) /* sky luminance perez model */ double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]) { + float x[5][4]; int i,j,num_lin; double c_perez[5]; if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) ) { - fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez!\n"); + fprintf(stderr,"Error: epsilon out of range in function calc_rel_lum_perez!\n"); exit(1); } @@ -1085,14 +1283,15 @@ double calc_rel_lum_perez(double dzeta,double gamma,do { if ( Delta < 0.2 ) Delta = 0.2; } - + + num_lin = get_numlin(epsilon); - + for (i=0;i<5;i++) for (j=0;j<4;j++) { x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j); - /* printf("x %d %d vaut %f\n",i,j,x[i][j]); */ + /* fprintf(stderr,"x %d %d vaut %f\n",i,j,x[i][j]); */ } @@ -1126,7 +1325,7 @@ void coeff_lum_perez(double Z, double epsilon, double if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) ) { - fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n"); + fprintf(stderr,"Error: epsilon out of range in function coeff_lum_perez!\n"); exit(1); } @@ -1135,10 +1334,11 @@ void coeff_lum_perez(double Z, double epsilon, double { if ( Delta < 0.2 ) Delta = 0.2; } - + + num_lin = get_numlin(epsilon); - //fprintf(stderr,"numlin %d\n", num_lin); + /*fprintf(stderr,"numlin %d\n", num_lin);*/ for (i=0;i<5;i++) for (j=0;j<4;j++) @@ -1170,18 +1370,21 @@ void coeff_lum_perez(double Z, double epsilon, double } + /* degrees into radians */ double radians(double degres) { return degres*M_PI/180.0; } + /* radian into degrees */ double degres(double radians) { return radians/M_PI*180.0; } + /* calculation of the angles dzeta and gamma */ void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z) { @@ -1191,7 +1394,7 @@ void theta_phi_to_dzeta_gamma(double theta,double phi, else if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1.1 ) { printf("error in calculation of gamma (angle between point and sun"); - exit(3); + exit(1); } else *gamma = acos(cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)); @@ -1199,38 +1402,21 @@ void theta_phi_to_dzeta_gamma(double theta,double phi, -/********************************************************************************/ -/* Fonction: integ_lv */ -/* */ -/* In: float *lv,*theta */ -/* int sun_pos */ -/* */ -/* Out: double */ -/* */ -/* Update: 29/08/93 */ -/* */ -/* Rem: */ -/* */ -/* But: calcul l'integrale de luminance relative sans la dir. du soleil */ -/* */ -/********************************************************************************/ double integ_lv(float *lv,float *theta) { int i; double buffer=0.0; - + for (i=0;i<145;i++) + { buffer += (*(lv+i))*cos(radians(*(theta+i))); - + } + return buffer*2*M_PI/144; - } - - - /* enter day number(double), return E0 = square(R0/R): eccentricity correction factor */ double get_eccentricity() @@ -1243,7 +1429,6 @@ double get_eccentricity() 0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle); return (E0); - } @@ -1251,13 +1436,12 @@ double get_eccentricity() double air_mass() { double m; - if (sunzenith>90) { - fprintf(stderr, "Solar zenith angle larger than 90 degrees in function air_mass()\n"); - exit(1); + if(suppress_warnings==0) + { fprintf(stderr, "Warning: air mass has reached the maximal value\n"); } + sunzenith=90; } - m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) ); return(m); }