-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSunRiseSet.cpp
331 lines (155 loc) · 6.81 KB
/
SunRiseSet.cpp
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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
#define PI 3.1415926
#include<math.h>
#include<iostream>
using namespace std;
int days_of_month_1[]={31,28,31,30,31,30,31,31,30,31,30,31};
int days_of_month_2[]={31,29,31,30,31,30,31,31,30,31,30,31};
long double h=-0.833;
//定义全局变量
void input_date(int c[]){
int i;
cout<<"Enter the date (form: 2009 03 10):"<<endl;
for(i=0;i<3;i++){
cin>>c[i];
}
}
//输入日期
void input_glat(int c[]){
int i;
cout<<"Enter the degree of latitude(range: 0°- 60°,form: 40 40 40 (means 40°40′40″)):"<<endl;
for(i=0;i<3;i++){
cin>>c[i];
}
}
//输入纬度
void input_glong(int c[]){
int i;
cout<<"Enter the degree of longitude(west is negativ,form: 40 40 40 (means 40°40′40″)):"<<endl;
for(i=0;i<3;i++){
cin>>c[i];
}
}
//输入经度
int leap_year(int year){
if(((year%400==0) || (year%100!=0) && (year%4==0))) return 1;
else return 0;
}
//判断是否为闰年:若为闰年，返回1；若非闰年，返回0
int days(int year, int month, int date){
int i,a=0;
for(i=2000;i<year;i++){
if(leap_year(i)) a=a+366;
else a=a+365;
}
if(leap_year(year)){
for(i=0;i<month-1;i++){
a=a+days_of_month_2[i];
}
}
else {
for(i=0;i<month-1;i++){
a=a+days_of_month_1[i];
}
}
a=a+date;
return a;
}
//求从格林威治时间公元2000年1月1日到计算日天数days
long double t_century(int days, long double UTo){
return ((long double)days+UTo/360)/36525;
}
//求格林威治时间公元2000年1月1日到计算日的世纪数t
long double L_sun(long double t_century){
return (280.460+36000.770*t_century);
}
//求太阳的平黄径
long double G_sun(long double t_century){
return (357.528+35999.050*t_century);
}
//求太阳的平近点角
long double ecliptic_longitude(long double L_sun,long double G_sun){
return (L_sun+1.915*sin(G_sun*PI/180)+0.02*sin(2*G_sun*PI/180));
}
//求黄道经度
long double earth_tilt(long double t_century){
return (23.4393-0.0130*t_century);
}
//求地球倾角
long double sun_deviation(long double earth_tilt, long double ecliptic_longitude){
return (180/PI*asin(sin(PI/180*earth_tilt)*sin(PI/180*ecliptic_longitude)));
}
//求太阳偏差
long double GHA(long double UTo, long double G_sun, long double ecliptic_longitude){
return (UTo-180-1.915*sin(G_sun*PI/180)-0.02*sin(2*G_sun*PI/180)+2.466*sin(2*ecliptic_longitude*PI/180)-0.053*sin(4*ecliptic_longitude*PI/180));
}
//求格林威治时间的太阳时间角GHA
long double e(long double h, long double glat, long double sun_deviation){
return 180/PI*acos((sin(h*PI/180)-sin(glat*PI/180)*sin(sun_deviation*PI/180))/(cos(glat*PI/180)*cos(sun_deviation*PI/180)));
}
//求修正值e
long double UT_rise(long double UTo, long double GHA, long double glong, long double e){
return (UTo-(GHA+glong+e));
}
//求日出时间
long double UT_set(long double UTo, long double GHA, long double glong, long double e){
return (UTo-(GHA+glong-e));
}
//求日落时间
long double result_rise(long double UT, long double UTo, long double glong, long double glat, int year, int month, int date){
long double d;
if(UT>=UTo) d=UT-UTo;
else d=UTo-UT;
if(d>=0.1) {
UTo=UT;
UT=UT_rise(UTo,GHA(UTo,G_sun(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo)))),glong,e(h,glat,sun_deviation(earth_tilt(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo))))));
result_rise(UT,UTo,glong,glat,year,month,date);
}
return UT;
}
//判断并返回结果（日出）
long double result_set(long double UT, long double UTo, long double glong, long double glat, int year, int month, int date){
long double d;
if(UT>=UTo) d=UT-UTo;
else d=UTo-UT;
if(d>=0.1){
UTo=UT;
UT=UT_set(UTo,GHA(UTo,G_sun(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo)))),glong,e(h,glat,sun_deviation(earth_tilt(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo))))));
result_set(UT,UTo,glong,glat,year,month,date);
}
return UT;
}
//判断并返回结果（日落）
int Zone(long double glong){
if(glong>=0) return (int)((int)(glong/15.0)+1);
else return (int)((int)(glong/15.0)-1);
}
//求时区
void output(long double rise, long double set, long double glong){
if((int)(60*(rise/15+Zone(glong)-(int)(rise/15+Zone(glong))))<10)
cout<<"The time at which the sun rises is "<<(int)(rise/15+Zone(glong))<<":0"<<(int)(60*(rise/15+Zone(glong)-(int)(rise/15+Zone(glong))))<<" .\n";
else cout<<"The time at which the sun rises is "<<(int)(rise/15+Zone(glong))<<":"<<(int)(60*(rise/15+Zone(glong)-(int)(rise/15+Zone(glong))))<<" .\n";
if((int)(60*(set/15+Zone(glong)-(int)(set/15+Zone(glong))))<10)
cout<<"The time at which the sun sets is "<<(int)(set/15+Zone(glong))<<": "<<(int)(60*(set/15+Zone(glong)-(int)(set/15+Zone(glong))))<<" .\n";
else cout<<"The time at which the sun sets is "<<(int)(set/15+Zone(glong))<<":"<<(int)(60*(set/15+Zone(glong)-(int)(set/15+Zone(glong))))<<" .\n";
}
//打印结果
int main(){
long double UTo=180.0;
int year,month,date;
long double glat,glong;
int c[3];
input_date(c);
year=c[0];
month=c[1];
date=c[2];
input_glat(c);
glat=c[0]+c[1]/60+c[2]/3600;
input_glong(c);
glong=c[0]+c[1]/60+c[2]/3600;
long double rise,set;
rise=result_rise(UT_rise(UTo,GHA(UTo,G_sun(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo)))),glong,e(h,glat,sun_deviation(earth_tilt(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo)))))),UTo,glong,glat,year,month,date);
set=result_set(UT_set(UTo,GHA(UTo,G_sun(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo)))),glong,e(h,glat,sun_deviation(earth_tilt(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo)))))),UTo,glong,glat,year,month,date);
output(rise,set,glong);
system("pause");
return 0;
}