--- ray/src/gen/gendaylit.c 2013/08/09 16:44:19 2.12 +++ ray/src/gen/gendaylit.c 2018/08/31 16:01:45 2.17 @@ -4,8 +4,10 @@ * 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 */ +#define _USE_MATH_DEFINES #include #include #include @@ -19,7 +21,7 @@ double normsc(); -/*static char *rcsid="$Header: /usr/local/cvs/radiance/ray/src/gen/gendaylit.c,v 2.12 2013/08/09 16:44:19 greg Exp $";*/ +/*static char *rcsid="$Header: /usr/local/cvs/radiance/ray/src/gen/gendaylit.c,v 2.17 2018/08/31 16:01:45 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, @@ -52,10 +54,13 @@ 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}; -/* 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(); @@ -88,7 +93,7 @@ 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); @@ -97,27 +102,20 @@ FILE * frlibopen(char* fname); double get_eccentricity(); double air_mass(); -double solar_sunset(); -double solar_sunrise(); +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 = 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.0; /* 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; @@ -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; @@ -168,9 +167,9 @@ float timeinterval = 0; char *progname; char errmsg[128]; +double st; - int main(int argc, char** argv) { int i; @@ -181,7 +180,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); @@ -189,13 +188,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++) @@ -213,6 +212,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; @@ -236,8 +265,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': @@ -269,12 +299,6 @@ int main(int argc, char** argv) globalirradiance = atof(argv[++i]); break; - /* - case 'l': - sunaltitude_border = atof(argv[++i]); - break; - */ - case 'i': timeinterval = atof(argv[++i]); break; @@ -282,46 +306,32 @@ int main(int argc, char** argv) 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); /* 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(); - - if(*(c_perez+1)>0) - { - fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1)); - print_error_sky(); - exit(1); - } - printsky(); return 0; + } - - - void computesky() { @@ -338,7 +348,7 @@ void computesky() if (month) { /* from date and time */ int jd; - double sd, st; + double sd; jd = jdate(month, day); /* Julian date */ sd = sdec(jd); /* solar declination */ @@ -346,13 +356,7 @@ void computesky() st = hour; else st = hour + stadj(jd); - - if(stsolar_sunset(month,day)) { - print_error_sky(); - exit(1); - } - if(timeinterval) { @@ -360,19 +364,38 @@ void computesky() fprintf(stderr, "time interval negative\n"); exit(1); } - - if(fabs(solar_sunrise(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); @@ -384,18 +407,6 @@ void computesky() - /* if loop for the -l option. W.Sprenger (01/2013) */ - /* - 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(1); - } - */ - if (!cloudy && altitude > 87.*M_PI/180.) { @@ -407,6 +418,8 @@ void computesky() altitude = 87.*M_PI/180.; } + + sundir[0] = -sin(azimuth)*cos(altitude); sundir[1] = -cos(azimuth)*cos(altitude); sundir[2] = sin(altitude); @@ -414,8 +427,7 @@ void computesky() /* 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) /* P */ @@ -439,8 +451,9 @@ void computesky() 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*/ @@ -464,7 +477,7 @@ void computesky() 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 { @@ -504,7 +517,6 @@ void computesky() { if (suppress_warnings==0) fprintf(stderr, "Warning: global irradiance is higher than the time-dependent solar constant s0\n"); - globalirradiance=erbs_s0*0.999; } @@ -517,7 +529,7 @@ void computesky() 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!"); + printf("# WARNING: the -E option is only recommended for a rough estimation!\n"); directirradiance=directirradiance/sin(altitude); @@ -540,7 +552,7 @@ void computesky() - else {fprintf(stderr,"error at the input arguments"); exit(1);} + else { fprintf(stderr,"error at the input arguments"); exit(1); } @@ -562,8 +574,6 @@ void computesky() - - /*calculation of the modelled luminance */ for (j=0;j<145;j++) { @@ -609,62 +619,53 @@ void computesky() 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() -{ - 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]); -} - - - - - - double solar_sunset(int month,int day) { float W; @@ -674,6 +675,8 @@ double solar_sunset(int month,int day) } + + double solar_sunrise(int month,int day) { float W; @@ -685,15 +688,13 @@ double solar_sunrise(int month,int day) +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); - - - - - -void printsky() /* print out sky */ -{ if (dosun&&(skyclearness>1)) { printf("\nvoid light solar\n"); @@ -710,18 +711,57 @@ void printsky() /* print out sky */ printf("0\n0\n"); printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle); } +/* 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("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],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]); + } +} + + +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 %.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]); - + 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); @@ -735,7 +775,9 @@ 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\n", progname, msg); @@ -752,12 +794,13 @@ void userror(char* msg) /* print usage error and qui 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.3 (2013/08/08) \n\n"); + fprintf(stderr, " gendaylit version 2.5 (2018/04/18) \n\n"); exit(1); } + double normsc() /* compute normalization factor (E0*F2/L0) */ { static double nfc[2][5] = { @@ -781,6 +824,8 @@ double normsc() /* compute normalization fac + + void printhead(int ac, char** av) /* print command header */ { putchar('#'); @@ -794,6 +839,8 @@ void printhead(int ac, char** av) /* print command he + + /* Perez models */ /* Perez global horizontal luminous efficacy model */ @@ -865,7 +912,6 @@ double glob_h_effi_PEREZ() - for (i=1; i<=category_total_number; i++) { if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) ) @@ -879,6 +925,8 @@ double glob_h_effi_PEREZ() } + + /* global horizontal diffuse efficacy model, according to PEREZ */ double glob_h_diffuse_effi_PEREZ() { @@ -886,17 +934,12 @@ 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_PEREZ \n"); */ + 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; @@ -952,7 +995,6 @@ double glob_h_diffuse_effi_PEREZ() - category_number = -1; for (i=1; i<=category_total_number; i++) { @@ -962,9 +1004,9 @@ double glob_h_diffuse_effi_PEREZ() if (category_number == -1) { if (suppress_warnings==0) - fprintf(stderr, "ERROR: Model parameters out of range, skyclearness = %lf \n", skyclearness); + fprintf(stderr, "Warning: sky clearness (= %.3f) too high, printing error sky\n", skyclearness); print_error_sky(); - exit(1); + exit(0); } @@ -977,6 +1019,9 @@ double glob_h_diffuse_effi_PEREZ() + + + /* direct normal efficacy model, according to PEREZ */ double direct_n_effi_PEREZ() @@ -987,8 +1032,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) */ @@ -1062,28 +1107,30 @@ 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){ - if (suppress_warnings==0) - /* fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */ - skyclearness=skyclearsup-0.1; + /* 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); */ + /* if (suppress_warnings==0) + fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */ skybrightness=skybrigsup; } @@ -1092,37 +1139,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); } } @@ -1181,8 +1255,8 @@ 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(); @@ -1210,8 +1284,6 @@ while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(d skybrightness = sky_brightness(); skyclearness = sky_clearness(); check_parametrization(); - - /*fprintf(stderr,"skyclearness = %lf, skybrightness = %lf, directirradiance = %lf, diffuseirradiance = %lf\n",skyclearness, skybrightness, directirradiance, diffuseirradiance);*/ } @@ -1248,7 +1320,7 @@ double calc_rel_lum_perez(double dzeta,double gamma,do 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); } @@ -1299,7 +1371,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); } @@ -1368,7 +1440,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)); @@ -1412,8 +1484,9 @@ 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);