Arcadie Cracan / Helios

Sun position calculation library. Adaptation of Hannes Hassler's Helios class.

Dependents:   sunTracker weather_station_proj weather_station_project weather_station_proj_v1_2

Helios.cpp@0:ad31da30ae64, 2018-06-24 (annotated)

Committer:
acracan
Date:
Sun Jun 24 12:10:33 2018 +0000
Revision:
0:ad31da30ae64
Adaptation of Hannes Hassler's Helios class

Who changed what in which revision?

UserRevisionLine numberNew contents of line
acracan 0:ad31da30ae64 1 /*
acracan 0:ad31da30ae64 2 Helios.cpp-
acracan 0:ad31da30ae64 3 Copyright (c) 2011 Hannes Hassler. All rights reserved.
acracan 0:ad31da30ae64 4
acracan 0:ad31da30ae64 5 This library is free software; you can redistribute it and/or
acracan 0:ad31da30ae64 6 modify it under the terms of the GNU Lesser General Public
acracan 0:ad31da30ae64 7 License as published by the Free Software Foundation; either
acracan 0:ad31da30ae64 8 version 2.1 of the License, or (at your option) any later version.
acracan 0:ad31da30ae64 9
acracan 0:ad31da30ae64 10 This library is distributed in the hope that it will be useful,
acracan 0:ad31da30ae64 11 but WITHOUT ANY WARRANTY; without even the implied warranty of
acracan 0:ad31da30ae64 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
acracan 0:ad31da30ae64 13 Lesser General Public License for more details.
acracan 0:ad31da30ae64 14
acracan 0:ad31da30ae64 15 You should have received a copy of the GNU Lesser General Public
acracan 0:ad31da30ae64 16 License along with this library; if not, write to the Free Software
acracan 0:ad31da30ae64 17 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
acracan 0:ad31da30ae64 18
acracan 0:ad31da30ae64 19 This library can be used for calculating the solar position on Arduino.
acracan 0:ad31da30ae64 20
acracan 0:ad31da30ae64 21 The algorithm is an adaption from
acracan 0:ad31da30ae64 22 the "PSA" solar positioning algorithm, as documented in:
acracan 0:ad31da30ae64 23
acracan 0:ad31da30ae64 24 Blanco-Muriel et al.: Computing the Solar Vector. Solar Energy Vol 70 No 5 pp 431-441.
acracan 0:ad31da30ae64 25 http://dx.doi.org/10.1016/S0038-092X(00)00156-0
acracan 0:ad31da30ae64 26
acracan 0:ad31da30ae64 27 According to the paper, "The algorithm allows .. the true solar vector
acracan 0:ad31da30ae64 28 to be determined with an accuracy of 0.5
acracan 0:ad31da30ae64 29 minutes of arc for the period 1999–2015.
acracan 0:ad31da30ae64 30 The original code has been downloaded from
acracan 0:ad31da30ae64 31 http://www.psa.es/sdg/sunpos.htm
acracan 0:ad31da30ae64 32
acracan 0:ad31da30ae64 33 Adaptions:
acracan 0:ad31da30ae64 34 Modified calculation of number of Days since 1.Jan 2000 (dJulianDate-2451545.0)
acracan 0:ad31da30ae64 35 Neccessary because of the limited double precision on Arduino
acracan 0:ad31da30ae64 36 (double has the same precision as float on the current Arduino (2011))
acracan 0:ad31da30ae64 37 It should be used only for dates between 1.1.2000 and 31.12.2100
acracan 0:ad31da30ae64 38 (PSA itself has garantueed accuracy only until 2015)
acracan 0:ad31da30ae64 39
acracan 0:ad31da30ae64 40 */
acracan 0:ad31da30ae64 41
acracan 0:ad31da30ae64 42 #include "Helios.h"
acracan 0:ad31da30ae64 43 #include <math.h>
acracan 0:ad31da30ae64 44 #include <mbed.h>
acracan 0:ad31da30ae64 45
acracan 0:ad31da30ae64 46
acracan 0:ad31da30ae64 47 Helios::Helios(double latitude, double longitude, int tzOffset)
acracan 0:ad31da30ae64 48 {
acracan 0:ad31da30ae64 49 setLocalLatitude(latitude);
acracan 0:ad31da30ae64 50 setLocalLongitude(longitude);
acracan 0:ad31da30ae64 51 setLocalTimeZoneOffset(tzOffset);
acracan 0:ad31da30ae64 52 }
acracan 0:ad31da30ae64 53
acracan 0:ad31da30ae64 54
acracan 0:ad31da30ae64 55 void Helios::updatePosition()
acracan 0:ad31da30ae64 56 {
acracan 0:ad31da30ae64 57 // rtc time
acracan 0:ad31da30ae64 58 time_t rtcTime;
acracan 0:ad31da30ae64 59 struct tm *rtcTimeInfo;
acracan 0:ad31da30ae64 60
acracan 0:ad31da30ae64 61 double dElapsedJulianDays;
acracan 0:ad31da30ae64 62 double dDecimalHours;
acracan 0:ad31da30ae64 63 double dEclipticLongitude;
acracan 0:ad31da30ae64 64 double dEclipticObliquity;
acracan 0:ad31da30ae64 65 double dRightAscension;
acracan 0:ad31da30ae64 66 double dDeclination;
acracan 0:ad31da30ae64 67
acracan 0:ad31da30ae64 68 double dSin_EclipticLongitude;
acracan 0:ad31da30ae64 69 double dMeanLongitude;
acracan 0:ad31da30ae64 70 double dMeanAnomaly;
acracan 0:ad31da30ae64 71 double dOmega;
acracan 0:ad31da30ae64 72
acracan 0:ad31da30ae64 73 double dGreenwichMeanSiderealTime;
acracan 0:ad31da30ae64 74 double dLocalMeanSiderealTime;
acracan 0:ad31da30ae64 75
acracan 0:ad31da30ae64 76 double dHourAngle;
acracan 0:ad31da30ae64 77 double dCos_HourAngle;
acracan 0:ad31da30ae64 78 double dParallax;
acracan 0:ad31da30ae64 79 double dZenithAngle;
acracan 0:ad31da30ae64 80
acracan 0:ad31da30ae64 81 // get current time from RTC
acracan 0:ad31da30ae64 82 time(&rtcTime);
acracan 0:ad31da30ae64 83 rtcTimeInfo = localtime(&rtcTime);
acracan 0:ad31da30ae64 84
acracan 0:ad31da30ae64 85 // Calculate difference in days between the current Julian Day
acracan 0:ad31da30ae64 86 // and JD 2451545.0, which is noon 1 January 2000 Universal Time
acracan 0:ad31da30ae64 87
acracan 0:ad31da30ae64 88 // Calculate time of the day in UT decimal hours
acracan 0:ad31da30ae64 89 dDecimalHours = (rtcTimeInfo->tm_hour - iTzOffset) + (rtcTimeInfo->tm_min
acracan 0:ad31da30ae64 90 + rtcTimeInfo->tm_sec / 60.0 ) / 60.0;
acracan 0:ad31da30ae64 91
acracan 0:ad31da30ae64 92 // Calculate current Julian Day
acracan 0:ad31da30ae64 93 long int iYfrom2000=rtcTimeInfo->tm_year-100;
acracan 0:ad31da30ae64 94 long int iA=(14-(rtcTimeInfo->tm_mon+1))/12;
acracan 0:ad31da30ae64 95 long int iM=(rtcTimeInfo->tm_mon+1)+12*iA-3;
acracan 0:ad31da30ae64 96
acracan 0:ad31da30ae64 97 long int liAux3=(153*iM+2)/5;
acracan 0:ad31da30ae64 98 long int liAux4=365*(iYfrom2000-iA);
acracan 0:ad31da30ae64 99 long int liAux5=(iYfrom2000-iA)/4;
acracan 0:ad31da30ae64 100
acracan 0:ad31da30ae64 101
acracan 0:ad31da30ae64 102 dElapsedJulianDays=(double)(rtcTimeInfo->tm_mday+liAux3+liAux4+liAux5+59)+
acracan 0:ad31da30ae64 103 -0.5+dDecimalHours/24.0;
acracan 0:ad31da30ae64 104
acracan 0:ad31da30ae64 105
acracan 0:ad31da30ae64 106 // Calculate ecliptic coordinates (ecliptic longitude and obliquity of the
acracan 0:ad31da30ae64 107 // ecliptic in radians but without limiting the angle to be less than 2*Pi
acracan 0:ad31da30ae64 108 // (i.e., the result may be greater than 2*Pi)
acracan 0:ad31da30ae64 109 dOmega=2.1429-0.0010394594*dElapsedJulianDays;
acracan 0:ad31da30ae64 110 dMeanLongitude = 4.8950630+ 0.017202791698*dElapsedJulianDays; // Radians
acracan 0:ad31da30ae64 111 dMeanAnomaly = 6.2400600+ 0.0172019699*dElapsedJulianDays;
acracan 0:ad31da30ae64 112 dEclipticLongitude = dMeanLongitude + 0.03341607*sin( dMeanAnomaly )
acracan 0:ad31da30ae64 113 + 0.00034894*sin( 2*dMeanAnomaly )-0.0001134
acracan 0:ad31da30ae64 114 -0.0000203*sin(dOmega);
acracan 0:ad31da30ae64 115 dEclipticObliquity = 0.4090928 - 6.2140e-9*dElapsedJulianDays
acracan 0:ad31da30ae64 116 +0.0000396*cos(dOmega);
acracan 0:ad31da30ae64 117
acracan 0:ad31da30ae64 118 // Calculate celestial coordinates ( right ascension and declination ) in radians
acracan 0:ad31da30ae64 119 // but without limiting the angle to be less than 2*Pi (i.e., the result may be
acracan 0:ad31da30ae64 120 // greater than 2*Pi)
acracan 0:ad31da30ae64 121 dSin_EclipticLongitude= sin( dEclipticLongitude );
acracan 0:ad31da30ae64 122 double dY1 = cos( dEclipticObliquity ) * dSin_EclipticLongitude;
acracan 0:ad31da30ae64 123 double dX1 = cos( dEclipticLongitude );
acracan 0:ad31da30ae64 124 dRightAscension = atan2( dY1,dX1 );
acracan 0:ad31da30ae64 125 if( dRightAscension < 0.0 ) dRightAscension = dRightAscension + twopi;
acracan 0:ad31da30ae64 126 dDeclination = asin( sin( dEclipticObliquity )*dSin_EclipticLongitude );
acracan 0:ad31da30ae64 127
acracan 0:ad31da30ae64 128 // Calculate local coordinates ( azimuth and zenith angle ) in degrees
acracan 0:ad31da30ae64 129 dGreenwichMeanSiderealTime = 6.6974243242 +
acracan 0:ad31da30ae64 130 0.0657098283*dElapsedJulianDays
acracan 0:ad31da30ae64 131 + dDecimalHours;
acracan 0:ad31da30ae64 132
acracan 0:ad31da30ae64 133 dLocalMeanSiderealTime = (dGreenwichMeanSiderealTime*15
acracan 0:ad31da30ae64 134 + dLongitude)*rad;
acracan 0:ad31da30ae64 135 dHourAngle = dLocalMeanSiderealTime - dRightAscension;
acracan 0:ad31da30ae64 136
acracan 0:ad31da30ae64 137 dCos_HourAngle= cos( dHourAngle );
acracan 0:ad31da30ae64 138 dZenithAngle = (acos( dCos_Latitude*dCos_HourAngle
acracan 0:ad31da30ae64 139 *cos(dDeclination) + sin( dDeclination )*dSin_Latitude));
acracan 0:ad31da30ae64 140 double dY = -sin( dHourAngle );
acracan 0:ad31da30ae64 141 double dX = tan( dDeclination )*dCos_Latitude - dSin_Latitude*dCos_HourAngle;
acracan 0:ad31da30ae64 142 dAzimuth=atan2( dY, dX );
acracan 0:ad31da30ae64 143 if ( dAzimuth < 0.0 )
acracan 0:ad31da30ae64 144 dAzimuth = dAzimuth + twopi;
acracan 0:ad31da30ae64 145 dAzimuth = dAzimuth/rad;
acracan 0:ad31da30ae64 146 // Parallax Correction
acracan 0:ad31da30ae64 147 dParallax=(dEarthMeanRadius/dAstronomicalUnit)
acracan 0:ad31da30ae64 148 *sin(dZenithAngle);
acracan 0:ad31da30ae64 149 dZenithAngle=(dZenithAngle
acracan 0:ad31da30ae64 150 + dParallax)/rad;
acracan 0:ad31da30ae64 151 dElevation=90-dZenithAngle;
acracan 0:ad31da30ae64 152 }
acracan 0:ad31da30ae64 153
acracan 0:ad31da30ae64 154 void Helios::setLocalLatitude(double latitude)
acracan 0:ad31da30ae64 155 {
acracan 0:ad31da30ae64 156 double dLatitudeInRadians;
acracan 0:ad31da30ae64 157
acracan 0:ad31da30ae64 158 dLatitude = latitude;
acracan 0:ad31da30ae64 159 dLatitudeInRadians = dLatitude*rad;
acracan 0:ad31da30ae64 160 dCos_Latitude = cos( dLatitudeInRadians );
acracan 0:ad31da30ae64 161 dSin_Latitude = sin( dLatitudeInRadians );
acracan 0:ad31da30ae64 162 }
acracan 0:ad31da30ae64 163
acracan 0:ad31da30ae64 164 void Helios::setLocalLongitude(double longitude)
acracan 0:ad31da30ae64 165 {
acracan 0:ad31da30ae64 166 dLongitude = longitude;
acracan 0:ad31da30ae64 167 }
acracan 0:ad31da30ae64 168
acracan 0:ad31da30ae64 169 void Helios::setLocalTimeZoneOffset(int tzOffset)
acracan 0:ad31da30ae64 170 {
acracan 0:ad31da30ae64 171 iTzOffset = tzOffset;
acracan 0:ad31da30ae64 172 }
acracan 0:ad31da30ae64 173
acracan 0:ad31da30ae64 174 double Helios::azimuth()
acracan 0:ad31da30ae64 175 {
acracan 0:ad31da30ae64 176 return dAzimuth;
acracan 0:ad31da30ae64 177 }
acracan 0:ad31da30ae64 178
acracan 0:ad31da30ae64 179 double Helios::elevation()
acracan 0:ad31da30ae64 180 {
acracan 0:ad31da30ae64 181 return dElevation;
acracan 0:ad31da30ae64 182 }