| 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 |
|
Parent Directory
| 
Revision Log