blob: 44d9fa1ea98c84bb554e0e55ba7291e21ce583ed (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
|
#include <stdio.h>
#include <string.h>
#include <stdbool.h>
#include <math.h>
#include <time.h>
#include "../include/sun.h"
double deg_to_rad(double angle) {
return angle * (M_PI / 180);
}
double rad_to_deg(double angle) {
return angle * (180 / M_PI);
}
double sun_calculation(struct tm curr_time, double longitude, double latitude) {
// Calulate day of the year
int day_in_year = floor(275 * (curr_time.tm_mon + 1) / 9) -
((floor(((curr_time.tm_mon + 1) + 9) / 12)) *
(1 + floor(((curr_time.tm_year + 1900) - 4 * floor((curr_time.tm_year + 1900) / 4) + 2) / 3))) +
curr_time.tm_mday - 30;
// Longitude to hour, get approx time
double long_hour = longitude / 15;
double approx_time;
if (curr_time.tm_hour < 12)
approx_time = day_in_year + ((6 - long_hour) / 24);
else
approx_time = day_in_year + ((18 - long_hour) / 24);
// Sun's mean anomaly
double mean_anomaly = (0.9856 * approx_time) - 3.289;
// Sun's true longitude
double true_longitude = mean_anomaly + (1.916 * sin(deg_to_rad(mean_anomaly)) +
(0.02 * sin(deg_to_rad(2 * mean_anomaly)) + 282.634));
// Keep longitude between [0, 360)
if (true_longitude >= 360)
true_longitude -= 360;
else if (true_longitude < 0)
true_longitude += 360;
// Sun's right ascension
double right_ascension = rad_to_deg(atan(0.91764 * tan(deg_to_rad(true_longitude))));
// Keep Right Ascension between [0, 360)
if (right_ascension >= 360)
right_ascension -= 360;
else if (right_ascension < 0)
right_ascension += 360;
// Quadrant calibration
int longitude_quadrant = (floor(true_longitude / 90)) * 90;
int right_ascension_quadrant = (floor(right_ascension / 90)) * 90;
right_ascension += (longitude_quadrant - right_ascension_quadrant);
right_ascension /= 15; // convert to hours
// Calculate Sun's declination
double sin_declination = 0.39782 * sin(deg_to_rad(true_longitude));
double cos_declination = cos(asin(sin_declination));
// Calculate Sun's local hour angle
double sun_hour;
double cos_hour = (-0.01454 - (sin_declination * sin(deg_to_rad(latitude)))) /
(cos_declination * cos(deg_to_rad(latitude)));
if (curr_time.tm_hour < 12)
sun_hour = 360 - rad_to_deg(acos(cos_hour));
else
sun_hour = rad_to_deg(acos(cos_hour));
sun_hour /= 15; // convert to hours
// Local time
double local_mean_time = sun_hour + right_ascension - (0.06571 * approx_time) - 6.622;
// Adjust to UTC
double sun_UTC = local_mean_time - long_hour;
// Keep Sun UTC between [0, 24)
if (sun_UTC >= 24)
sun_UTC -= 24;
else if (sun_UTC < 0)
sun_UTC += 24;
double sun_time = sun_UTC + LOCAL_OFFSET;
if (sun_time < 0)
return (24 + sun_time);
else if (sun_time > 24)
return (sun_time - 24);
else
return sun_time;
}
|
