--- ray/src/gen/gendaylit.c	2018/08/31 16:01:45	2.17
+++ ray/src/gen/gendaylit.c	2021/01/28 19:03:15	2.21
@@ -5,6 +5,7 @@
  * 				*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
@@ -21,7 +22,7 @@
 
 double  normsc();
 
-/*static	char *rcsid="$Header: /usr/local/cvs/radiance/ray/src/gen/gendaylit.c,v 2.17 2018/08/31 16:01:45 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,
@@ -104,8 +105,6 @@ 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);
 
 const double	AU = 149597890E3;
 const double 	solar_constant_e = 1367;    /* solar constant W/m^2 */
@@ -122,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 */
@@ -134,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;
@@ -200,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]);
@@ -321,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();
@@ -347,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) {
@@ -693,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)) 
@@ -781,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");
@@ -1420,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);
 }
 
 
@@ -1458,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.);
 }