--- ray/src/gen/sun.c	2003/02/22 02:07:24	2.3
+++ ray/src/gen/sun.c	2019/11/07 23:15:07	2.7
@@ -1,25 +1,41 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: sun.c,v 2.3 2003/02/22 02:07:24 greg Exp $";
+static const char	RCSid[] = "$Id: sun.c,v 2.7 2019/11/07 23:15:07 greg Exp $";
 #endif
 /*
  *           SOLAR CALCULATIONS
  *
  *               3/31/87
  *
+ *		Michalsky algorithm added October 2019:
+ *		"The Astronomical Almanac's Algorithm for Approximate
+ *		Solar Position (1950-2050)" by Joseph J. Michalsky,
+ *		published in 1988 in Solar Energy, Vol. 40, No. 3.
+ *		Also added correction to sdec() (365 was 368 originally)
+ *
  */
 
 #include  <math.h>
+#include  "sun.h"
 
-#define  PI  3.141592654
+#ifdef M_PI
+#define  PI	M_PI
+#else
+#define  PI	3.14159265358979323846
+#endif
+#undef DEG
+#define	DEG	(PI/180.)
 
-double  s_latitude = 0.66;	/* site latitude (radians) */
-double  s_longitude = 2.13;	/* site longitude (radians) */
-double  s_meridian = 2.0944;	/* standard meridian (radians) */
 
+double  s_latitude = 0.66;	/* site latitude (radians north of equator) */
+double  s_longitude = 2.13;	/* site longitude (radians west of Greenwich) */
+double  s_meridian = 120.*DEG;	/* standard meridian (radians west) */
 
+
 int
-jdate(month, day)		/* Julian date (days into year) */
-int  month, day;
+jdate(		/* Julian date (days into year) */
+	int month,
+	int day
+)
 {
 	static short  mo_da[12] = {0,31,59,90,120,151,181,212,243,273,304,334};
 	
@@ -28,37 +44,90 @@ int  month, day;
 
 
 double
-stadj(jd)		/* solar time adjustment from Julian date */
-int  jd;
+stadj(		/* solar time adjustment from Julian date */
+	int  jd
+)
 {
-	return( 0.170 * sin( (4*PI/373) * (jd - 80) ) -
-		0.129 * sin( (2*PI/355) * (jd - 8) ) +
-		12 * (s_meridian - s_longitude) / PI );
+	return( 0.170 * sin( (4.*PI/373.) * (jd - 80) ) -
+		0.129 * sin( (2.*PI/355.) * (jd - 8) ) +
+		(12./PI) * (s_meridian - s_longitude) );
 }
 
 
 double
-sdec(jd)		/* solar declination angle from Julian date */
-int  jd;
+sdec(		/* solar declination angle from Julian date */
+	int  jd
+)
 {
-	return( 0.4093 * sin( (2*PI/368) * (jd - 81) ) );
+	return( 0.4093 * sin( (2.*PI/365.) * (jd - 81) ) );
 }
 
 
 double
-salt(sd, st)	/* solar altitude from solar declination and solar time */
-double  sd, st;
+salt(	/* solar altitude from solar declination and solar time */
+	double sd,
+	double st
+)
 {
 	return( asin( sin(s_latitude) * sin(sd) -
-			cos(s_latitude) * cos(sd) * cos(st*(PI/12)) ) );
+			cos(s_latitude) * cos(sd) * cos(st*(PI/12.)) ) );
 }
 
 
 double
-sazi(sd, st)	/* solar azimuth from solar declination and solar time */
-double  sd, st;
+sazi(	/* solar azimuth from solar declination and solar time */
+	double sd,
+	double st
+)
 {
-	return( -atan2( cos(sd)*sin(st*(PI/12)),
+	return( -atan2( cos(sd)*sin(st*(PI/12.)),
  			-cos(s_latitude)*sin(sd) -
- 			sin(s_latitude)*cos(sd)*cos(st*(PI/12)) ) );
+ 			sin(s_latitude)*cos(sd)*cos(st*(PI/12.)) ) );
+}
+
+
+/****************** More accurate Michalsky algorithm ****************/
+
+
+/* circle normalization */
+static double
+norm_cir(double r, const double p)
+{
+	while (r < 0) r += p;
+	while (r >= p) r -= p;
+	return(r);
+}
+
+
+/* Almanac Julian date relative to noon UT on Jan 1, 2000 (fractional days) */
+double
+mjdate(int year, int month, int day, double hour)
+{
+	int	jd = jdate(month, day);
+	jd += (month > 2) & !(year&3);
+	jd += (year - 1949)*365 + (year - 1949)/4;
+	hour += s_meridian*(12./PI);
+	return(jd + hour*(1./24.) + (2432916.5-2451545.));
+}
+
+
+/* Solar declination (and solar time) from Almanac Julian date (fractional) */
+double
+msdec(double mjd, double *stp)
+{					/* Ecliptic coordinates (radians) */
+	double	L = norm_cir(280.460*DEG + 0.9856474*DEG*mjd, 2.*PI);
+	double	g = norm_cir(357.528*DEG + 0.9856003*DEG*mjd, 2.*PI);
+	double	l = L + 1.915*DEG*sin(g) + 0.020*DEG*sin(2.*g);
+	double	ep = 23.439*DEG - 4e-7*DEG*mjd;
+	double	sin_l = sin(l);
+
+	if (stp) {			/* solar time requested, also? */
+		double	ra = atan2(sin_l*cos(ep), cos(l));
+		double	utime = 24.*(mjd - floor(mjd)) + 12.;
+		double	gmst = 6.697375 + 0.0657098242*mjd + utime;
+		double	lmst = gmst - s_longitude*(12./PI);
+
+		*stp = norm_cir(lmst - ra*(12./PI) + 12., 24.);
+	}
+	return(asin(sin(ep)*sin_l));	/* return solar declination angle */
 }