--- ray/src/gen/gensky.c	2004/09/10 18:19:24	2.23
+++ ray/src/gen/gensky.c	2020/07/25 19:18:01	2.28
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: gensky.c,v 2.23 2004/09/10 18:19:24 greg Exp $";
+static const char	RCSid[] = "$Id: gensky.c,v 2.28 2020/07/25 19:18:01 greg Exp $";
 #endif
 /*
  *  gensky.c - program to generate sky functions.
@@ -10,20 +10,13 @@ static const char	RCSid[] = "$Id: gensky.c,v 2.23 2004
  *     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];
@@ -96,13 +86,13 @@ void printdefaults(void);
 void userror(char  *msg);
 double normsc(void);
 int cvthour(char  *hs);
-void printhead(register int  ac, register char  **av);
 
 
 int
-main(argc, argv)
-int  argc;
-char  *argv[];
+main(
+	int  argc,
+	char  *argv[]
+)
 {
 	int  got_meridian = 0;
 	int  i;
@@ -134,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;
@@ -180,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();
@@ -200,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);
@@ -361,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);
@@ -380,8 +383,8 @@ 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++;
@@ -415,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');
 }