地球日升日落时间算法原理及源代码

地球日升日落时间算法原理及源代码 篇一:地球日升日落时间算法原理及源代码 地球日升日落时间算法原理及源代码 * 本算法已考虑太阳的椭圆运动、轨道的长期变化、光行差、黄赤交角的长期变化 * 本算法适用于理想气温、气压、地形 * 本算法的计算结果与《2008年中国天文年历》的计算结果完全相符 在地球上的人们，看到每天太阳绕地球自转轴做圆周运动(周日运动)，而且不会觉查地球在运动。有一种古代的仪器叫日晷，也称作日规。我们就称他为“日圆规”吧。如果日圆规平行于赤道面，就可以生动的展示出太阳的这种动动。 太阳的这种运动，在站心天球中扫过的轨迹是基本上是平行于赤道的一个圆(实际是一个不封闭的准圆)。扫过一圈所用的时间约为24小时。在这个圆圈中，一般可以找到两个特殊点:这两个点正好在地平圈上。太阳经过这两个点时间，就是日出(升点)日没(降点)时间。考虑大气折射因素，这两个点不妨定义在:在地平圈下方50角分。为了述叙方便，今设升点为A，降点为B，日出时间为Ta，日落时间为Tb。 不好意思，昨天没有画图，现在补上 计算步骤: 1、计算出 T 时刻太阳赤道坐标，即赤经α和赤纬δ，最初的T 可选择当日的12点，顺便算出太阳时角Ho，即子午圈相对太阳的经度。严格说是午半圈相对太阳的角度为时角Ho 2、当赤纬δ已知，就用δ作为A点和B点的赤纬的估值。由于T与Ta(或T与Tb)的时间差最多也就几个小时，所以太阳赤纬δ基本不变，所以用δ作为A点或B点的初始估值，也是很准确的。 3、当A点的的赤纬已知，那么经过A点的平行圈就是确定的。反过来，在这个平行圈上找出地平圈下方50角分的点，这个点就是A或B。如何计算A点和B点的坐标，我们可以使用球面坐标得到，这个公式与黄道转赤道的公式非常相似。A(或B)的坐标一般描述为“赤纬δ，时角H”，所谓A点的时角，就是子午圈相对于A的经度。这里用到球面公式，读者不妨自行推导一下，这一推导虽然有点麻烦，但确实只需高中的数学知识就可以实现。 4、太阳的时角(即Ho)与点A的时角(H)的差值ΔH = Ho - H就是太阳在平行圈上跑过头的角度，为了追朔日出时间(太阳经过A的时间)，太阳就应该退回 ΔH角度，我们又知道，太阳在平行圈上运行的速度大约是360度/天，所以，所以时间就应退回ΔH/(360度/天)，即太阳经过A点的时间为 Ta = T - ΔH/360 5、以Ta代替为T做为初始值，重新按照以上四个步骤再计算一遍，就可以精确得到Ta 附1:子午圈的经度使用恒星时计算，那么太阳时角就是 GST-L-α，式中GST-L是本地恒星时，所以GST-L-α就是本地子午圈相对太阳的时角，本文简称太阳时角。 附2:赤经α和赤纬δ应用当日分点坐标(视赤经视赤纬)，这样，与恒星时的计算使用的坐标系统才是相同的。 附3:如果计算月出与月没: 如果没有特别的精度要求，与日出日没的算法一模一样。月亮有时比较大，有时比较小，所以给计算带来一些麻烦，但是，如果你不要求精确到1秒，就把月亮看成与太阳一样大。另外，气温、气压、地形，均对日出有影响。还应注意到，月亮离我们很近，站在地心与站在地面看月亮，月亮的位置是不同的，站在地面看月亮，月亮的高度角下降，所以不能取50角分，可参考《天文算法》取值。 附4:以上算法中的第5点，实际上是逐步逼近的思想方法。这类思想方法与高等数学有关，相关著作请参考《数值方法》。 本算法应用于《寿星万年历》。 <script language=javascript> /****************** 太阳升起时间计算。xjw01，2009.6.21 ******************/ //===============数学工具================== function sin(x) { return Math.sin(x);} function cos(x) { return Math.cos(x);} function rad2rrad(v){//对超过-PI到PI的角度转为-PI到PI v=v % (2*Math.PI); if(v<=-Math.PI) return v+2*Math.PI; if(v>Math.PI) return v-2*Math.PI; return v; } function timeStr(jd){ //取儒略日数中的时间 var T=jd+0.5; T=(T-Math.floor(T))*24; var h,m,s; h=Math.floor(T); T=(T-h)*60; m=Math.floor(T); T=(T-m)*60; s=Math.floor(T+0.5); return h+':'+m+':'+s; } var rad = 180*3600/Math.PI; //===============太阳升起================== function sunLat(t){ //太阳黄经,t为儒略世纪数UT t += (32*(t+1.8)*(t+1.8)-20)/86400/36525; //儒略世纪年数,力学 时var 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; } function sheng(jd,L,fa,TZ){ //jd儒略日平午，L地理经度，fa地理 纬度，TZ时区，返回太阳升起时间 jd-=TZ; var T=jd/36525; var J = sunLat(T), sinJ=sin(J), cosJ=cos(J); //太阳黄经以及它的正 余弦值 var gst = 2*Math.PI*(0.7790572732640 + 1.00273781191135448*jd) //恒星时(子午圈位置) + (0.014506 + 4612.15739966*T + 1.39667721*T*T)/rad; var E = (84381.4060 -46.836769*T)/rad; //黄赤交角 var A=Math.atan2( sinJ*cos(E), cosJ); //太阳赤经 var D=Math.asin ( sin(E)*sinJ );//太阳赤纬 var cosH0 = (sin(-50*60/rad) - sin(fa)*sin(D) ) / ( cos(fa)*cos(D) ); //日出的时角计算 if(cosH0>=1||cosH0<=-1) return 0; var H0 = -Math.acos( cosH0 ); //升点时角(日出) var H = gst-L-A; //太阳时角 return jd - rad2rrad(H-H0)/6.28+TZ; } //===============测试代码================== //北京的地标 -116?23' 纬度 + 39?54' L = -(116+23/60)/180*Math.PI; fa= +( 39+54/60)/180*Math.PI; for(i=0;i<30;i++){ //计算30天测试 T=i; //T从2000年1月1.5日起算 T = sheng(T,L,fa,8/24); T = sheng(T,L,fa,8/24); document.write('T=' + i +' '+ timeStr(T) + '<br>'); } </script> 以上程序，2000年1月1日(T=0)至2000年1月30日(T=29)的计算结果，北京: T=0 7:35:58 T=1 7:36:6 T=2 7:36:13 T=3 7:36:17 T=4 7:36:19 T=5 7:36:19 T=6 7:36:17 T=7 7:36:12 T=8 7:36:6 T=9 7:35:57 T=10 7:35:45 T=11 7:35:32 T=12 7:35:17 T=13 7:34:59 T=14 7:34:39 T=15 7:34:17 T=16 7:33:53 T=17 7:33:26 T=18 7:32:58 T=19 7:32:28 T=20 7:31:55 T=21 7:31:20 T=22 7:30:44 T=23 7:30:5 T=24 7:29:25 T=25 7:28:42 T=26 7:27:58 T=27 7:27:12 T=28 7:26:24 T=29 7:25:34 篇二:日出日落时间计算程序(C语言) //日出日落时间计算C语言程序 #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){ 篇三:怎样用经纬度计算日出日落的时间 怎样用经纬度计算日出日落的时间 下面是一种随经纬度变化的日出日落时间计算方法,我成功运用在一智能路灯控制器中,希望对需要的朋友有帮助。 已知:日出日落时太阳的位置h,-0.833?,要计算地的地理位置,经度Long,纬度G1at,时区zone,UTo为上次计算的日出日落时间,第一次计算时UTo,180?。 (1)先计算出从格林威治时间公元2000年1月1日到计算日天数days; (2)计算从格林威治时间公元2000年1月1日到计算日的世纪数t, 则 t,(days+UTo,360),36525; (3)计算太阳的平黄径 L=280.460+36000.770×t; (4)计算太阳的平近点角 G,357.528+35999.050×t (5)计算太阳的黄道经度 λ,L+1.915×sinG+0.020xsin(2G); (6)计算地球的倾角 ε,23.4393-0.0130×t; (7)计算太阳的偏差 δ,arcsin(sinε×sinλ); (8)计算格林威治时间的太阳时间角GHA: GHA=UTo-180-1.915×sinG-0.020×sin(2G)+2.466×sin(2λ)-0.053×sin(4λ) (9)计算修正值e: e=arcos{[sinh-sin(Glat)sin(δ)]/cos(Glat)cos(δ)} (10)计算新的日出日落时间 UT,UTo-(GHA+Long?e); 其中“+”示计算日出时间,“-”表示计算日落时间; (11)比较UTo和UT之差的绝对值,如果大于0.1?即0.007小时,把UT作为新的日出日落时间值,重新从第(2)步开始进行迭代计算,如果UTo和UT之差的绝对值小于0.007小时,则UT即为所求的格林威治日出日落时间; (12)上面的计算以度为单位,即180?=12小时,因此需要转化为以小时表示的时间,再加上所在的时区数Zone,即要计算地的日出日 落时间为 T=UT/15+Zone 上面的计算日出日落时间方法适用于小于北纬60?和南纬60?之间的区域,如果计算位置为西半球时,经度Long为负数。