| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // Compute times of sunise and sunset at a specified latitude and longitude. | ||
| 2 | // | ||
| 3 | // This software minimizes computational work by performing the full calculation | ||
| 4 | // of the solar position three times, at the beginning, middle, and end of the | ||
| 5 | // period of interest. Three point interpolation is used to predict the position | ||
| 6 | // for each hour, and the arithmetic mean is used to predict the half-hour positions. | ||
| 7 | // | ||
| 8 | // The full computational burden is negligible on modern computers, but the | ||
| 9 | // algorithm is effective and still useful for small embedded systems. | ||
| 10 | // | ||
| 11 | // This software was originally adapted to javascript by Stephen R. Schmitt | ||
| 12 | // from a BASIC program from the 'Astronomical Computing' column of Sky & Telescope, | ||
| 13 | // August 1994, page 84, written by Roger W. Sinnott. | ||
| 14 | // | ||
| 15 | // Subsequently adapted from Stephen R. Schmitt's javascript to c++ for the Arduino | ||
| 16 | // by Cyrus Rahman, this work is subject to Stephen Schmitt's copyright: | ||
| 17 | // | ||
| 18 | // Copyright 2007 Stephen R. Schmitt | ||
| 19 | // Subsequent work Copyright 2020 Cyrus Rahman | ||
| 20 | // You may use or modify this source code in any way you find useful, provided | ||
| 21 | // that you agree that the author(s) have no warranty, obligations or liability. You | ||
| 22 | // must determine the suitability of this source code for your use. | ||
| 23 | // | ||
| 24 | // Redistributions of this source code must retain this copyright notice. | ||
| 25 | |||
| 26 | #include <math.h> | ||
| 27 | #include "SunRise.h" | ||
| 28 | |||
| 29 | #define K1 15*(M_PI/180)*1.0027379 | ||
| 30 | |||
| 31 | struct skyCoordinates { | ||
| 32 | double RA; // Right ascension | ||
| 33 | double declination; // Declination | ||
| 34 | }; | ||
| 35 | |||
| 36 | // Determine the nearest sun rise or set event previous, and the nearest | ||
| 37 | // sun rise or set event subsequent, to the specified time in seconds since the | ||
| 38 | // Unix epoch (January 1, 1970) and at the specified latitude and longitude in | ||
| 39 | // degrees. | ||
| 40 | // | ||
| 41 | // We look for events from SR_WINDOW/2 hours in the past to SR_WINDOW/2 hours | ||
| 42 | // in the future. | ||
| 43 | void | ||
| 44 | 325 | SunRise::calculate(double latitude, double longitude, time_t t, TwilightType type) { | |
| 45 | skyCoordinates sunPosition[3]; | ||
| 46 | double offsetDays; | ||
| 47 | |||
| 48 |
1/2✓ Branch 0 taken 325 times.
✗ Branch 1 not taken.
|
325 | initClass(); |
| 49 | 325 | queryTime = t; | |
| 50 |
1/2✓ Branch 0 taken 325 times.
✗ Branch 1 not taken.
|
325 | offsetDays = julianDate(t) - 2451545L; // Days since Jan 1, 2000, 1200UTC. |
| 51 | // Begin testing (SR_WINDOW / 2) hours before requested time. | ||
| 52 | 325 | offsetDays -= (double)SR_WINDOW / (2 * 24) ; | |
| 53 | |||
| 54 | // Calculate coordinates at start, middle, and end of search period. | ||
| 55 |
2/2✓ Branch 0 taken 975 times.
✓ Branch 1 taken 325 times.
|
1300 | for (int i = 0; i < 3; i ++) { |
| 56 |
1/2✓ Branch 0 taken 975 times.
✗ Branch 1 not taken.
|
975 | sunPosition[i] = sun(offsetDays + i * (double)SR_WINDOW / (2 * 24)); |
| 57 | } | ||
| 58 | |||
| 59 | // If the RA wraps around during this period, unwrap it to keep the | ||
| 60 | // sequence smooth for interpolation. | ||
| 61 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 325 times.
|
325 | if (sunPosition[1].RA <= sunPosition[0].RA) |
| 62 | ✗ | sunPosition[1].RA += 2 * M_PI; | |
| 63 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 325 times.
|
325 | if (sunPosition[2].RA <= sunPosition[1].RA) |
| 64 | ✗ | sunPosition[2].RA += 2 * M_PI; | |
| 65 | |||
| 66 | // Initialize interpolation array. | ||
| 67 | skyCoordinates spWindow[3]; | ||
| 68 | 325 | spWindow[0].RA = sunPosition[0].RA; | |
| 69 | 325 | spWindow[0].declination = sunPosition[0].declination; | |
| 70 | |||
| 71 |
2/2✓ Branch 0 taken 7800 times.
✓ Branch 1 taken 325 times.
|
8125 | for (int k = 0; k < SR_WINDOW; k++) { // Check each interval of search period |
| 72 | 7800 | float ph = (float)(k + 1)/SR_WINDOW; | |
| 73 | |||
| 74 |
1/2✓ Branch 0 taken 7800 times.
✗ Branch 1 not taken.
|
7800 | spWindow[2].RA = interpolate(sunPosition[0].RA, |
| 75 | sunPosition[1].RA, | ||
| 76 | sunPosition[2].RA, ph); | ||
| 77 |
1/2✓ Branch 0 taken 7800 times.
✗ Branch 1 not taken.
|
7800 | spWindow[2].declination = interpolate(sunPosition[0].declination, |
| 78 | sunPosition[1].declination, | ||
| 79 | sunPosition[2].declination, ph); | ||
| 80 | |||
| 81 | // Look for sunrise/set events during this interval. | ||
| 82 |
1/2✓ Branch 0 taken 7800 times.
✗ Branch 1 not taken.
|
7800 | testSunRiseSet(k, offsetDays, latitude, longitude, spWindow, type); |
| 83 | |||
| 84 | 7800 | spWindow[0] = spWindow[2]; // Advance to next interval. | |
| 85 | } | ||
| 86 | 325 | } | |
| 87 | |||
| 88 | // Look for sun rise or set events during an hour. | ||
| 89 | void | ||
| 90 | 7800 | SunRise::testSunRiseSet(int k, double offsetDays, double latitude, double longitude, | |
| 91 | skyCoordinates *sp, TwilightType type) { | ||
| 92 | double ha[3], VHz[3]; | ||
| 93 | double lSideTime; | ||
| 94 | |||
| 95 | // Get (local_sidereal_time - SR_WINDOW / 2) hours in radians. | ||
| 96 |
1/2✓ Branch 0 taken 7800 times.
✗ Branch 1 not taken.
|
7800 | lSideTime = localSiderealTime(offsetDays, longitude) * 2* M_PI / 360; |
| 97 | |||
| 98 | // Calculate Hour Angle. | ||
| 99 | 7800 | ha[0] = lSideTime - sp[0].RA + k*K1; | |
| 100 | 7800 | ha[2] = lSideTime - sp[2].RA + k*K1 + K1; | |
| 101 | |||
| 102 | // Hour Angle and declination at half hour. | ||
| 103 | 7800 | ha[1] = (ha[2] + ha[0])/2; | |
| 104 | 7800 | sp[1].declination = (sp[2].declination + sp[0].declination)/2; | |
| 105 | |||
| 106 | 7800 | double s = sin(M_PI / 180 * latitude); | |
| 107 | 7800 | double c = cos(M_PI / 180 * latitude); | |
| 108 | |||
| 109 | // refraction + sun semidiameter at horizon | ||
| 110 | 7800 | double z = 0; | |
| 111 | |||
| 112 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7800 times.
|
7800 | if (type == TT_ASTRONOMICAL) |
| 113 | ✗ | z = cos(M_PI / 180 * 108.0); | |
| 114 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7800 times.
|
7800 | else if (type == TT_NAUTICAL) |
| 115 | ✗ | z = cos(M_PI / 180 * 102.0); | |
| 116 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7800 times.
|
7800 | else if (type == TT_CIVIL) |
| 117 | ✗ | z = cos(M_PI / 180 * 96.0); | |
| 118 | else | ||
| 119 | 7800 | z = cos(M_PI / 180 * 90.833); | |
| 120 | |||
| 121 | 7800 | VHz[0] = s * sin(sp[0].declination) + c * cos(sp[0].declination) * cos(ha[0]) - z; | |
| 122 | 7800 | VHz[2] = s * sin(sp[2].declination) + c * cos(sp[2].declination) * cos(ha[2]) - z; | |
| 123 | |||
| 124 |
2/2✓ Branch 0 taken 7220 times.
✓ Branch 1 taken 580 times.
|
7800 | if (signbit(VHz[0]) == signbit(VHz[2])) |
| 125 | 7220 | goto noevent; // No event this hour. | |
| 126 | |||
| 127 | 580 | VHz[1] = s * sin(sp[1].declination) + c * cos(sp[1].declination) * cos(ha[1]) - z; | |
| 128 | |||
| 129 | double a, b, d, e, time; | ||
| 130 | 580 | a = 2 * VHz[2] - 4 * VHz[1] + 2 * VHz[0]; | |
| 131 | 580 | b = 4 * VHz[1] - 3 * VHz[0] - VHz[2]; | |
| 132 | 580 | d = b * b - 4 * a * VHz[0]; | |
| 133 | |||
| 134 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 580 times.
|
580 | if (d < 0) |
| 135 | ✗ | goto noevent; // No event this hour. | |
| 136 | |||
| 137 | 580 | d = sqrt(d); | |
| 138 | 580 | e = (-b + d) / (2 * a); | |
| 139 |
4/4✓ Branch 0 taken 576 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 286 times.
✓ Branch 3 taken 290 times.
|
580 | if ((e < 0) || (e > 1)) |
| 140 | 290 | e = (-b - d) / (2 * a); | |
| 141 | 580 | time = k + e + 1 / 120; // Time since k=0 of event (in hours). | |
| 142 | |||
| 143 | // The time we started searching + the time from the start of the search to the | ||
| 144 | // event is the time of the event. Add (time since k=0) - window/2 hours. | ||
| 145 | time_t eventTime; | ||
| 146 | 580 | eventTime = queryTime + (time - SR_WINDOW / 2) *60 *60; | |
| 147 | |||
| 148 | double hz, nz, dz, az; | ||
| 149 | 580 | hz = ha[0] + e * (ha[2] - ha[0]); // Azimuth of the sun at the event. | |
| 150 | 580 | nz = -cos(sp[1].declination) * sin(hz); | |
| 151 | 580 | dz = c * sin(sp[1].declination) - s * cos(sp[1].declination) * cos(hz); | |
| 152 | 580 | az = atan2(nz, dz) / (M_PI / 180); | |
| 153 |
2/2✓ Branch 0 taken 290 times.
✓ Branch 1 taken 290 times.
|
580 | if (az < 0) |
| 154 | 290 | az += 360; | |
| 155 | |||
| 156 | // If there is no previously recorded event of this type, save this event. | ||
| 157 | // | ||
| 158 | // If this event is previous to queryTime, and is the nearest event to queryTime | ||
| 159 | // of events of its type previous to queryType, save this event, replacing the | ||
| 160 | // previously recorded event of its type. Events subsequent to queryTime are | ||
| 161 | // treated similarly, although since events are tested in chronological order | ||
| 162 | // no replacements will occur as successive events will be further from | ||
| 163 | // queryTime. | ||
| 164 | // | ||
| 165 | // If this event is subsequent to queryTime and there is an event of its type | ||
| 166 | // previous to queryTime, then there is an event of the other type between the | ||
| 167 | // two events of this event's type. If the event of the other type is | ||
| 168 | // previous to queryTime, then it is the nearest event to queryTime that is | ||
| 169 | // previous to queryTime. In this case save the current event, replacing | ||
| 170 | // the previously recorded event of its type. Otherwise discard the current | ||
| 171 | // event. | ||
| 172 | // | ||
| 173 |
3/4✓ Branch 0 taken 290 times.
✓ Branch 1 taken 290 times.
✓ Branch 2 taken 290 times.
✗ Branch 3 not taken.
|
580 | if ((VHz[0] < 0) && (VHz[2] > 0)) { |
| 174 | ✗ | if (!hasRise || | |
| 175 | ✗ | (signbit(riseTime - queryTime) == signbit(eventTime - queryTime) && | |
| 176 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 290 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 290 times.
✗ Branch 5 not taken.
|
290 | fabs(riseTime - queryTime) > fabs(eventTime - queryTime)) || |
| 177 | ✗ | (signbit(riseTime - queryTime) != signbit(eventTime - queryTime) && | |
| 178 | ✗ | (hasSet && | |
| 179 | ✗ | signbit(riseTime - queryTime) == signbit(setTime - queryTime)))) { | |
| 180 | 290 | riseTime = eventTime; | |
| 181 | 290 | riseAz = az; | |
| 182 | 290 | hasRise = true; | |
| 183 | } | ||
| 184 | } | ||
| 185 |
3/4✓ Branch 0 taken 290 times.
✓ Branch 1 taken 290 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 290 times.
|
580 | if ((VHz[0] > 0) && (VHz[2] < 0)) { |
| 186 | ✗ | if (!hasSet || | |
| 187 | ✗ | (signbit(setTime - queryTime) == signbit(eventTime - queryTime) && | |
| 188 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 290 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 290 times.
✗ Branch 5 not taken.
|
290 | fabs(setTime - queryTime) > fabs(eventTime - queryTime)) || |
| 189 | ✗ | (signbit(setTime - queryTime) != signbit(eventTime - queryTime) && | |
| 190 | ✗ | (hasRise && | |
| 191 | ✗ | signbit(setTime - queryTime) == signbit(riseTime - queryTime)))) { | |
| 192 | 290 | setTime = eventTime; | |
| 193 | 290 | setAz = az; | |
| 194 | 290 | hasSet = true; | |
| 195 | } | ||
| 196 | } | ||
| 197 | |||
| 198 | 290 | noevent: | |
| 199 | // There are obscure cases in the polar regions that require extra logic. | ||
| 200 |
3/4✓ Branch 0 taken 2484 times.
✓ Branch 1 taken 5316 times.
✓ Branch 2 taken 2484 times.
✗ Branch 3 not taken.
|
7800 | if (!hasRise && !hasSet) |
| 201 | 2484 | isVisible = !signbit(VHz[2]); | |
| 202 |
3/4✓ Branch 0 taken 5316 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4568 times.
✓ Branch 3 taken 748 times.
|
5316 | else if (hasRise && !hasSet) |
| 203 | 4568 | isVisible = (queryTime > riseTime); | |
| 204 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 748 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
748 | else if (!hasRise && hasSet) |
| 205 | ✗ | isVisible = (queryTime < setTime); | |
| 206 | else | ||
| 207 |
3/6✓ Branch 0 taken 748 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 748 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 748 times.
|
748 | isVisible = ((riseTime < setTime && riseTime < queryTime && setTime > queryTime) || |
| 208 | ✗ | (riseTime > setTime && (riseTime < queryTime || setTime > queryTime))); | |
| 209 | |||
| 210 | 15600 | return; | |
| 211 | } | ||
| 212 | |||
| 213 | // Sun position using fundamental arguments | ||
| 214 | // (Van Flandern & Pulkkinen, 1979) | ||
| 215 | skyCoordinates | ||
| 216 | 975 | SunRise::sun(double dayOffset) { | |
| 217 | 975 | double centuryOffset = dayOffset / 36525 + 1; // Centuries from 1900.0 | |
| 218 | |||
| 219 | 975 | double l = 0.779072 + 0.00273790931 * dayOffset; | |
| 220 | 975 | double g = 0.993126 + 0.00273777850 * dayOffset; | |
| 221 | |||
| 222 | 975 | l = 2 * M_PI * (l - floor(l)); | |
| 223 | 975 | g = 2 * M_PI * (g - floor(g)); | |
| 224 | |||
| 225 | double v, u, w; | ||
| 226 | 975 | v = 0.39785 * sin(l) | |
| 227 | 975 | - 0.01000 * sin(l - g) | |
| 228 | 975 | + 0.00333 * sin(l + g) | |
| 229 | 975 | - 0.00021 * centuryOffset * sin(l); | |
| 230 | |||
| 231 | 975 | u = 1 | |
| 232 | 975 | - 0.03349 * cos(g) | |
| 233 | 975 | - 0.00014 * cos(2*l) | |
| 234 | 975 | + 0.00008 * cos(l); | |
| 235 | |||
| 236 | 975 | w = -0.00010 | |
| 237 | 975 | - 0.04129 * sin(2*l) | |
| 238 | 975 | + 0.03211 * sin(g) | |
| 239 | 975 | + 0.00104 * sin(2*l - g) | |
| 240 | 975 | - 0.00035 * sin(2*l + g) | |
| 241 | 975 | - 0.00008 * centuryOffset * sin(g); | |
| 242 | |||
| 243 | double s; | ||
| 244 | skyCoordinates sc; | ||
| 245 | 975 | s = w / sqrt(u - v*v); // Right ascension | |
| 246 | 975 | sc.RA = l + atan(s / sqrt(1 - s*s)); | |
| 247 | |||
| 248 | 975 | s = v / sqrt(u); // Declination | |
| 249 | 975 | sc.declination = atan(s / sqrt(1 - s*s)); | |
| 250 | 975 | return(sc); | |
| 251 | } | ||
| 252 | |||
| 253 | // 3-point interpolation | ||
| 254 | double | ||
| 255 | 15600 | SunRise::interpolate(double f0, double f1, double f2, double p) { | |
| 256 | 15600 | double a = f1 - f0; | |
| 257 | 15600 | double b = f2 - f1 - a; | |
| 258 | 15600 | return(f0 + p * (2*a + b * (2*p - 1))); | |
| 259 | } | ||
| 260 | |||
| 261 | // Determine Julian date from Unix time. | ||
| 262 | // Provides marginally accurate results with Arduino 4-byte double. | ||
| 263 | double | ||
| 264 | 325 | SunRise::julianDate(time_t t) { | |
| 265 | 325 | return (t / 86400.0L + 2440587.5); | |
| 266 | } | ||
| 267 | |||
| 268 | #if __ISO_C_VISIBLE < 1999 | ||
| 269 | // Arduino compiler is missing this function as of 6/2020. | ||
| 270 | // | ||
| 271 | // The Arduino ATmega platforms (including the Uno) are also missing rint(). | ||
| 272 | // This can be worked around by inserting "#define rint(x) (double)lrint(x)" here, | ||
| 273 | // but since these platforms use only four bytes for double precision - which is | ||
| 274 | // insufficient for the correct performance of the required calculations - you | ||
| 275 | // should instead upgrade to an Arduino Due or better. | ||
| 276 | // | ||
| 277 | #define remainder(x, y) ((double)((double)x - (double)y * rint((double)x / (double)y))) | ||
| 278 | |||
| 279 | // double remainder(double x, double y) { | ||
| 280 | // return(x - (y * rint(x / y))); | ||
| 281 | // } | ||
| 282 | #endif | ||
| 283 | |||
| 284 | // Local Sidereal Time | ||
| 285 | // Provides local sidereal time in degrees, requires longitude in degrees | ||
| 286 | // and time in fractional Julian days since Jan 1, 2000, 1200UTC (e.g. the | ||
| 287 | // Julian date - 2451545). | ||
| 288 | // cf. USNO Astronomical Almanac and | ||
| 289 | // https://astronomy.stackexchange.com/questions/24859/local-sidereal-time | ||
| 290 | double | ||
| 291 | 7800 | SunRise::localSiderealTime(double offsetDays, double longitude) { | |
| 292 | 7800 | double lSideTime = (15.0L * (6.697374558L + 0.06570982441908L * offsetDays + | |
| 293 | 7800 | remainder(offsetDays, 1) * 24 + 12 + | |
| 294 | 7800 | 0.000026 * (offsetDays / 36525) * (offsetDays / 36525)) | |
| 295 | 7800 | + longitude) / 360; | |
| 296 | 7800 | lSideTime -= floor(lSideTime); | |
| 297 | 7800 | lSideTime *= 360; // Convert to degrees. | |
| 298 | 7800 | return(lSideTime); | |
| 299 | } | ||
| 300 | |||
| 301 | // Class initialization. | ||
| 302 | void | ||
| 303 | 325 | SunRise::initClass() { | |
| 304 | 325 | queryTime = 0; | |
| 305 | 325 | riseTime = 0; | |
| 306 | 325 | setTime = 0; | |
| 307 | 325 | riseAz = 0; | |
| 308 | 325 | setAz = 0; | |
| 309 | 325 | hasRise = false; | |
| 310 | 325 | hasSet = false; | |
| 311 | 325 | isVisible = false; | |
| 312 | 325 | } | |
| 313 |