在写这篇文章之前的5个小时之内我并没哟把程序和历法联系起来的打算。但是现在我竟然开始写了。
这要源于我的一个同学问我的一个问题。
问题大概是这样的:把现在的一个日期转换为GPS周的某一天。
当时同学在电话中提到GPS周,GPS天,儒略日等等词语。前面两个我还勉强可以猜出是什么意思(毕竟我写过两年GPS方面的程序)。但是儒略日就不知道了。所以闲暇时Baidu了一下。在此记录算是学习笔记吧
先介绍几个概念
1、儒略历法
(英语:Julian calendar)是现今国际通用的公历的前身。西方国家16世纪大多采用它。公元前46年,罗马统帅盖厄斯·儒略·凯撒在埃及亚历山大的希腊数学家兼天文学家索西琴尼的帮助下制订的,并在公元前46年1月1日起执行实行,取代旧罗马历法的一种历法。所以人们就把这一历法称为《儒略历》。
2、格里历法
格里历,即公历或格列高利历,是现行国际通行的历法,属于阳历的一种,通称阳历,其前身是奥古斯都历,而奥古斯都历的前身是儒略历。其历年为一个回归年(365.2425日),划分为12个历月。是教皇格里高利13世(也译格雷果里)在公元1582年改革儒略历制定的历法,儒略历的回归年为365.25,与实际的回归365.2422相差甚多,当时儒略历和地球实际位置的误差已达14天,格里历将误差纠正,确定所有整数世纪年除了可被400整除的外一律不设闰年,同时规定1582年10月4日之后的那天为1582年10月15日,但原有星期不变。新颁布的历法理论上可以达到两万年内误差不超过一天,但由于地球自转的变化,实际到公元4909年误差就可达一天。与儒略历相比,公历是新制定的历法,所以有时候,包括中华民国《中国国家标准》CNS 7648《数据元及交换格式–信息交换–日期及时间的表示法》,又称新历。
3、儒略日
从天文计算的角度,只要以日为单位连续计时就可以了,这就是儒略日。年月的设置更多地是为了生活生产的需要,对于天文计算并不是必须的。
现在广泛使用的是法国学者史伽利日(Joseph Justus Scaliger,1540-1609)提出的儒略日系统。之所以叫做儒略日,与上面讲过的儒略历并不相干,而是因为史伽利日的父亲,意大利学者Julius Caesar Scaliger(1484-1558)与颁布实施儒略历的罗马统治者儒略•恺撒同名。系统以公元前4713年1月1日正午为起点,向后连续计日,简记为JD。积累到现在,已经是一个很大的数字。例如2000年1月1日地球力学时12时的儒略日记法就是JD 2451545.0,这是一个很重要的时刻,特别记为 J 2000.0。由于儒略日数字位数太多,国际天文学联合会于1973 年采用简化儒略日(MJD),其定义为 MJD = JD - 2400000.5。MJD 相应的起点是1858 年11月17日世界时0时。
那么为什么儒略日的起始时间定为公元前4713年1月1日正午呢?
原因是这样的:原来史伽利日构造这个系统时考虑了三种周期:阳历日期与星期会合的28年 (365.25×4×7日)周期,阳历与阴历会合的19年(29.53×(12×19+7)≈365.25×19=6939.75日,这个数字等到我们后 面讲农历的时候再来解释)周期,以及罗马政府为征税登记财产的15年周期。取3者的最小公倍数7,980年为儒略周期,然后向上推算,得到公元前4713 年1月1日为这三个周期同时开始的历元。换句话说,这一天既是公历元旦,又是农历初一和星期日,还是罗马政府(假如有的话)登记财产的日子。这就是儒略日起算历元的由来。 儒略日的引入提出了两个需要解决的问题,就是怎样在公历年月日序数和儒略日序数之间进行互换。下面先讨论如何由公历年月日序数 y,m,d 推算当天的儒略日序数 J(y,m,d)。这个问题稍许有些难度,我们不妨后退一步,先把需要考虑的时间范围限制在一年之内,考虑如何由月日序数 m,d 推算当天的积日 S(m,d),从中或许可以找到解决问题的思路。
4、积日
一年里某日积日概念的引入,想法与儒略日的设立是一致的。即忽略月的存在,从年初第一天起连续计日,直至年末的最后一天。只不过它的作用范围只限于一年之 内,到了下一年再重新开始计算,不像儒略周期那样长达7,980年。这样,年,特别是闰年的因素就不必考虑了。1月1日的积日为1,12月31日的积日平 年时为365,闰年时为366。 积日的计算原本并无难处,却由于2月的特殊情形而变得复杂了。凯撒把闰月设在2月,奥古斯都从2月减去一日,都缘于2月是古罗马处决死囚的凶月,颇具人道 意味。但由于闰月设在2月,2月的日数变得不唯一,必须按平年和闰年加以调节,这就影响到以后十个月份的积日也不唯一。设想一下,如果把闰月设在年末的最 后一天,那么除了这一天之外,一年中所有其他日子的积日数都是唯一不变的,包括每月0日的积日,问题就变得简单了。 但是,闰月在二月的设置尽管有不便之处,却不是随便可以更改的。我们最好换一个思路考虑,能否在计算的过程中,把3月当作当年的第一个月,而把次年的2月 当作当年的最后一个月?只要把次年的1月和2月称为当年的13月和14月,每年由3—14共十二个月组成,这样一来,闰月落在了一年的最后一个月,只要计 算完毕后再恢复原来年月的归属关系,困难就迎刃而解了。 这样安排之下,对于月序数为 m,日序数为 d 的一天,对应的积日可以表示成: S(m,d)=S0(m)+d (3<=m<=14, d<=31) (1) 式中S0(m)是这个月0日的积日,叫做月首积日。日序数为0的日子原本是不存在的,引入它只是为了计算的方便。由4月1日前推1天,就是3月31日,也 就是4月0日,3月0日应是上年14月的最后一天,依上年是否闰年而分别是14月29日和28日,其余依此类推。月首积日的计算有下面的公式可用: S0(m)=Floor(30.6*(m+1))-122 (3<=m<=14) (2) 公式中的符号Floor是实数向下取整运算,Floor(x)代表不大于括号中实数x的最大整数,数学书中常用[x]表示。我们这里由于接着要讲编程,故直接借用了程序语言的表达式。(2)式的正确性,读者可以自行验证。结合(1),(2)两式,就得到了计算积日的公式: S(m,d)=Floor(30.6*(m+1))+d-122 (3<=m<=14, d<=31) (3) 以上的思考过程带给我们如下的启发:类似于求取积日的做法,如果我们能求出当年3月0日的儒略日J0(y),那么当天的儒略日J(y,m,d)就可以表示成 J(y,m,d)=J0(y)+S(m,d) (4) 问题归结为如何求取年首儒略日。
5、GPS周
GPS周(时)是GPS系统内部所采用的时间系统。 时间零点定义的为:1980年1月5日与1980年1月6日晨之间的之夜。最大时间单位是周(一周:604800秒)。 表示方法:从1980年1月6日0时开始起算的周数加上周内时间的秒数(从每周周六/周日之夜开始起算的秒数,例如:1980年1月6日0时0分0秒的GPS周:第0周,第0秒 2004年5月1日10时5分15秒的GPS周:第1268周,第554715秒
6、UTC时间
UTC是协调世界时(Universal Time Coordinated)的英文缩写,是由国际无线电咨询委员会规定和推荐,并由国际时间局(BIH)负责保持的以秒为基础的时间标度。UTC相当于本初子午线(即经度0°)上的平均太阳时,过去曾用格林威治平均时(GMT)来表示。
那么UTC与世界各地的时间应如何换算呢?它是将全世界分为24个时区,地球的东、西经各180°(共360°)被24个时区平分,每个时区各占15°。以经度0°(即本初子午线)为基准,东经7°30′与西经7°30′之间的区域为零时区;东经和西经的7°30′与22°30′之间的区域分别为东一区和西一区;以此类推。从零时区起,向东每增加一个时区时间加1小时,向西每增加一个时区减1小时。UTC与零时区时间相同,以2004年7月15日0000UTC(即本初子午线上2004年7月15日零点整)为例,美国旧金山位于西八区,比零时区晚8小时,故此时旧金山时间为2004年7月14日16点整;而北京位于东八区,比零时区早8小时,此时北京时间为2004年7月15日8点整。
根据国际电信联盟(ITU)规定,为了统一起见,使用一个统一的时间,在国际无线电通信中除另有指明外,均应使用UTC,并用4位数字表示。
有点罗嗦了,说了这么多概念。接下来咱们就说一下程序怎么处理。还是围绕我同学的那个问题。公历时间转换为GPS周。根据上面的定义我想了两种方法:
1、根据GPS时间的定义之间计算输入日期距GPS起始时间的天数。然后就可以轻松的得到输入的时间是GPS的第几周的第几天了
2、先把输入时间换算为儒略日,然后根据儒略日与GPS的时间关系得到输入的时间是GPS的第几周的第几天
这两种方法的关键点是如何计算闰年,这个算法算好了剩余的就是小菜一碟了。
下面就介绍闰年的计算
每个可被 4 整除的年是一个闰年。 不过,可被 100 整除的年不是闰年。 但是,可以被 400 整除的年还是闰年。算法见GPS周算法程序的第一个函数
有的朋友会问第二种算法是不是太罗嗦了。其实之所以会有第二种算法,原因是在默写编程语言中提供了儒略日计算的函数。例如Delphi中就提供了DateTimeToJulianDate函数。这样就可以直接使用这个函数算出自然天了。这样比第一种算法要简单很多,但是在没有提供这种算法的语言中我建议使用第一种方法。
下面是我用VB写的一种处理方法。抛砖引玉。大家不要见笑
建立一个VB的工程,添加一个按钮控件和5个文本框。
形状如下:
复制一下代码到代码区编译即可!
''计算输入年份到GPS初始年之间的闰年的年数,因为每个闰年需要多加一天
Private Function leap(ye As Long) As Long
Dim k As Long
k = 0
Do While (ye >= 1980)
If (ye Mod 4 = 0 And ye Mod 100 <> 0 Or ye Mod 400 = 0) Then
k = k + 1
End If
ye = ye - 1
Loop
leap = k
End Function
Private Function tianshu(ByVal y As Long, ByVal m As Long, ByVal d As Long) As Long
Dim total As Long
total = 0
total = (y - 1980) * 365
total = total + d - 6
If m >= 3 Then
total = total + leap(y)
Else
If y Mod 4 = 0 And y Mod 100 <> 0 Or y Mod 400 = 0 Then
total = total + leap(y) - 1
Else
total = total + leap(y)
End If
End If
m = m - 1
Do While (m >= 1)
If (m = 1 Or m = 3 Or m = 5 Or m = 7 Or m = 8 Or m = 10) Then
total = total + 31
Else
If (m = 4 Or m = 6 Or m = 9 Or m = 11) Then
total = total + 30
Else
total = total + 28
End If
End If
m = m - 1
Loop
tianshu = total
End Function
Private Sub Command1_Click()
Dim year As Long
Dim manth As Long
Dim day As Long
Dim i As Long
i = 0
year = Val(Text1.Text)
i = 0
manth = Val(Text2.Text)
If manth <= 0 Or manth > 12 Then
i = MsgBox("输入的月份非法", 0 + vbExclamation, "提示")
End If
day = Val(Text3.Text)
If i <> 1 Then
i = 0
If (manth = 1 Or manth = 3 Or manth = 5 Or manth = 7 Or manth = 8 Or manth = 10 Or manth = 12) Then
If day <= 0 Or day > 31 Then
i = MsgBox("输入的日期非法", 0 + vbExclamation, "提示")
Text3.SetFocus
End If
Else
If (manth = 4 Or manth = 6 Or manth = 9 Or manth = 11) Then
If day <= 0 Or day > 30 Then
i = MsgBox("输入的日期非法", 0 + vbExclamation, "提示")
Text3.SetFocus
End If
Else
If (year Mod 4 = 0 And year Mod 100 <> 0 Or year Mod 400 = 0) Then
If day <= 0 Or day > 29 Then
i = MsgBox("输入的日期非法", 0 + vbExclamation, "提示")
Text3.SetFocus
End If
Else
If day <= 0 Or day > 28 Then
i = MsgBox("输入的日期非法", 0 + vbExclamation, "提示")
Text3.SetFocus
End If
End If
End If
End If
Else
Text2.SetFocus
End If
If i <> 1 Then
Text4.Text = tianshu(year, manth, day) \ 7
Text5.Text = tianshu(year, manth, day) Mod 7
End If
End Sub
|