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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
|
// Encoding: UTF-8
#include "calendar.h"
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <iomanip>
// 常量
#define day_secs 86400
#define delta_time 1e-11
/* calendar.h规定之量 */
Julian julian;
parameter p;
double epsilon;
FILE *fp;
pair<double, double> getSunPos(double t, double flag = false) {
pair<double, double> sun = p.sun_longitude(t);
double alpha, delta; // 太阳赤经赤纬
double lambda = sun.first, beta = sun.second; // 太阳黄经黄纬
// 黄赤交角epsilon
epsilon = p.get_epsilon(t);
delta =
asin(sin(epsilon) * sin(lambda) * cos(beta) + cos(epsilon) * sin(beta));
alpha = atan2(
(cos(epsilon) * sin(lambda) * cos(beta) - sin(epsilon) * sin(beta)),
(cos(lambda) * cos(beta)));
// 时角
double H = p.getHourAngle(t, LONGITUDE, LATITUDE, alpha);
double A, h; //方位角和高度角
h = asin(sin(LATITUDE) * sin(delta) + cos(LATITUDE) * cos(delta) * cos(H));
A = atan2(-cos(delta) * sin(H),
cos(LATITUDE) * sin(delta) - sin(LATITUDE) * cos(delta) * cos(H));
if (flag) {
fprintf(
fp, "%s\t%s\t%s\t%s\t%s\t%s\n", radToDMS(0, 360, lambda).c_str(),
radToDMS(-90, 90, beta).c_str(), radToDMS(0, 360, alpha).c_str(),
radToDMS(-90, 90, delta).c_str(), radToDMS(0, 360, A).c_str(),
radToDMS(-90, 90, h).c_str());
}
return make_pair(A, h + M_PI / 216);
}
void binarySearch(string message, time_t begin, Date date) {
double t = julian.getJulianKiloYear(begin);
double t_end = julian.getJulianKiloYear(begin + day_secs / 2);
double middle;
time_t res;
pair<double, double> sunPos;
while (t_end - t >= delta_time) {
middle = (t + t_end) / 2;
sunPos = getSunPos(middle, false);
if (sunPos.second * getSunPos(t).second > 0) {
t = middle;
} else {
t_end = middle;
}
}
res = julian.kiloYearToTime(middle, date.tm_year + 1900);
Date *tm = localtime(&res);
printf("%s\t%4d-%02d-%02d %02d:%02d:%02d\tAzimuth\t%s\n", message.c_str(),
tm->tm_year + 1900, tm->tm_mon + 1, tm->tm_mday, tm->tm_hour,
tm->tm_min, tm->tm_sec, radToDMS(0, 360, sunPos.first).c_str());
fprintf(fp, "%s\t%4d-%02d-%02d %02d:%02d:%02d\tAzimuth\t%s\n",
message.c_str(), tm->tm_year + 1900, tm->tm_mon + 1, tm->tm_mday,
tm->tm_hour, tm->tm_min, tm->tm_sec,
radToDMS(0, 360, sunPos.first).c_str());
}
void getSunTime(Date date) {
time_t time = mktime(&date);
binarySearch("Sunrize", time, date);
binarySearch("Sunset", time + day_secs / 2, date);
}
void showSunPos(Date date) {
time_t time = mktime(&date);
double t;
fprintf(fp, "\n\nSun Trace Table\n\n");
fprintf(fp, "--------------------------------------------------------------"
"--------------------------------------------------------------"
"----------\n");
fprintf(fp, "Time |\tlongitude\t\tlatitude\t\tright "
"ascension\tdeclination\tAzimuth\t\tAltitude\n");
fprintf(fp, "--------------------------------------------------------------"
"--------------------------------------------------------------"
"----------\n");
for (int i = 0; i < 24; i++, time += 3600, date.tm_hour++) {
t = julian.getJulianKiloYear(time);
fprintf(fp, "%02d:%02d:%02d |\t", date.tm_hour, date.tm_min,
date.tm_sec);
getSunPos(t, true);
}
fprintf(fp, "--------------------------------------------------------------"
"--------------------------------------------------------------"
"----------\n");
}
int main() {
Date date;
cout << "Input a day you want(in YYYY-MM-DD): ";
scanf("%d-%d-%d", &date.tm_year, &date.tm_mon, &date.tm_mday);
date.tm_year -= 1900;
date.tm_mon -= 1;
date.tm_hour = date.tm_min = date.tm_sec = 0;
fp = fopen("sunPos.txt", "w");
if (fp == NULL) {
printf("File open failed\n");
exit(1);
}
fprintf(fp, "Sun trace File (Date: %4d-%02d-%02d)\n\n\n",
date.tm_year + 1900, date.tm_mon + 1, date.tm_mday);
getSunTime(date);
showSunPos(date);
fclose(fp);
return 0;
}
|
