先贴在这了，后面应该用得着

 

1. #define PI 3.1415926
3. #include<math.h>
5. #include<iostream>
7. using namespace std;
10. int days_of_month_1[]={
    31,28,31,30,31,30,31,31,30,31,30,31};
12. int days_of_month_2[]={
    31,29,31,30,31,30,31,31,30,31,30,31};
14. long double h=-0.833;
16. //定义全局变量
19. void input_date(int c[]){
21. int i;
23. cout<<"Enter the date (form: 2009 03 10):"<<endl;
25. for(i=0;i<3;i++){
27. cin>>c[i];
29. }
31. }
33. //输入日期
36. void input_glat(int c[]){
38. int i;
40. cout<<"Enter the degree of latitude(range: 0°- 60°,form: 40 40 40 (means 40°40′40″)):"<<endl;
42. for(i=0;i<3;i++){
44. cin>>c[i];
46. }
48. }
50. //输入纬度
53. void input_glong(int c[]){
55. int i;
57. cout<<"Enter the degree of longitude(west is negativ,form: 40 40 40 (means 40°40′40″)):"<<endl;
59. for(i=0;i<3;i++){
61. cin>>c[i];
63. }
65. }
67. //输入经度
70. int leap_year(int year){
72. if(((year%400==0) || (year%100!=0) && (year%4==0))) return 1;
74. else return 0;
76. }
78. //判断是否为闰年:若为闰年，返回1；若非闰年，返回0
81. int days(int year, int month, int date){
83. int i,a=0;
85. for(i=2000;i<year;i++){
87. if(leap_year(i)) a=a+366;
89. else a=a+365;
91. }
93. if(leap_year(year)){
95. for(i=0;i<month-1;i++){
97. a=a+days_of_month_2[i];
99. }
101. }
103. else {
105. for(i=0;i<month-1;i++){
107. a=a+days_of_month_1[i];
109. }
111. }
113. a=a+date;
115. return a;
117. }
119. //求从格林威治时间公元2000年1月1日到计算日天数days
122. long double t_century(int days, long double UTo){
124. return ((long double)days+UTo/360)/36525;
126. }
128. //求格林威治时间公元2000年1月1日到计算日的世纪数t
131. long double L_sun(long double t_century){
133. return (280.460+36000.770*t_century);
135. }
137. //求太阳的平黄径
140. long double G_sun(long double t_century){
142. return (357.528+35999.050*t_century);
144. }
146. //求太阳的平近点角
149. long double ecliptic_longitude(long double L_sun,long double G_sun){
151. return (L_sun+1.915*sin(G_sun*PI/180)+0.02*sin(2*G_sun*PI/180));
153. }
155. //求黄道经度
158. long double earth_tilt(long double t_century){
160. return (23.4393-0.0130*t_century);
162. }
164. //求地球倾角
167. long double sun_deviation(long double earth_tilt, long double ecliptic_longitude){
169. return (180/PI*asin(sin(PI/180*earth_tilt)*sin(PI/180*ecliptic_longitude)));
171. }
173. //求太阳偏差
176. long double GHA(long double UTo, long double G_sun, long double ecliptic_longitude){
178. 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));
180. }
182. //求格林威治时间的太阳时间角GHA
185. long double e(long double h, long double glat, long double sun_deviation){
187. 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)));
189. }
191. //求修正值e
194. long double UT_rise(long double UTo, long double GHA, long double glong, long double e){
196. return (UTo-(GHA+glong+e));
198. }
200. //求日出时间
203. long double UT_set(long double UTo, long double GHA, long double glong, long double e){
205. return (UTo-(GHA+glong-e));
207. }
209. //求日落时间
212. long double result_rise(long double UT, long double UTo, long double glong, long double glat, int year, int month, int date){
214. long double d;
216. if(UT>=UTo) d=UT-UTo;
218. else d=UTo-UT;
220. if(d>=0.1) {
222. UTo=UT;
224. 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))))));
226. result_rise(UT,UTo,glong,glat,year,month,date);
228. }
230. return UT;
232. }
234. //判断并返回结果（日出）
237. long double result_set(long double UT, long double UTo, long double glong, long double glat, int year, int month, int date){
239. long double d;
241. if(UT>=UTo) d=UT-UTo;
243. else d=UTo-UT;
245. if(d>=0.1){
247. UTo=UT;
249. 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))))));
251. result_set(UT,UTo,glong,glat,year,month,date);
253. }
255. return UT;
257. }
259. //判断并返回结果（日落）
262. int Zone(long double glong){
264. if(glong>=0) return (int)((int)(glong/15.0)+1);
266. else return (int)((int)(glong/15.0)-1);
268. }
270. //求时区
273. void output(long double rise, long double set, long double glong){
275. if((int)(60*(rise/15+Zone(glong)-(int)(rise/15+Zone(glong))))<10)
277. 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";
279. 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";
281. if((int)(60*(set/15+Zone(glong)-(int)(set/15+Zone(glong))))<10)
283. 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";
285. 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";
287. }
289. //打印结果
293. int main(){
295. long double UTo=180.0;
297. int year,month,date;
299. long double glat,glong;
301. int c[3];
303. input_date(c);
305. year=c[0];
307. month=c[1];
309. date=c[2];
311. input_glat(c);
313. glat=c[0]+c[1]/60+c[2]/3600;
315. input_glong(c);
317. glong=c[0]+c[1]/60+c[2]/3600;
319. long double rise,set;
321. 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);
323. 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);
325. output(rise,set,glong);
327. system("pause");
329. return 0;
331. }