ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/sun.c
(Generate patch)

Comparing ray/src/gen/sun.c (file contents):
Revision 2.6 by greg, Fri Aug 9 16:51:15 2013 UTC vs.
Revision 2.7 by greg, Thu Nov 7 23:15:07 2019 UTC

# Line 6 | Line 6 | static const char      RCSid[] = "$Id$";
6   *
7   *               3/31/87
8   *
9 + *              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   */
16  
17   #include  <math.h>
# Line 14 | Line 20 | static const char      RCSid[] = "$Id$";
20   #ifdef M_PI
21   #define  PI     M_PI
22   #else
23 < #define  PI     3.141592653589793
23 > #define  PI     3.14159265358979323846
24   #endif
25 + #undef DEG
26 + #define DEG     (PI/180.)
27  
20 – double  s_latitude = 0.66;      /* site latitude (radians) */
21 – double  s_longitude = 2.13;     /* site longitude (radians) */
22 – double  s_meridian = 2.0944;    /* standard meridian (radians) */
28  
29 + 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  
33 +
34   int
35   jdate(          /* Julian date (days into year) */
36          int month,
# Line 39 | Line 48 | stadj(         /* solar time adjustment from Julian date */
48          int  jd
49   )
50   {
51 <        return( 0.170 * sin( (4*PI/373) * (jd - 80) ) -
52 <                0.129 * sin( (2*PI/355) * (jd - 8) ) +
53 <                12 * (s_meridian - s_longitude) / PI );
51 >        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   }
55  
56  
# Line 50 | Line 59 | sdec(          /* solar declination angle from Julian date */
59          int  jd
60   )
61   {
62 <        return( 0.4093 * sin( (2*PI/368) * (jd - 81) ) );
62 >        return( 0.4093 * sin( (2.*PI/365.) * (jd - 81) ) );
63   }
64  
65  
# Line 61 | Line 70 | salt(  /* solar altitude from solar declination and sol
70   )
71   {
72          return( asin( sin(s_latitude) * sin(sd) -
73 <                        cos(s_latitude) * cos(sd) * cos(st*(PI/12)) ) );
73 >                        cos(s_latitude) * cos(sd) * cos(st*(PI/12.)) ) );
74   }
75  
76  
# Line 71 | Line 80 | sazi(  /* solar azimuth from solar declination and sola
80          double st
81   )
82   {
83 <        return( -atan2( cos(sd)*sin(st*(PI/12)),
83 >        return( -atan2( cos(sd)*sin(st*(PI/12.)),
84                          -cos(s_latitude)*sin(sd) -
85 <                        sin(s_latitude)*cos(sd)*cos(st*(PI/12)) ) );
85 >                        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.460*DEG + 0.9856474*DEG*mjd, 2.*PI);
119 >        double  g = norm_cir(357.528*DEG + 0.9856003*DEG*mjd, 2.*PI);
120 >        double  l = L + 1.915*DEG*sin(g) + 0.020*DEG*sin(2.*g);
121 >        double  ep = 23.439*DEG - 4e-7*DEG*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   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines