--- ray/src/gen/gendaylit.c 2013/11/18 18:07:16 2.15 +++ ray/src/gen/gendaylit.c 2021/01/28 19:03:15 2.21 @@ -4,6 +4,8 @@ * Centre de Valbonne, 500 route des Lucioles, 06565 Sophia Antipolis Cedex, France * *BOUYGUES * 1 Avenue Eugene Freyssinet, Saint-Quentin-Yvelines, France +* print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018 +* version 2.6 (2021/01/29): dew point dependency added according to Perez publication 1990 (W -> atm_preci_water=exp(0.07*Td-0.075) ). by J. Wienold, EPFL */ #define _USE_MATH_DEFINES @@ -20,7 +22,7 @@ double normsc(); -/*static char *rcsid="$Header: /usr/local/cvs/radiance/ray/src/gen/gendaylit.c,v 2.15 2013/11/18 18:07:16 greg Exp $";*/ +/*static char *rcsid="$Header: /usr/local/cvs/radiance/ray/src/gen/gendaylit.c,v 2.21 2021/01/28 19:03:15 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, @@ -53,6 +55,9 @@ float defangle_phi[] = { 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}; +/* default values for Berlin */ +float locus[] = { +-4.843e9,2.5568e6,0.24282e3,0.23258,-4.843e9,2.5568e6,0.24282e3,0.23258,-1.2848,1.7519,-0.093786}; @@ -100,15 +105,7 @@ double air_mass(); 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; - const double AU = 149597890E3; const double solar_constant_e = 1367; /* solar constant W/m^2 */ const double solar_constant_l = 127500; /* solar constant lux */ @@ -124,6 +121,7 @@ const double skybrigsup = 0.6; /* required values */ +int year = 0; /* year (optional) */ int month, day; /* date */ double hour; /* time */ int tsolar; /* 0=standard, 1=solar */ @@ -136,7 +134,7 @@ double skyclearness = 0; double skybrightness = 0; double solarradiance; double diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance; -double sunzenith, daynumber, atm_preci_water=2; +double sunzenith, daynumber, atm_preci_water, Td=10.97353115; /*double sunaltitude_border = 0;*/ double diffnormalization = 0; @@ -146,7 +144,7 @@ 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 color_output=0; int suppress_warnings=0; /* default values */ @@ -157,6 +155,7 @@ double betaturbidity = 0.1; double gprefl = 0.2; int S_INTER=0; + /* computed values */ double sundir[3]; double groundbr = 0; @@ -201,10 +200,18 @@ int main(int argc, char** argv) for (i = 4; i < argc; i++) if (argv[i][0] == '-' || argv[i][0] == '+') switch (argv[i][1]) { + case 'd': + Td = atof(argv[++i]); + if (Td < -40 || Td > 40) { + Td=10.97353115; } + break; case 's': cloudy = 0; dosun = argv[i][0] == '+'; break; + case 'y': + year = atoi(argv[++i]); + break; case 'R': u_solar = argv[i][1] == 'R' ? -1 : 1; solarbr = atof(argv[++i]); @@ -213,6 +220,36 @@ int main(int argc, char** argv) cloudy = argv[i][0] == '+' ? 2 : 1; dosun = 0; break; + case 'C': + if (argv[i][2] == 'I' && argv[i][3] == 'E' ) { + locus[0] = -4.607e9; + locus[1] = 2.9678e6; + locus[2] = 0.09911e3; + locus[3] = 0.244063; + locus[4] = -2.0064e9; + locus[5] = 1.9018e6; + locus[6] = 0.24748e3; + locus[7] = 0.23704; + locus[8] = -3.0; + locus[9] = 2.87; + locus[10] = -0.275; + }else{ color_output = 1; + } + break; + case 'l': + locus[0] = atof(argv[++i]); + locus[1] = atof(argv[++i]); + locus[2] = atof(argv[++i]); + locus[3] = atof(argv[++i]); + locus[4] = locus[0]; + locus[5] = locus[1]; + locus[6] = locus[2]; + locus[7] = locus[3]; + locus[8] = atof(argv[++i]); + locus[9] = atof(argv[++i]); + locus[10] = atof(argv[++i]); + break; + case 't': betaturbidity = atof(argv[++i]); break; @@ -292,6 +329,7 @@ int main(int argc, char** argv) { fprintf(stderr,"Out of memory error in function main"); return 1; } + atm_preci_water=exp(0.07*Td-0.075); printhead(argc, argv); computesky(); printsky(); @@ -318,17 +356,22 @@ void computesky() /* compute solar direction */ if (month) { /* from date and time */ - int jd; double sd; - jd = jdate(month, day); /* Julian date */ - sd = sdec(jd); /* solar declination */ - if (tsolar) /* solar time */ - st = hour; - else - st = hour + stadj(jd); - - + st = hour; + if (year) { /* Michalsky algorithm? */ + double mjd = mjdate(year, month, day, hour); + if (tsolar) + sd = msdec(mjd, NULL); + else + sd = msdec(mjd, &st); + } else { + int jd = jdate(month, day); /* Julian date */ + sd = sdec(jd); /* solar declination */ + if (!tsolar) /* get solar time? */ + st = hour + stadj(jd); + } + if(timeinterval) { if(timeinterval<0) { @@ -664,6 +707,7 @@ 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); + printf("# epsilon, delta, atmospheric precipitable water content : %.4f %.4f %.4f \n", skyclearness, skybrightness,atm_preci_water ); if (dosun&&(skyclearness>1)) @@ -682,14 +726,31 @@ void printsky() printf("0\n0\n"); 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"); +/* print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018 */ + if (color_output==1 && skyclearness < 4.5 && skyclearness >1.065 ) + { + fprintf(stderr, " warning: sky clearness(epsilon)= %f \n",skyclearness); + fprintf(stderr, " warning: intermediate sky!! \n"); + fprintf(stderr, " warning: color model for intermediate sky pending \n"); + fprintf(stderr, " warning: no color output ! \n"); + color_output=0; + } + if (color_output==1) + { + printf("\nvoid colorfunc skyfunc\n"); + printf("4 skybright_r skybright_g skybright_b perezlum_c.cal\n"); printf("0\n"); - printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr, + printf("22 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n", diffnormalization, groundbr, *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4), - sundir[0], sundir[1], sundir[2]); + sundir[0], sundir[1], sundir[2],skyclearness,locus[0],locus[1],locus[2],locus[3],locus[4],locus[5],locus[6],locus[7],locus[8],locus[9],locus[10]); + }else{ + 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]); + } } @@ -735,12 +796,12 @@ void usage_error(char* msg) /* print usage error and { if (msg != NULL) 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, "Usage: %s month day hour [-y year] [...]\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, " 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"); @@ -748,7 +809,7 @@ void usage_error(char* msg) /* print usage error and 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.4 (2013/09/04) \n\n"); + fprintf(stderr, " gendaylit version 2.5 (2018/04/18) \n\n"); exit(1); } @@ -1374,14 +1435,14 @@ void coeff_lum_perez(double Z, double epsilon, double /* degrees into radians */ double radians(double degres) { - return degres*M_PI/180.0; + return degres*(M_PI/180.); } /* radian into degrees */ double degres(double radians) { - return radians/M_PI*180.0; + return radians*(180./M_PI); } @@ -1412,7 +1473,7 @@ double integ_lv(float *lv,float *theta) buffer += (*(lv+i))*cos(radians(*(theta+i))); } - return buffer*2*M_PI/144; + return buffer*(2.*M_PI/145.); }