儒略日计算

Python014

儒略日计算,第1张

求出给定年(I),月(J),日(K)的儒略日:

儒略日= K - 32075 + 1461 * (I + 4800 + (J-14)/12)/4+367*(J-2-(J-14)/12*12)/12-3*((I+4900+(J-14)/12)/100)/4

例如1986年12月25日的儒略日(I=1986,J=12,K=25)=2446791

儒略日(Julian day)是指由公元前4713年1月1日,协调世界时中午12时开始所经过的天数,多为天文学家采用,用以作为天文学的单一历法,把不同历法的年表统一起来。

儒略日是一种不用年月的长期纪日法,简写为JD。是由法国学者Joseph Justus Scliger(1540–1609)在1583年所创,这名称是为了纪念他的父亲——意大利学者Julius Caesar Scaliger(1484–1558)。

儒略日的起点订在公元前4713年(天文学上记为 -4712 年)1月1日格林威治时间平午(世界时12:00),即JD 0 指定为 4713 B.C. 1月1日12:00 UT到4713 B.C. 1月2日12:00 UT的24小时。每一天赋予了一个唯一的数字,顺数而下,如:1996年1月1日12:00:00的儒略日是2450084。这个日期是考虑了太阳、月亮的运行周期,以及当时收税的间隔而订出来的。Joseph Scliger定义儒略周期为7980年,是因28、19、15的最小公倍数为28×19×15=7980。其中:

28年为一太阳周期(solar cycle),经过一太阳周期,则星期的日序与月的日序会重复。

19年为一太阴周期,或称默冬章(Metonic cycle),因235朔望月=19回归年,经过一太阴周期则阴历月年的日序重复。

15年为一小纪(indiction cycle),此为罗马皇帝君士坦丁(Constantine)所颁,每15年评定财产价值以供课税,成为古罗马用的一个纪元单位,

故以7980年为一儒略周期,而所选的起点公元前4713年,则是这三个循环周期同时开始的最近年份。

以儒略日计日是为方便计算年代相隔久远或不同历法的两事件所间隔的日数。

由于儒略日数字位数太多,国际天文学联合会于1973年采用简化儒略日(MJD),其定义为 MJD = JD– 2400000.5。MJD相应的起点是1858年11月17日世界时0时。

月球位置计算

月球的位置计算是所有天体计算里比较复杂的部分。往往为了提高一点计算的精确度,必须使用数百个扰动修正项,目前这些修正项在许多专业的天文书籍里都找得到。由於月球每天在天空移动的数度很快,因此计算时间的正确性要求的比较多。之前计算太阳与行星时,都采用地方平均太阳时,但是计算月球位置必须采用力学时来减低误差。这里提供的方法,月球位置的误差约小於1/15度。虽然在许多观测上,这个精密度可能稍嫌不足,但是却能够简单地推算月球的位置也算是实用的方法。

虽然月球绕行地球公转,但是严格说起来,月球应该看成绕太阳运行。由於受到地球的扰动,所以月球以蛇行的方式环绕太阳。月球的轨道面会产生很大的位移,这也是他与其他行星最大的不同点。因为迅速移动的轨道面,造成实际计算非常不易。大致来说,月球相对於黄道面的轨道进动周期为18.6年(如同地球的岁差,我们又称为升交点进动)。月球的近日点也会周期性的运动,主要的项次约为8.85年。另外月球绕地运行的轨道半径每年约增加数公分之多。

计算时间

时间 JD1= JD-JD2000.0

说明:这里的JD是使用力学时(ET)来计算,相关的资料请参考「时间」的网页。

距离2000年的儒略世纪 T=(JD-JD2000)/ 36525

太阳的轨道资料

轨道周期 P=1.00004回归年

距离春分点平均离角 ε=280.466457+0.985647358 JD1+0.000304 T2

近地点离角 ω=282.937348+0.00004707624 JD1+0.0004569 T2

轨道偏心率 E=0.01670862-0.00004204 T

轨道半径 A=1.00001161 AU

计算结果

近日点平均离角 M=ε-ω

利用级数展开的克卜勒方程式(Kepler’s Equation)求得真实的近日点离角

ν=M+360 E sin(M)/π+900 E2 sin(2M)/4π-180E3sin(M)/ 4π

太阳的黄道经度 λ=ν+ω

