|
发表于 2014-6-18 15:43:47
|
显示全部楼层
本帖最后由 dachun 于 2014-6-18 15:49 编辑
最近正好有个项目需要这个功能,从CSDN下载的C代码,经过验证计算日出日落时间和网络查询偏差不超5分钟,也算精确了,可以结贴了吧!呵呵
- #include <stdio.h>
- #include <stdlib.h>
- #include <string.h>
- #include "math.h"
- #include "time.h"
- #define M_PI 3.14
- static double RAD = 180.0 * 3600 /M_PI;
- static double richu;
- static double midDayTime;
- static double dawnTime;
- static double jd;
- static double wd;
- static double jd_degrees;
- static double jd_seconds;
- static double wd_degrees;
- static double wd_seconds;
- /*************************
- * 儒略日的计算
- *
- * @param y 年
- *
- * @param M 月
- *
- * @param d 日
- *
- * @param h 小时
- *
- * @param m 分
- *
- * @param s秒
- *
- * @return int
- ***************************/
- static double timeToDouble(int y, int M, double d)
- {
- // double A=0;
- double B=0;
- double jd=0;
- //设Y为给定年份,M为月份,D为该月日期(可以带小数)。
- //若M > 2,Y和M不变,若 M =1或2,以Y–1代Y,以M+12代M,换句话说,如果日期在1月或2月,则被看作是在前一年的13月或14月。
- //对格里高利历有 :A = INT(Y/100) B = 2 - A + INT(A/4)
- //对儒略历,取 B = 0
- //JD = INT(365.25(Y+4716))+INT(30.6001(M+1))+D+B-1524.5 (7.1)
- B=-13;
- jd=floor(365.25 * (y + 4716))+ floor(30.60001 * (M + 1)) +B+ d- 1524.5;
- return jd;
- }
- static void doubleToStr(double time,char *str)
- {
- double t;
- int h,m,s;
- t = time + 0.5;
- t = (t - (int) t) * 24;
- h = (int) t;
- t = (t - h) * 60;
- m = (int) t;
- t = (t - m) * 60;
- s = (int) t;
- sprintf(str,"%02d:%02d:%02d",h,m,s);
- }
- /****************************
- * @param t 儒略世纪数
- *
- * @return 太阳黄经
- *****************************/
- static double sunHJ(double t)
- {
- double j;
- t = t + (32.0 * (t + 1.8) * (t + 1.8) - 20) / 86400.0 / 36525.0;
- // 儒略世纪年数,力学时
- j = 48950621.66 + 6283319653.318 * t + 53 * t * t - 994 + 334166 *cos(4.669257 + 628.307585 * t) + 3489 * cos(4.6261 + 1256.61517 * t) + 2060.6 * cos(2.67823 + 628.307585 * t) * t;
- return (j / 10000000);
- }
- static double mod(double num1, double num2)
- {
- num2 = fabs(num2);
- // 只是取决于Num1的符号
- return num1 >= 0 ?(num1 - (floor(num1 / num2)) * num2 ): ((floor(fabs(num1) / num2)) * num2 - fabs(num1));
- }
- /********************************
- * 保证角度∈(-π,π)
- *
- * @param ag
- * @return ag
- ***********************************/
- static double degree(double ag)
- {
- ag = mod(ag, 2 * M_PI);
- if(ag<=-M_PI){
- ag=ag+2*M_PI;
- }
- else if(ag>M_PI){
- ag=ag-2*M_PI;
- }
- return ag;
- }
- /***********************************
- *
- * @param date 儒略日平午
- *
- * @param lo 地理经度
- *
- * @param la 地理纬度
- *
- * @param tz 时区
- *
- * @return 太阳升起时间
- *************************************/
- double sunRiseTime(double date, double lo, double la, double tz)
- {
- double t,j,sinJ,cosJ,gst,E,a,D,cosH0,cosH1,H0,H1,H;
- date = date - tz;
- // 太阳黄经以及它的正余弦值
- t = date / 36525;
- j = sunHJ(t);
- // 太阳黄经以及它的正余弦值
- sinJ = sin(j);
- cosJ = cos(j);
- // 其中2*M_PI*(0.7790572732640 + 1.00273781191135448*jd)恒星时(子午圈位置)
- gst = 2 * M_PI * (0.779057273264 + 1.00273781191135 * date) + (0.014506 + 4612.15739966 * t + 1.39667721 * t * t) / RAD;
- E = (84381.406 - 46.836769 * t) / RAD; // 黄赤交角
- a = atan2(sinJ * cos(E), cosJ);// '太阳赤经
- D = asin(sin(E) * sinJ); // 太阳赤纬
- cosH0 = (sin(-50 * 60 / RAD) - sin(la) * sin(D)) / (cos(la) * cos(D)); // 日出的时角计算,地平线下50分
- cosH1 = (sin(-6 * 3600 / RAD) - sin(la) * sin(D)) / (cos(la) * cos(D)); // 天亮的时角计算,地平线下6度,若为航海请改为地平线下12度
- // 严格应当区分极昼极夜,本程序不算
- if (cosH0 >= 1 || cosH0 <= -1){
- return -0.5;// 极昼
- }
- H0 = -acos(cosH0); // 升点时角(日出)若去掉负号 就是降点时角,也可以利用中天和升点计算
- H1 = -acos(cosH1);
- H = gst - lo - a; // 太阳时角
- midDayTime = date - degree(H) / M_PI / 2 + tz; // 中天时间
- dawnTime = date - degree(H - H1) / M_PI / 2 + tz;// 天亮时间
- return date - degree(H - H0) / M_PI / 2 + tz; // 日出时间,函数返回值
- }
- int main()
- {
- time_t curtime;
- struct tm * timeinfo;
- int i,j,k;
- char timestr[20];
- int year;
- int month;
- double day;
- int hour;
- int min;
- int sec;
- int tz;
- FILE *fp;
- //上海东经121度29分,北纬31度14分
- jd_degrees=121;
- jd_seconds=28;
- wd_degrees=31;
- wd_seconds=14;
- //乌鲁木齐东经87度35分,北纬43度48分
- //jd_degrees=87;
- //jd_seconds=35;
- //wd_degrees=43;
- //wd_seconds=48;
- //广州 E113°16' N23°06'
- //jd_degrees=113;
- //jd_seconds=16;
- //wd_degrees=23;
- //wd_seconds=6;
- tz=8;
- //得到当前时间
- time( &curtime );
- timeinfo=localtime(&curtime);
- printf ( "The current date/time is: %s", asctime(timeinfo));
- //step 1
- jd=-(jd_degrees+jd_seconds/60)/180*M_PI;
- wd=(wd_degrees+wd_seconds/60)/180*M_PI;
- printf("jd=%f\r\n",jd);
- printf("wd=%f\r\n",wd);
- //step 2
- year=timeinfo->tm_year+1900;
- month=timeinfo->tm_mon+1;
- day=timeinfo->tm_mday;
- hour=timeinfo->tm_hour;
- min=timeinfo->tm_min;
- sec=timeinfo->tm_sec;
- richu = timeToDouble(year,month,day) - 2451544.5;
- printf ( "mjd:%f\r\n", richu);
- printf("year=%d,month=%d,day=%f\r\n",year,month,day);
- for (i = 0; i < 10; i++){
- richu = sunRiseTime(richu, jd, wd, tz/24.0);// 逐步逼近法算10次
- }
- doubleToStr(richu,timestr);
- printf("日出时间 %s\r\n",timestr);
- doubleToStr(midDayTime + midDayTime - richu, timestr);
- printf("日落时间 %s\r\n",timestr);
- printf ( "richu:%f,midDayTime:%f\r\n", richu,midDayTime);
- doubleToStr(midDayTime,timestr);
- printf("中天时间 %s\r\n",timestr);
- doubleToStr(dawnTime,timestr);
- printf("天亮时间 %s\r\n",timestr);
- doubleToStr(midDayTime + midDayTime - dawnTime,timestr);
- printf("天黑时间 %s\r\n",timestr);
- fp=fopen("richu.txt","w");
- if(fp==NULL){
- printf("file open error\r\n");
- }
- fclose(fp);
- printf("Hello world!\n");
- return 0;
- }
复制代码 |
|