|
| 1 | +#define PI 3.1415926 |
| 2 | + |
| 3 | + #include<math.h> |
| 4 | + |
| 5 | + #include<iostream> |
| 6 | + |
| 7 | + using namespace std; |
| 8 | + |
| 9 | + |
| 10 | + int days_of_month_1[]={31,28,31,30,31,30,31,31,30,31,30,31}; |
| 11 | + |
| 12 | + int days_of_month_2[]={31,29,31,30,31,30,31,31,30,31,30,31}; |
| 13 | + |
| 14 | + long double h=-0.833; |
| 15 | + |
| 16 | + //定义全局变量 |
| 17 | + |
| 18 | + |
| 19 | +void input_date(int c[]){ |
| 20 | + |
| 21 | + int i; |
| 22 | + |
| 23 | + cout<<"Enter the date (form: 2009 03 10):"<<endl; |
| 24 | + |
| 25 | + for(i=0;i<3;i++){ |
| 26 | + |
| 27 | + cin>>c[i]; |
| 28 | + |
| 29 | + } |
| 30 | + |
| 31 | + } |
| 32 | + |
| 33 | + //输入日期 |
| 34 | + |
| 35 | + |
| 36 | +void input_glat(int c[]){ |
| 37 | + |
| 38 | + int i; |
| 39 | + |
| 40 | + cout<<"Enter the degree of latitude(range: 0°- 60°,form: 40 40 40 (means 40°40′40″)):"<<endl; |
| 41 | + |
| 42 | + for(i=0;i<3;i++){ |
| 43 | + |
| 44 | + cin>>c[i]; |
| 45 | + |
| 46 | + } |
| 47 | + |
| 48 | + } |
| 49 | + |
| 50 | + //输入纬度 |
| 51 | + |
| 52 | + |
| 53 | +void input_glong(int c[]){ |
| 54 | + |
| 55 | + int i; |
| 56 | + |
| 57 | + cout<<"Enter the degree of longitude(west is negativ,form: 40 40 40 (means 40°40′40″)):"<<endl; |
| 58 | + |
| 59 | + for(i=0;i<3;i++){ |
| 60 | + |
| 61 | + cin>>c[i]; |
| 62 | + |
| 63 | + } |
| 64 | + |
| 65 | + } |
| 66 | + |
| 67 | + //输入经度 |
| 68 | + |
| 69 | + |
| 70 | +int leap_year(int year){ |
| 71 | + |
| 72 | + if(((year%400==0) || (year%100!=0) && (year%4==0))) return 1; |
| 73 | + |
| 74 | + else return 0; |
| 75 | + |
| 76 | + } |
| 77 | + |
| 78 | + //判断是否为闰年:若为闰年,返回1;若非闰年,返回0 |
| 79 | + |
| 80 | + |
| 81 | + int days(int year, int month, int date){ |
| 82 | + |
| 83 | + int i,a=0; |
| 84 | + |
| 85 | + for(i=2000;i<year;i++){ |
| 86 | + |
| 87 | + if(leap_year(i)) a=a+366; |
| 88 | + |
| 89 | + else a=a+365; |
| 90 | + |
| 91 | + } |
| 92 | + |
| 93 | + if(leap_year(year)){ |
| 94 | + |
| 95 | + for(i=0;i<month-1;i++){ |
| 96 | + |
| 97 | + a=a+days_of_month_2[i]; |
| 98 | + |
| 99 | + } |
| 100 | + |
| 101 | + } |
| 102 | + |
| 103 | + else { |
| 104 | + |
| 105 | + for(i=0;i<month-1;i++){ |
| 106 | + |
| 107 | + a=a+days_of_month_1[i]; |
| 108 | + |
| 109 | + } |
| 110 | + |
| 111 | + } |
| 112 | + |
| 113 | + a=a+date; |
| 114 | + |
| 115 | + return a; |
| 116 | + |
| 117 | + } |
| 118 | + |
| 119 | + //求从格林威治时间公元2000年1月1日到计算日天数days |
| 120 | + |
| 121 | + |
| 122 | + long double t_century(int days, long double UTo){ |
| 123 | + |
| 124 | + return ((long double)days+UTo/360)/36525; |
| 125 | + |
| 126 | + } |
| 127 | + |
| 128 | + //求格林威治时间公元2000年1月1日到计算日的世纪数t |
| 129 | + |
| 130 | + |
| 131 | + long double L_sun(long double t_century){ |
| 132 | + |
| 133 | + return (280.460+36000.770*t_century); |
| 134 | + |
| 135 | + } |
| 136 | + |
| 137 | + //求太阳的平黄径 |
| 138 | + |
| 139 | + |
| 140 | +long double G_sun(long double t_century){ |
| 141 | + |
| 142 | + return (357.528+35999.050*t_century); |
| 143 | + |
| 144 | + } |
| 145 | + |
| 146 | + //求太阳的平近点角 |
| 147 | + |
| 148 | + |
| 149 | +long double ecliptic_longitude(long double L_sun,long double G_sun){ |
| 150 | + |
| 151 | + return (L_sun+1.915*sin(G_sun*PI/180)+0.02*sin(2*G_sun*PI/180)); |
| 152 | + |
| 153 | + } |
| 154 | + |
| 155 | + //求黄道经度 |
| 156 | + |
| 157 | + |
| 158 | +long double earth_tilt(long double t_century){ |
| 159 | + |
| 160 | + return (23.4393-0.0130*t_century); |
| 161 | + |
| 162 | + } |
| 163 | + |
| 164 | + //求地球倾角 |
| 165 | + |
| 166 | + |
| 167 | +long double sun_deviation(long double earth_tilt, long double ecliptic_longitude){ |
| 168 | + |
| 169 | + return (180/PI*asin(sin(PI/180*earth_tilt)*sin(PI/180*ecliptic_longitude))); |
| 170 | + |
| 171 | + } |
| 172 | + |
| 173 | + //求太阳偏差 |
| 174 | + |
| 175 | + |
| 176 | +long double GHA(long double UTo, long double G_sun, long double ecliptic_longitude){ |
| 177 | + |
| 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)); |
| 179 | + |
| 180 | + } |
| 181 | + |
| 182 | + //求格林威治时间的太阳时间角GHA |
| 183 | + |
| 184 | + |
| 185 | + long double e(long double h, long double glat, long double sun_deviation){ |
| 186 | + |
| 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))); |
| 188 | + |
| 189 | + } |
| 190 | + |
| 191 | + //求修正值e |
| 192 | + |
| 193 | + |
| 194 | + long double UT_rise(long double UTo, long double GHA, long double glong, long double e){ |
| 195 | + |
| 196 | + return (UTo-(GHA+glong+e)); |
| 197 | + |
| 198 | + } |
| 199 | + |
| 200 | + //求日出时间 |
| 201 | + |
| 202 | + |
| 203 | +long double UT_set(long double UTo, long double GHA, long double glong, long double e){ |
| 204 | + |
| 205 | + return (UTo-(GHA+glong-e)); |
| 206 | + |
| 207 | + } |
| 208 | + |
| 209 | + //求日落时间 |
| 210 | + |
| 211 | + |
| 212 | +long double result_rise(long double UT, long double UTo, long double glong, long double glat, int year, int month, int date){ |
| 213 | + |
| 214 | + long double d; |
| 215 | + |
| 216 | + if(UT>=UTo) d=UT-UTo; |
| 217 | + |
| 218 | + else d=UTo-UT; |
| 219 | + |
| 220 | + if(d>=0.1) { |
| 221 | + |
| 222 | + UTo=UT; |
| 223 | + |
| 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)))))); |
| 225 | + |
| 226 | + result_rise(UT,UTo,glong,glat,year,month,date); |
| 227 | + |
| 228 | + } |
| 229 | + |
| 230 | + return UT; |
| 231 | + |
| 232 | + } |
| 233 | + |
| 234 | + //判断并返回结果(日出) |
| 235 | + |
| 236 | + |
| 237 | +long double result_set(long double UT, long double UTo, long double glong, long double glat, int year, int month, int date){ |
| 238 | + |
| 239 | + long double d; |
| 240 | + |
| 241 | + if(UT>=UTo) d=UT-UTo; |
| 242 | + |
| 243 | + else d=UTo-UT; |
| 244 | + |
| 245 | + if(d>=0.1){ |
| 246 | + |
| 247 | + UTo=UT; |
| 248 | + |
| 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)))))); |
| 250 | + |
| 251 | + result_set(UT,UTo,glong,glat,year,month,date); |
| 252 | + |
| 253 | + } |
| 254 | + |
| 255 | + return UT; |
| 256 | + |
| 257 | + } |
| 258 | + |
| 259 | + //判断并返回结果(日落) |
| 260 | + |
| 261 | + |
| 262 | +int Zone(long double glong){ |
| 263 | + |
| 264 | + if(glong>=0) return (int)((int)(glong/15.0)+1); |
| 265 | + |
| 266 | + else return (int)((int)(glong/15.0)-1); |
| 267 | + |
| 268 | + } |
| 269 | + |
| 270 | + //求时区 |
| 271 | + |
| 272 | + |
| 273 | +void output(long double rise, long double set, long double glong){ |
| 274 | + |
| 275 | + if((int)(60*(rise/15+Zone(glong)-(int)(rise/15+Zone(glong))))<10) |
| 276 | + |
| 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"; |
| 278 | + |
| 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"; |
| 280 | + |
| 281 | + if((int)(60*(set/15+Zone(glong)-(int)(set/15+Zone(glong))))<10) |
| 282 | + |
| 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"; |
| 284 | + |
| 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"; |
| 286 | + |
| 287 | + } |
| 288 | + |
| 289 | + //打印结果 |
| 290 | + |
| 291 | + |
| 292 | + |
| 293 | +int main(){ |
| 294 | + |
| 295 | + long double UTo=180.0; |
| 296 | + |
| 297 | + int year,month,date; |
| 298 | + |
| 299 | + long double glat,glong; |
| 300 | + |
| 301 | + int c[3]; |
| 302 | + |
| 303 | + input_date(c); |
| 304 | + |
| 305 | + year=c[0]; |
| 306 | + |
| 307 | + month=c[1]; |
| 308 | + |
| 309 | + date=c[2]; |
| 310 | + |
| 311 | + input_glat(c); |
| 312 | + |
| 313 | + glat=c[0]+c[1]/60+c[2]/3600; |
| 314 | + |
| 315 | + input_glong(c); |
| 316 | + |
| 317 | + glong=c[0]+c[1]/60+c[2]/3600; |
| 318 | + |
| 319 | + long double rise,set; |
| 320 | + |
| 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); |
| 322 | + |
| 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); |
| 324 | + |
| 325 | + output(rise,set,glong); |
| 326 | + |
| 327 | + system("pause"); |
| 328 | + |
| 329 | + return 0; |
| 330 | + |
| 331 | + } |
0 commit comments