/[ascend]/trunk/models/johnpye/grena/sunpos_grena.c
ViewVC logotype

Contents of /trunk/models/johnpye/grena/sunpos_grena.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2471 - (show annotations) (download) (as text)
Sat May 28 10:24:10 2011 UTC (11 years, 1 month ago) by jpye
File MIME type: text/x-csrc
File size: 5292 byte(s)
Adding initial bindings to sun position calculation, some debugging completed. Testing is next.
1 /* ASCEND modelling environment
2 Copyright (C) 2011 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 *//** @FILE
19 This file's code is based on the source code given in the publication
20 R Grena (2008), An algorithm for the computation of the solar position,
21 Solar Energy (82), pp 462-470.
22
23 The original code was in C++ and returned several intermediate results.
24 This modified version returns only the zenith and azimuth angles for
25 given date/time.
26 */
27
28 #include "sunpos_grena.h"
29
30 /* we have converted this file from C++ to C. not so classy any more ;-) */
31
32 void SunPos_calc_time(SunPos *S, double UT, int Day, int Month, int Year, double Delta_t){
33
34 // calculation of JD and JDE
35 double dYear, dMonth;
36 if(Month <= 2){
37 dYear = (double)Year - 1.;
38 dMonth = (double)Month + 12.;
39 }else{
40 dYear = (double)Year;
41 dMonth = (double)Month;
42 }
43
44 // universal time
45 S->t_G = (double)trunc(365.25 * (dYear - 2000))
46 + (double)trunc(30.6001 * (dMonth + 1))
47 + (double)Day + UT/24. - 1158.5;
48
49 S->Delta_t = Delta_t;
50 }
51
52
53 void SunPos_set_time(SunPos *S, double t_G, double Delta_t){
54 S->t_G = t_G;
55 S->Delta_t = Delta_t;
56 }
57
58
59 void SunPos_set_lat_long(SunPos *S, double latitude, double longitude){
60 S->latitude = latitude;
61 S->longitude = longitude;
62 }
63
64
65 void SunPos_set_pressure_temp(SunPos *S, double p, double T){
66 S->p = p;
67 S->T = T;
68 }
69
70 /*
71 Note: there is some confusion in nomenclature in the paper regarding
72 'JD_t' and 't_G' which are assumed to be the same time. This seems to John
73 to be correct.
74 */
75 void SunPos_calc_zen_azi(SunPos *S, double *zenith, double *azimuth){
76 double HourAngle;
77 double TopocRightAscension;
78 double TopocDeclination;
79 double TopocHourAngle;
80 double Elevation_no_refrac;
81 double RefractionCorrection;
82
83 // HELIOCENTRIC LONGITUDE
84
85 double t = S->t_G + S->Delta_t / 86400;
86
87 // linear increase + annual harmonic
88
89 double ang = 1.72019e-2 * t - 0.0563;
90 double HeliocLongitude = 1.740940 + 1.7202768683e-2 * t + 3.34118e-2 * sin(ang) + 3.488e-4 * sin(2*ang);
91
92 // moon perturbation
93
94 HeliocLongitude += 3.13e-5 * sin(2.127730e-1*t - 0.585);
95
96 // harmonic correction
97
98 HeliocLongitude += 1.26e-5 * sin(4.243e-3 * t + 1.46)
99 + 2.35e-5 * sin(1.0727e-2 * t + 0.72)
100 + 2.76e-5 * sin(1.5799e-2 * t + 2.35)
101 + 2.75e-5 * sin(2.1551e-2 * t - 1.98)
102 + 1.26e-5 * sin(3.1490e-2 * t - 0.80);
103
104 // polynomial correction
105
106 double t2 = t/1000;
107 HeliocLongitude += ((( -2.30796e-7 * t2 + 3.7976e-6) * t2 - 2.0458e-5) * t + 3.976e-5) * t2*t2;
108
109 // to obtain obtain Heliocentric longitude in the range [0,2pi] uncomment:
110 // HeliocLongitude = fmod(HeliocLongitude, 2*PI);
111
112 // END HELIOCENTRIC LONGITUDE CALCULATION
113
114 // Correction to longitude due to nutation
115
116 double delta_psi = 8.33e-5 * sin(9.252e-4 * t - 1.173);
117
118 // Earth axis inclination
119
120 double epsilon = -6.21e-9 * t + 0.409086 + 4.46e-5 * sin(9.252e-4 * t + 0.397);
121
122 // Geocentric global solar coordinates:
123
124 double GeocSolarLongitude = HeliocLongitude + PI + delta_psi - 9.932e-5;
125
126 double s_lambda = sin(GeocSolarLongitude);
127
128 double RightAscension = atan2(s_lambda * cos(epsilon), cos(GeocSolarLongitude));
129
130 double Declination = asin(sin(epsilon) * s_lambda);
131
132 // local hour angle of the sun
133
134 HourAngle = 6.30038809903 * S->t_G + 4.8824623 + delta_psi * 0.9174
135 + S->longitude - RightAscension;
136
137 // to obtain the local hour angle in the range [0,2pi] uncomment:
138 // HourAngle = fmod(HourAngle,2*PI);
139
140 double c_lat = cos(S->latitude);
141 double s_lat = sin(S->latitude);
142 double c_H = cos(HourAngle);
143 double s_H = sin(HourAngle);
144
145 // parallax correction to right ascension
146
147 double d_alpha = -4.26e-5 * c_lat * s_H;
148 TopocRightAscension = RightAscension + d_alpha;
149 TopocHourAngle = HourAngle - d_alpha;
150
151 // parallax correction to declination:
152
153 TopocDeclination = Declination - 4.26e-5 * (s_lat - Declination * c_lat);
154 double s_delta_corr = sin(TopocDeclination);
155 double c_delta_corr = cos(TopocDeclination);
156 double c_H_corr = c_H + d_alpha * s_H;
157 double s_H_corr = s_H - d_alpha * c_H;
158
159 // solar elevation angle, without refraction correction
160 Elevation_no_refrac = asin(s_lat * s_delta_corr + c_lat * c_delta_corr * c_H_corr);
161
162 // refraction correction: it is calculated only
163 // if Elevation_no_refract > elev_min
164
165 const double elev_min = -0.01;
166
167 if(Elevation_no_refrac > elev_min){
168 RefractionCorrection = 0.084217 * S->p / (273 + S->T)
169 / tan(Elevation_no_refrac + 0.0031376 / (Elevation_no_refrac + 0.089186));
170 }else{
171 RefractionCorrection = 0;
172 }
173
174 // local coordinates of the sun
175
176 *zenith = PI/2 - Elevation_no_refrac - RefractionCorrection;
177
178 *azimuth = atan2(s_H_corr, c_H_corr*s_lat - s_delta_corr/c_delta_corr*c_lat);
179 }
180

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22