日出日落时间计算程序(C语言)

日出日落时间计算程序(C语言)//日出日落时间计算C语言程序#define PI 3.1415926#include<math.h>#include<iostream>using namespace std;intdays_of_month_1[]={31,28,31,30,31,30,31,3 1,30,31,30,31};intdays_of_month_2[]={31,29,31,30,31,30,31,3 1,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；若非闰年,返回0int 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日到计算日天数dayslong double t_century(int days, long double UTo){return ((long double)days+UTo/360)/36525;}//求格林威治时间公元2000年1月1日到计算日的世纪数tlong 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(P I/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_lon gitude*PI/180)-0.053*sin(4*ecliptic_longi tude*PI/180));}//求格林威治时间的太阳时间角GHAlong double e(long double h, long double glat, long double sun_deviation){return180/PI*acos((sin(h*PI/180)-sin(glat*PI/180)*sin(sun_deviation*PI/180))/(cos(glat*P I/180)*cos(sun_deviation*PI/180)));}//求修正值elong 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(da ys(year,month,date),UTo)),ecliptic_longit ude(L_sun(t_century(days(year,month,date) ,UTo)),G_sun(t_century(days(year,month,da te),UTo)))),glong,e(h,glat,sun_deviation( earth_tilt(t_century(days(year,month,date ),UTo)),ecliptic_longitude(L_sun(t_centur y(days(year,month,date),UTo)),G_sun(t_cen tury(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(day s(year,month,date),UTo)),ecliptic_longitu de(L_sun(t_century(days(year,month,date), UTo)),G_sun(t_century(days(year,month,dat e),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_cent ury(days(year,month,date),UTo))))));result_set(UT,UTo,glong,glat,year,month,d ate);}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)(ri se/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+Z one(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+Zo ne(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_su n(t_century(days(year,month,date),UTo)),e cliptic_longitude(L_sun(t_century(days(ye ar,month,date),UTo)),G_sun(t_century(days (year,month,date),UTo)))),glong,e(h,glat, sun_deviation(earth_tilt(t_century(days(y ear,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)),ecli ptic_longitude(L_sun(t_century(days(year, month,date),UTo)),G_sun(t_century(days(ye ar,month,date),UTo)))),glong,e(h,glat,sun _deviation(earth_tilt(t_century(days(year ,month,date),UTo)),ecliptic_longitude(L_s un(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;}。