月球的轨道资料(以黄道面来定义)

月球的平均经度 L=218.316646+13.17639647564 JD1-0.0014664 T2

月球的平均近日点离角 Mm=L-(83.353243+0.11140352394 JD1-0.0103217 T2)

升交点引数 N=125.044555-0.05295376277 JD1+0.0020756 T2

月面轨道倾角 I=5.15668983

扰动修正项

Ev=1.2739 sin( 2L-2λ-Mm )

Ae=0.1858 sin( M )

A3=0.37 sin( M )

修正月球的平均近日点离角

Mm'=Mm+Ev-Ae-A3

Ec=6.2886 sin( Mm' )

A4=0.214 sin( 2Mm' )

修正月球的平均经度

L'=L+Ev+Ec-Ae+A4

V=0.6583 sin( 2L'-2λ)

再次修正月球的平均经度 L''=L'+V

修正升交点引数 N'=N-0.16 sin( M )

计算月球的黄道经度与纬度

sin(θ)=sin( L''-N' ) cos( I )

cos(θ)=cos( L''-N' )

判断θ的象限后,求得θ。

月球的黄道经度 λm=θ+N'

月球的黄道纬度 βm=sin-1[ sin( L''-N' ) sin( I ) ]

若需要月球的赤道经度与纬度,可参考座标系统的变换部分。将黄道座标转成赤道座标。或者更进一步转成地平方为角。

范例

2000年9月1日 AM10:30 月球的黄道经纬度 、赤道经纬度 与地平方为角

查得当年力学时与平均太阳时的差值 ΔT=64秒

ET=UT+64秒

计算ET儒略日得到 JD=2451788.60492

JD1=243.60492

儒略世纪 T=0.00667

太阳轨道资料

ε=520.492123

ω=282.948812

E=0.01671

太阳黄道位置计算结果

M=237.543311

λ=158.8947 β=0

月球的轨道资料(以黄道面来定义)

月球的平均经度 L=188.3759761

月球的平均近日点离角 Mm=77.88239044

升交点引数 N=112.14385644

月面轨道倾角 I=5.15668983

扰动修正项

Ev=-0.41717013

Ae=-0.15695074

A3=-0.31254992

修正月球的平均近日点离角

Mm'=77.93472097

Ec=6.14968536

A4=0.08748689

修正月球的平均经度

L'=194.35292896

V=0.62138281

再次修正月球的平均经度 L''=194.97431177

修正升交点引数 N'=112.27901316

计算月球的黄道经度与纬度

sin(θ)=0.98786949

cos(θ)=0.127146

判断θ的象限后,求得θ=82.66593614

月球的黄道经度 λm=194.9449493

月球的黄道纬度 βm=5.11472622

经过黄道座标与赤道座标系统后

α=13h 2m56s

δ=-1°10' 7''

经过地平座标转换后

A=107°

h=30°

白天也可以看到月球!

若是经过多项的摄动计算后,月球的位置会更精确,它的赤道座标位置:

α=13h 1m24s

δ=-1° 2' 31''

计算月球位置的c语言程序

////////////////////////////////////////////////////////////////////////////////////////////////////////

// 名称:月球位置计算

// 作者:胡铂(http://hubble.lamost.org)

// 日期:2004-09-29

// 说明:根据北京天文同好会提供的《Astronomy Algrithms》翻译版实现,并得到了

// zjuglr的帮助。

//

///////////////////////////////////////////////////////////////////////////////////////////////////////

#include "math.h"

#include "stdio.h"

#define DE 3.141592654/180

//////////////////////////////////计算儒略日历书时//////////////////////////////////////////////////////////////

float jde(int Y,int M,int D,int hour,int min,int sec)

{

int f,g

double mid1,mid2,J,JDE,A

if(M>=3)

{

f=Y

g=M

}

if(M==1||M==2)

{

f=Y-1

g=M+12

}

mid1=floor(365.25*f)

mid2=floor(30.6001*(g+1))

A=2-floor(f/100)+floor(f/400)

J=mid1+mid2+D+A+1720994.5

JDE=J+hour/24+min/1440+sec/86400

return JDE

////////////////////////////////////////////////////////////////////////////////////////////////////////////////

}

void main(void)

{

/////////////////////////////////////////变量定义///////////////////////////////////////////////////////////////////

int i,year,month,day,hour,min,sec

double JDE,T,L1,D,M,M1,F,A1,A2,A3,E,SUML,lamda,SUMB,beta,SUMR,SIN1,SIN2,COS1,Dist

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////数据/////////////////////////////////////////////////////////////////

static double La[60]={0,2,2,0,0,0,2,2,2,2,0,1,0,2,0,0,4,0,4,2,2,1,1,2,2,4,2,0,2,2,1,2,0,0,2,2,2,4,0,3,2,4,0,2,2,2,4,0,4,1,2,0,1,3,4,2,0,1,2,2}

static double Lb[60]={0,0,0,0,1,0,0,-1,0,-1,1,0,1,0,0,0,0,0,0,1,1,0,1,-1,0,0,0,1,0,-1,0,-2,1,2,-2,0,0,-1,0,0,1,-1,2,2,1,-1,0,0,-1,0,1,0,1,0,0,-1,2,1,0,0}

static double Lc[60]={1,-1,0,2,0,0,-2,-1,1,0,-1,0,1,0,1,1,-1,3,-2,-1,0,-1,0,1,2,0,-3,-2,-1,-2,1,0,2,0,-1,1,0,-1,2,-1,1,-2,-1,-1,-2,0,1,4,0,-2,0,2,1,-2,-3,2,1,-1,3,-1}

static double Ld[60]={0,0,0,0,0,2,0,0,0,0,0,0,0,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0,-2,2,0,2,0,0,0,0,0,0,-2,0,0,0,0,-2,-2,0,0,0,0,0,0,0,-2}

static double Sl[60]={6288774,1274027,658314,213618,-185116,-114332,58793,57066,53322,45758,-40923,-34720,-30383,15327,-12528,10980,10675,10034,8548,-7888,-6766,-5163,4987,4036,3994,3861,3665,-2689,-2602,2390,-2348,2236,-2120,-2069,2048,-1773,-1595,1215,-1110,-892,-810,759,-713,-700,691,596,549,537,520,-487,-399,-381,351,-340,330,327,-323,299,294,0}

static double Sr[60]={-20905355,-3699111,-2955968,-569925,48888,-3149,246158,-152138,-170733,-204586,-129620,108743,104755,10321,0,79661,-34782,-23210,-21636,24208,30824,-8379,-16675,-12831,-10445,-11650,14403,-7003,0,10056,6322,-9884,5751,0,-4950,4130,0,-3958,0,3258,2616,0,-2117,2354,0,0,0,0,0,0,0,-4421,0,0,0,0,1165,0,0,8752}

static double Sb[60]={5128122,280602,277693,173237,55413,46271,32573,17198,9266,8822,8216,4324,4200,-3359,2463,2211,2065,-1870,1828,-1794,-1749,-1565,-1491,-1475,-1410,-1344,-1335,1107,1021,833,777,671,607,596,491,-451,439,422,421,-366,-351,331,315,302,-283,-229,223,223,-220,-220,-185,181,-177,176,166,-164,132,-119,115,107}

static double Ba[60]={0,0,0,2,2,2,2,0,2,0,2,2,2,2,2,2,2,0,4,0,0,0,1,0,0,0,1,0,4,4,0,4,2,2,2,2,0,2,2,2,2,4,2,2,0,2,1,1,0,2,1,2,0,4,4,1,4,1,4,2}

static double Bb[60]={0,0,0,0,0,0,0,0,0,0,-1,0,0,1,-1,-1,-1,1,0,1,0,1,0,1,1,1,0,0,0,0,0,0,0,0,-1,0,0,0,0,1,1,0,-1,-2,0,1,1,1,1,1,0,-1,1,0,-1,0,0,0,-1,-2}

static double Bc[60]={0,1,1,0,-1,-1,0,2,1,2,0,-2,1,0,-1,0,-1,-1,-1,0,0,-1,0,1,1,0,0,3,0,-1,1,-2,0,2,1,-2,3,2,-3,-1,0,0,1,0,1,1,0,0,-2,-1,1,-2,2,-2,-1,1,1,-1,0,0}

static double Bd[60]={1,1,-1,-1,1,-1,1,1,-1,-1,-1,-1,1,-1,1,1,-1,-1,-1,1,3,1,1,1,-1,-1,-1,1,-1,1,-3,1,-3,-1,-1,1,-1,1,-1,1,1,1,1,-1,3,-1,-1,1,-1,-1,1,-1,1,-1,-1,-1,-1,-1,-1,1}

///////////////////////////////////////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////计算日期和时间///////////////////////////////////////////////////////////

year=2004

month=9

day=4

hour=0

min=0

sec=0

//////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////计算时间日期的儒略日历书时//////////////////////////////////////////////////

JDE=jde(year,month,day,hour,min,sec)

/////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////计算自J2000.0开始的儒略世纪数///////////////////////////////////////

T=(JDE-2451545)/36525

/////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////

L1=218.3164477+481267.88123421*T-0.0015786*T*T+T*T*T/538841-T*T*T*T/65194000//月亮的平黄经

D=297.8501921+445267.1114034*T-0.0018819*T*T+T*T*T/545868-T*T*T*T/113065000//月亮的平均太阳距角

M=357.5291092+35999.0502909*T-0.0001536*T*T+T*T*T/24490000//太阳的平近点角

M1=134.9633964+477198.8675055*T+0.0087414*T*T+T*T*T/69699-T*T*T*T/14712000//月亮的平近点角

F=93.2720950+483202.0175233*T-0.0036539*T*T-T*T*T/3526000+T*T*T*T/863310000//月亮的黄纬参量(由升交点起算的月球平均距离)

A1=119.75+131.849*T//金星的摄动

A2=53.09+479264.290*T//木星的摄动

A3=313.45+481266.484*T

E=1-0.002516*T-0.0000074*T*T//计算反映地球轨道偏心率变化

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////计算月球地心黄经周期项;////////////////////////////////////////////////

SUML=0

for(i=0i<=59i++)

{

SIN1=La[i]*D+Lb[i]*M+Lc[i]*M1+Ld[i]*F

SUML=SUML+Sl[i]*0.000001*sin(SIN1*DE)*pow(E,fabs(Lb[i]))

}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////

///////////////////////////////计算月球地心黄经///////////////////////////////////////////////////////////////

lamda=L1+SUML+(3958*sin(A1*DE)+1962*sin((L1-F)*DE)+318*sin(A2*DE))/1000000

lamda=fmod(lamda,360)

//////////////////////////////////////////////////////////////////////////////////////////////////////////////

/////////////////////////////计算月球地心黄纬周期项///////////////////////////////////////////////////////////

SUMB=0

for(i=0i<=59i++)

{

SIN2=Ba[i]*D+Bb[i]*M+Bc[i]*M1+Bd[i]*F

SUMB=SUMB+Sb[i]*0.000001*sin(SIN2*DE)*pow(E,fabs(Lb[i]))

}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////计算月球地心黄纬////////////////////////////////////////////////////

beta=SUMB+(-2235*sin(L1*DE) //(modified)

+382*sin(A3*DE)+175*sin((A1-F)*DE)

+175*sin((A1+F)*DE)+127*sin((L1-M1)*DE)

-115*sin((L1+M1)*DE))/1000000

//////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////计算月球地心距离周期项////////////////////////////////////////////////

SUMR=0

for(i=0i<=59i++)

{

COS1=La[i]*D+Lb[i]*M+Lc[i]*M1+Ld[i]*F

SUMR=SUMR+Sr[i]*0.001*cos(COS1*DE)*pow(E,fabs(Lb[i]))

}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////计算月球地心距离///////////////////////////////////////////////////////

//Dist=385000.56+SUMR/1000

Dist=385000.56+SUMR

////////////////////////////////////////////////////////////////////////////////////////////////////////////

printf("\n%d-%d-%d %d:%d:%d\n",year,month,day,hour,min,sec)

printf("JDE=%f\n",JDE)

printf("T=%f\n",T)

printf("L1=%f\n",L1)

printf("D=%f\n",D)

printf("M=%f\n",M)

printf("M1=%f\n",M1)

printf("F=%f\n",F)

printf("A1=%f\n",A1)

printf("A2=%f\n",A2)

printf("A3=%f\n",A3)

printf("E=%f\n",E)

printf("SUML=%f\n",SUML)

printf("SUMB=%f\n",SUMB)

printf("SUMR=%f\n",SUMR)

printf("lamda=%f\n beta=%f\n Dist=%f\n",lamda,beta,Dist)

可以吗