--- ray/src/gen/gensky.c	2003/11/16 10:29:38	2.22
+++ ray/src/gen/gensky.c	2023/03/30 20:19:05	2.29
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: gensky.c,v 2.22 2003/11/16 10:29:38 schorsch Exp $";
+static const char	RCSid[] = "$Id: gensky.c,v 2.29 2023/03/30 20:19:05 greg Exp $";
 #endif
 /*
  *  gensky.c - program to generate sky functions.
@@ -10,20 +10,13 @@ static const char	RCSid[] = "$Id: gensky.c,v 2.22 2003
  *     3/26/86
  */
 
-#include  <stdio.h>
+#include  "rtio.h"
 #include  <stdlib.h>
-#include  <string.h>
 #include  <math.h>
 #include  <ctype.h>
-
+#include  "sun.h"
 #include  "color.h"
 
-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);
-
 #ifndef  PI
 #define  PI		3.14159265358979323846
 #endif
@@ -38,10 +31,6 @@ extern double sazi(double sd, double st);
 #define  overcast	((skytype==S_OVER)|(skytype==S_UNIF))
 
 double  normsc();
-					/* sun calculation constants */
-extern double  s_latitude;
-extern double  s_longitude;
-extern double  s_meridian;
 
 #undef  toupper
 #define  toupper(c)	((c) & ~0x20)	/* ASCII trick to convert case */
@@ -69,6 +58,7 @@ struct {
 	{"", 0}
 };
 					/* required values */
+int  year = 0;					/* year (optional) */
 int  month, day;				/* date */
 double  hour;					/* time */
 int  tsolar;					/* 0=standard, 1=solar */
@@ -78,7 +68,7 @@ int  skytype = S_CLEAR;				/* sky type */
 int  dosun = 1;
 double  zenithbr = 0.0;
 int	u_zenith = 0;				/* -1=irradiance, 1=radiance */
-double  turbidity = 2.75;
+double  turbidity = 2.45;
 double  gprefl = 0.2;
 					/* computed values */
 double  sundir[3];
@@ -95,15 +85,16 @@ void printsky(void);
 void printdefaults(void);
 void userror(char  *msg);
 double normsc(void);
-void cvthour(char  *hs);
-void printhead(register int  ac, register char  **av);
+int cvthour(char  *hs);
 
 
 int
-main(argc, argv)
-int  argc;
-char  *argv[];
+main(
+	int  argc,
+	char  *argv[]
+)
 {
+	int  got_meridian = 0;
 	int  i;
 
 	progname = argv[0];
@@ -124,7 +115,7 @@ char  *argv[];
 		day = atoi(argv[2]);
 		if (day < 1 || day > 31)
 			userror("bad day");
-		cvthour(argv[3]);
+		got_meridian = cvthour(argv[3]);
 	}
 	for (i = 4; i < argc; i++)
 		if (argv[i][0] == '-' || argv[i][0] == '+')
@@ -133,6 +124,9 @@ char  *argv[];
 				skytype = S_CLEAR;
 				dosun = argv[i][0] == '+';
 				break;
+			case 'y':
+				year = atoi(argv[++i]);
+				break;
 			case 'r':
 			case 'R':
 				u_solar = argv[i][1]=='R' ? -1 : 1;
@@ -166,6 +160,10 @@ char  *argv[];
 				s_longitude = atof(argv[++i]) * (PI/180);
 				break;
 			case 'm':
+				if (got_meridian) {
+					++i;
+					break;		/* time overrides */
+				}
 				s_meridian = atof(argv[++i]) * (PI/180);
 				break;
 			default:
@@ -175,12 +173,17 @@ char  *argv[];
 		else
 			userror("bad option");
 
-	if (fabs(s_meridian-s_longitude) > 45*PI/180)
+	if (year && (year < 1950) | (year > 2050))
 		fprintf(stderr,
-	"%s: warning: %.1f hours btwn. standard meridian and longitude\n",
+			"%s: warning - year should be in range 1950-2050\n",
+			progname);
+	if (month && !tsolar && fabs(s_meridian-s_longitude) > 45*PI/180)
+		fprintf(stderr,
+	"%s: warning - %.1f hours btwn. standard meridian and longitude\n",
 			progname, (s_longitude-s_meridian)*12/PI);
 
