src/gen/sun.c

#ifndef lint
static const char       RCSid[] = "$Id: sun.c,v 2.6 2013/08/09 16:51:15 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"

#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 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(          /* 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};
        
        return(mo_da[month-1] + day);
}


double
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./PI) * (s_meridian - s_longitude) );
}


double
sdec(           /* solar declination angle from Julian date */
        int  jd
)
{
        return( 0.4093 * sin( (2.*PI/365.) * (jd - 81) ) );
}


double
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.)) ) );
}


double
sazi(   /* solar azimuth from solar declination and solar time */
        double sd,
        double st
)
{
        return( -atan2( cos(sd)*sin(st*(PI/12.)),
                        -cos(s_latitude)*sin(sd) -
                        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 */
}
Revision:	2.7
Committed:	Thu Nov 7 23:15:07 2019 UTC (4 years, 5 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R4, rad5R3, HEAD
Changes since 2.6:	+68 -12 lines
Log Message:	Added more accurate Michalsky solar position calculation and made correction to old IES handbook formula (thanks to Axel Jacobs)
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.7	static const char RCSid[] = "$Id: sun.c,v 2.6 2013/08/09 16:51:15 greg Exp $";
3	greg	1.1	#endif
4	greg	1.2	/*
5	greg	1.1	* SOLAR CALCULATIONS
6			*
7			* 3/31/87
8			*
9	greg	2.7	* Michalsky algorithm added October 2019:
10			* "The Astronomical Almanac's Algorithm for Approximate
11			* Solar Position (1950-2050)" by Joseph J. Michalsky,
12			* published in 1988 in Solar Energy, Vol. 40, No. 3.
13			* Also added correction to sdec() (365 was 368 originally)
14			*
15	greg	1.1	*/
16
17	greg	2.2	#include <math.h>
18	greg	2.6	#include "sun.h"
19	greg	2.2
20	greg	2.5	#ifdef M_PI
21			#define PI M_PI
22			#else
23	greg	2.7	#define PI 3.14159265358979323846
24	greg	2.5	#endif
25	greg	2.7	#undef DEG
26			#define DEG (PI/180.)
27
28	greg	1.1
29	greg	2.7	double s_latitude = 0.66; /* site latitude (radians north of equator) */
30			double s_longitude = 2.13; /* site longitude (radians west of Greenwich) */
31			double s_meridian = 120.DEG; / standard meridian (radians west) */
32	greg	1.1
33
34			int
35	schorsch	2.4	jdate( /* Julian date (days into year) */
36			int month,
37			int day
38			)
39	greg	1.1	{
40			static short mo_da[12] = {0,31,59,90,120,151,181,212,243,273,304,334};
41
42			return(mo_da[month-1] + day);
43			}
44
45
46			double
47	schorsch	2.4	stadj( /* solar time adjustment from Julian date */
48			int jd
49			)
50	greg	1.1	{
51	greg	2.7	return( 0.170 * sin( (4.PI/373.) (jd - 80) ) -
52			0.129 * sin( (2.PI/355.) (jd - 8) ) +
53			(12./PI) * (s_meridian - s_longitude) );
54	greg	1.1	}
55
56
57			double
58	schorsch	2.4	sdec( /* solar declination angle from Julian date */
59			int jd
60			)
61	greg	1.1	{
62	greg	2.7	return( 0.4093 * sin( (2.PI/365.) (jd - 81) ) );
63	greg	1.1	}
64
65
66			double
67	schorsch	2.4	salt( /* solar altitude from solar declination and solar time */
68			double sd,
69			double st
70			)
71	greg	1.1	{
72			return( asin( sin(s_latitude) * sin(sd) -
73	greg	2.7	cos(s_latitude) * cos(sd) * cos(st*(PI/12.)) ) );
74	greg	1.1	}
75
76
77			double
78	schorsch	2.4	sazi( /* solar azimuth from solar declination and solar time */
79			double sd,
80			double st
81			)
82	greg	1.1	{
83	greg	2.7	return( -atan2( cos(sd)sin(st(PI/12.)),
84	greg	1.1	-cos(s_latitude)*sin(sd) -
85	greg	2.7	sin(s_latitude)cos(sd)cos(st*(PI/12.)) ) );
86			}
87
88
89			/**************** More accurate Michalsky algorithm **************/
90
91
92			/* circle normalization */
93			static double
94			norm_cir(double r, const double p)
95			{
96			while (r < 0) r += p;
97			while (r >= p) r -= p;
98			return(r);
99			}
100
101
102			/* Almanac Julian date relative to noon UT on Jan 1, 2000 (fractional days) */
103			double
104			mjdate(int year, int month, int day, double hour)
105			{
106			int jd = jdate(month, day);
107			jd += (month > 2) & !(year&3);
108			jd += (year - 1949)*365 + (year - 1949)/4;
109			hour += s_meridian*(12./PI);
110			return(jd + hour*(1./24.) + (2432916.5-2451545.));
111			}
112
113
114			/* Solar declination (and solar time) from Almanac Julian date (fractional) */
115			double
116			msdec(double mjd, double *stp)
117			{ /* Ecliptic coordinates (radians) */
118			double L = norm_cir(280.460DEG + 0.9856474DEGmjd, 2.PI);
119			double g = norm_cir(357.528DEG + 0.9856003DEGmjd, 2.PI);
120			double l = L + 1.915DEGsin(g) + 0.020DEGsin(2.*g);
121			double ep = 23.439DEG - 4e-7DEG*mjd;
122			double sin_l = sin(l);
123
124			if (stp) { /* solar time requested, also? */
125			double ra = atan2(sin_l*cos(ep), cos(l));
126			double utime = 24.*(mjd - floor(mjd)) + 12.;
127			double gmst = 6.697375 + 0.0657098242*mjd + utime;
128			double lmst = gmst - s_longitude*(12./PI);
129
130			stp = norm_cir(lmst - ra(12./PI) + 12., 24.);
131			}
132			return(asin(sin(ep)sin_l)); / return solar declination angle */
133	greg	1.1	}