/[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 2470 - (show annotations) (download) (as text)
Sat May 28 09:12:38 2011 UTC (11 years, 1 month ago) by jpye
File MIME type: text/x-csrc
File size: 5050 byte(s)
Working on conversion of Grena code to C.
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 #include <math.h>
30
31 /* we have converted this file from C++ to C. not so classy any more ;-) */
32
33
34 double SunPos_calc_time(SunPos *S, double UT, int Day, int Month, int Year, double Delta_t){
35
36 // calculation of JD and JDE
37 double dYear, dMonth;
38 if(sunpos->Month <= 2){
39 dYear = (double)S->Year - 1.;
40 dMonth = (double)S->Month + 12.;
41 }else{
42 dYear = (double)S->Year;
43 dMonth = (double)S->Month;
44 }
45
46 double JD_t = (double)trunc(365.25 * (dYear - 2000))
47 + (double)trunc(30.6001 * (dMonth + 1))
48 + (double)S->Day + S->UT/24. - 1158.5;
49
50 double t = JD_t + sunpos->Delta_t/86400;
51
52 S->t = t;
53 }
54
55
56 void SunPos_set_lat_long(SunPos *S, double latitude, double longitude){
57 S->latitude = latitude;
58 S->longitude = longitude;
59 }
60
61
62 void SunPos_set_pressure_temp(SunPos *S, double p, double T){
63 S->p = p;
64 S->T = T;
65 }
66
67 void SunPos_set_time(SunPos *S, double t){
68 S->t = t;
69 }
70
71
72 void SunPos_calc_zen_azi(SunPos *S, double *zenith, double *azimuth){
73 double t = S->t;
74 double HourAngle;
75 double TopocRightAscension;
76 double TopocDeclination;
77 double TopocHourAngle;
78 double Elevation_no_refrac;
79 double RefractionCorrection;
80
81 // HELIOCENTRIC LONGITUDE
82
83 // linear increase + annual harmonic
84
85 double ang = 1.72019e-2 * t - 0.0563;
86 HeliocLongitude = 1.740940 + 1.7202768683e-2 * t + 3.34118e-2 * sin(ang) + 3.488e-4 * sin(2*ang);
87
88 // moon perturbation
89
90 HeliocLongitude += 3.13e-5 * sin(2.127730e-1*t - 0.585);
91
92 // harmonic correction
93
94 HeliocLongitude += 1.26e-5 * sin(4.243e-3 * t + 1.46)
95 + 2.35e-5 * sin(1.0727e-2 * t + 0.72)
96 + 2.76e-5 * sin(1.5799e-2 * t + 2.35)
97 + 2.75e-5 * sin(2.1551e-2 * t - 1.98)
98 + 1.26e-5 * sin(3.1490e-2 * t - 0.80);
99
100 // polynomial correction
101
102 double t2 = t/1000;
103 HeliocLongitude += ((( -2.30796e-7 * t2 + 3.7976e-6) * t2 - 2.0458e-5) * t + 3.976e-5) * t2*t2;
104
105 // to obtain obtain Heliocentric longitude in the range [0,2pi] uncomment:
106 // HeliocLongitude = fmod(HeliocLongitude, 2*PI);
107
108 // END HELIOCENTRIC LONGITUDE CALCULATION
109
110 // Correction to longitude due to nutation
111
112 double delta_psi = 8.33e-5 * sin(9.252e-4 * t - 1.173);
113
114 // Earth axis inclination
115
116 double epsilon = -6.21e-9 * t + 0.409086 + 4.46e-5 * sin(9.252e-4 * t + 0.397);
117
118 // Geocentric global solar coordinates:
119
120 GeocSolarLongitude = HeliocLongitude + PI + delta_psi - 9.932e-5;
121
122 double s_lambda = sin(GeocSolarLongitude);
123
124 RightAscension = atan2(s_lambda * cos(epsilon), cos(GeocSolarLongitude));
125
126 // local hour angle of the sun
127
128 HourAngle = 6.30038809903 * JD_t + 4.8824623 + delta_psi * 0.9174
129 + S->longitude - RightAscension;
130
131 // to obtain the local hour angle in the range [0,2pi] uncomment:
132 // HourAngle = fmod(HourAngle,2*PI);
133
134 double c_lat = cos(S->latitude);
135 double s_lat = sin(S->latitude);
136 double c_H = cos(HourAngle);
137 double s_H = sin(HourAngle);
138
139 // parallax correction to right ascension
140
141 double d_alpha = -4.26e-5 * c_lat * s_H;
142 TopocRightAscension = RightAscension + d_alpha;
143 TopocHourAngle = HourAngle - d_alpha;
144
145 // parallax correction to declination:
146
147 TopocDeclination = Declination - 4.26e-5 * (s_lat - Declination * c_lat);
148 double s_delta_corr = sin(TopocDeclination);
149 double c_delta_corr = cos(TopocDeclination);
150 double c_H_corr = c_H + d_alpha * s_H;
151 double s_H_corr = s_H - d_alpha * c_H;
152
153 // solar elevation angle, without refraction correction
154 Elevation_no_refrac = asin(s_lat * s_delta_corr + c_lat * c_delta_corr * c_H_corr);
155
156 // refraction correction: it is calculated only
157 // if Elevation_no_refract > elev_min
158
159 const double elev_min = -0.01;
160
161 if(Elevation_no_refrac > elev_min){
162 RefractionCorrection = 0.084217 * Pressure / (273 + Temperature)
163 / tan(Elevation_no_refrac + 0.0031376 / (Elevation_no_refrac + 0.089186));
164 }else{
165 RefractionCorrection = 0;
166
167 // local coordinates of the sun
168
169 *zenith = PI/2 - Elevation_no_refrac - RefractionCorrection;
170
171 *azimuth = atan2(s_H_corr, c_H_corr*s_lat - s_delta_corr/c_delta_corr*c_lat);
172 }
173

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