-	printhead(argc, argv);
+	fputs("# ", stdout);
+	printargs(argc, argv, stdout);
 
 	computesky();
 	printsky();
@@ -195,15 +198,20 @@ computesky(void)			/* compute sky parameters */
 	double	normfactor;
 					/* compute solar direction */
 	if (month) {			/* from date and time */
-		int  jd;
-		double  sd, st;
+		double  sd, st = hour;
 
-		jd = jdate(month, day);		/* Julian date */
-		sd = sdec(jd);			/* solar declination */
-		if (tsolar)			/* solar time */
-			st = hour;
-		else
-			st = hour + stadj(jd);
+		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);
+		}
 		altitude = salt(sd, st);
 		azimuth = sazi(sd, st);
 		printf("# Local solar time: %.2f\n", st);
@@ -281,7 +289,7 @@ printsky(void)			/* print out sky */
 	if (dosun) {
 		printf("\nvoid light solar\n");
 		printf("0\n0\n");
-		printf("3 %.2e %.2e %.2e\n", solarbr, solarbr, solarbr);
+		printf("3 %.3e %.3e %.3e\n", solarbr, solarbr, solarbr);
 		printf("\nsolar source sun\n");
 		printf("0\n0\n");
 		printf("4 %f %f %f 0.5\n", sundir[0], sundir[1], sundir[2]);
@@ -291,9 +299,9 @@ printsky(void)			/* print out sky */
 	printf("2 skybr skybright.cal\n");
 	printf("0\n");
 	if (overcast)
-		printf("3 %d %.2e %.2e\n", skytype, zenithbr, groundbr);
+		printf("3 %d %.3e %.3e\n", skytype, zenithbr, groundbr);
 	else
-		printf("7 %d %.2e %.2e %.2e %f %f %f\n",
+		printf("7 %d %.3e %.3e %.3e %f %f %f\n",
 				skytype, zenithbr, groundbr, F2,
 				sundir[0], sundir[1], sundir[2]);
 }
@@ -356,9 +364,9 @@ normsc(void)			/* compute normalization factor (E0*F2/
 				/* intermediate sky approx. */
 		{3.5556, -2.7152, -1.3081, 1.0660, 0.60227},
 	};
-	register double  *nf;
+	double  *nf;
 	double  x, nsc;
-	register int  i;
+	int  i;
 					/* polynomial approximation */
 	nf = nfc[skytype==S_INTER];
 	x = (altitude - PI/4.0)/(PI/4.0);
@@ -370,13 +378,13 @@ normsc(void)			/* compute normalization factor (E0*F2/
 }
 
 
-void
+int
 cvthour(			/* convert hour string */
 	char  *hs
 )
 {
-	register char  *cp = hs;
-	register int	i, j;
+	char  *cp = hs;
+	int	i, j;
 
 	if ( (tsolar = *cp == '+') ) cp++;		/* solar time? */
 	while (isdigit(*cp)) cp++;
@@ -388,7 +396,7 @@ cvthour(			/* convert hour string */
 	}
 	while (isdigit(*cp)) cp++;
 	if (!*cp)
-		return;
+		return(0);
 	if (tsolar || !isalpha(*cp)) {
 		fprintf(stderr, "%s: bad time format: %s\n", progname, hs);
 		exit(1);
@@ -400,7 +408,7 @@ cvthour(			/* convert hour string */
 				break;
 		if (!cp[j] && !tzone[i].zname[j]) {
 			s_meridian = tzone[i].zmer * (PI/180);
-			return;
+			return(1);
 		}
 	} while (tzone[i++].zname[0]);
 
@@ -410,19 +418,4 @@ cvthour(			/* convert hour string */
 		fprintf(stderr, " %s", tzone[i].zname);
 	putc('\n', stderr);
 	exit(1);
-}
-
-
-void
-printhead(		/* print command header */
-	register int  ac,
-	register char  **av
-)
-{
-	putchar('#');
-	while (ac--) {
-		putchar(' ');
-		fputs(*av++, stdout);
-	}
-	putchar('\n');
 }